sigma cutoff の濫用は、バイアスの原因となる。負の観測を切り捨てることで測定値を正にバイアスし、ノイズしかない場合でもシグナルがあるかのような統計値を生み出してしまう。このため、業界のコンセンサスとしては、-3 よりも大きくするのは良くないとなっている。
XDS では -3 が使用され FAQ - XDSwiki、HKL2000 付属の SCALEPACK でもそうだ(Scalepack Keywords | HKL Research Inc. の SIGMA CUTOFF の項を参照)。CrystFEL の process_hkl では、そもそも使用すべきでないとされている(CrystFEL best practices)。
この背景については、上で示した SCALEPACK のマニュアルに詳しい解説があるが、実際に R で確認してみた。
BG <- 5 # 背景ノイズの平均値 N <- 1000 # 観測の数 # あるシェルで有意なシグナルが全くないなら、どの反射を測定しても # 背景ノイズだけの値になるはず counts <- rpois(N, BG) counts_minus_BG <- counts - BG plot(density(counts_minus_BG), ylim=c(0, 0.5), xlab="Count", main="Distribution of background subtracted I") for (sigma_cutoff in c(-Inf, -4, -3, -2, -1, -0.5, 0)) { filtered <- counts_minus_BG[counts_minus_BG > mean(counts_minus_BG) + sigma_cutoff * sd(counts_minus_BG)] filtered_I_over_sigma <- mean(filtered) / sd(filtered) lines(density(filtered)) cat("sigma cutoff =", sigma_cutoff, "I/sigma =", filtered_I_over_sigma, "\n") }
結果は以下の通り。理論上は、I/sigma = 0.0 になるはずなのに、sigma cutoff を大きくするにつれて、どんどん値が増えてしまっている。
sigma cutoff = -Inf I/sigma = 0.00994426 sigma cutoff = -4 I/sigma = 0.00994426 sigma cutoff = -3 I/sigma = 0.00994426 sigma cutoff = -2 I/sigma = 0.02157932 sigma cutoff = -1 I/sigma = 0.2341382 sigma cutoff = -0.5 I/sigma = 0.5603886 sigma cutoff = 0 I/sigma = 1.59822
cutoff 適用後の反射強度の分布は次のようになる。離散的な値を density でプロットしてしまったので、変なコブができている。本来は hist を使うべき。
次に、CC1/2 について検討してみた。
BG <- 50 # 背景ノイズの平均値 Nref <- 500 # シェルにある反射の数 MULTI <- 20 # 多重度 N <- Nref * MULTI # あるシェルで有意なシグナルが全くないなら、どの反射を測定しても # 背景ノイズだけの値になるはず counts <- rpois(N, BG) counts_minus_BG <- counts - BG for (sigma_cutoff in c(-Inf, -4, -3, -2, -1, -0.5, 0)) { filtered <- counts_minus_BG filtered[counts_minus_BG < mean(counts_minus_BG) + sigma_cutoff * sd(counts_minus_BG)] <- NA # 反射 1 の観測 1, 2, 3, ..., MULTI/2 # 反射 2 の観測 1, 2, 3, ..., MULTI/2 # ... # 反射 Nref の観測 1, 2, 3, ..., 2/MULTI # 反射 1 の残りの観測 # 反射 2 の残りの観測 # ... # 反射 Nref の残りの観測 # という順で記録されているイメージ obs <- matrix(filtered, ncol=MULTI / 2) mergedI <- rowMeans(obs, na.rm=TRUE) cc_half <- cor(mergedI[1:Nref], mergedI[(Nref + 1): (Nref * 2)], use="complete.obs") cat("sigma cutoff =", sigma_cutoff, "CC1/2 =", cc_half, "\n") }
結果を見ると、sigma cutoff の影響を受けにくいようである。若干改善しているが単調増加ではないし、こんな統計値で有意なシグナルがあるとは誤認しないだろう。
sigma cutoff = -Inf CC1/2 = -0.01819305 sigma cutoff = -4 CC1/2 = -0.01819305 sigma cutoff = -3 CC1/2 = -0.01821197 sigma cutoff = -2 CC1/2 = 0.002451837 sigma cutoff = -1 CC1/2 = 0.04601235 sigma cutoff = -0.5 CC1/2 = 0.02382984 sigma cutoff = 0 CC1/2 = 0.02476777
閾値が -1 付近で CC が最大になるという挙動を理論的に説明できるのか気になっている。なお、パラメータを変えて色々試してみたが、同じような感じになった(←いい加減)。
