結晶学的 R factor と統計学におけるR^2

先日のラボ内勉強会で、決定係数 R^2 値(wikipedia)の話が出た。R^2 の定義は
$R^2 = 1 - {\sum_i(y_i - f_i)^2 \over \sum_i(y_i-\bar{y})^2}$
である。すなわち、モデルとデータの二乗残差をデータの分散で割って、1から引いたものである。

結晶学におけるモデルとデータの適合度評価には、R factor(wikipedia) を用いる。式は、
$R = {\sum_{hkl}||F_{calc}| - |F_{obs}|| \over \sum_{hkl}|F_{obs}|$
である。

完璧にモデルとデータがフィットする場合、R-factor は 0.0 になるが、R^2 は 1.0 になる。とはいえ、これは本質的な差ではない。それよりも気になるのは、なぜ二乗残差を使わずに絶対値和になっているのかということと、分母の違いである。

＃違いについてもっと考察

絶対値があると数学的に、あるいは最適化の際に扱いにくくないのか、という質問を受けたことがある。構造因子 F は複素数なので、内側や分母の絶対値は本質的ではない。外側が問題である。結晶学では R-factor が精密化における最適化目的関数ではないから、数学的に扱いやすいかどうかはそれほど問題にならない。実際、目的関数としては二乗残差やRice分布に基づく尤度が用いられている。

結晶学では R-factor の p値を求めたり、検定することはない。そのため、R-factor の理論分布を得るのが困難であることは問題になりにくい。もっとも、ランダムモデルの R-factor がいくつになるかは理論的に求められている(twin なしで 60% くらいだったっけ?)。また、Protein Data Bank にある約8万の結晶構造から、分解能ごとの R-factor の経験的分布は知られている。Phenix では "Compare statistics" ボタンを押すことで、自分のモデルの統計値が他の構造と比べてどうなのか確認することができる。直感的には、R work が分解能の10倍程度(2.0Å の構造なら 20%)、R free と R work の差が2-3%というのが目安で、分解能 2.0Å なのに R factor が 30% を切りません、となると頭の中で黄信号が灯るわけである。

＃統計学におけるモデル比較の問題と、結晶学におけるモデル比較の問題について考える
＃ F statistic, BIC, AIC, etc...
＃化学的常識により束縛をかけて、自由度を減ずること。その場合に、自由度の減少をどう定量すべきか
＃束縛をかけることは、ベイズ的には事前確率を導入すること??

低分子のX線結晶学では、R factor は 5% を切るようなこともザラで、ほとんど R sym に漸近するという。つまり、モデルが測定誤差の範囲内で完全に観測値を説明できるということになる。一方、生体高分子の結晶学では、R factor は 20% をかろうじて切ったあたりが普通である。この違いがどこに起因するのかは明らかではない。James Holton 氏が ccp4bb に書いていたところによると、氏の回折シミュレーションプログラム MLFSOM に実装されているノイズ源(photon counting statistic, シャッタのジッタ, etc) 以外にありそうだということである。ACA 2013 の Otwinowski氏らの発表の abstract を見ると、異方性振動・溶媒領域・表面側鎖のモデリングが課題とのことである。

これも、数週間前に書いた原稿を掘り起こしてきて、昨日の勉強会の内容を付け足したので、支離滅裂な内容になってしまった。