2011-12-02

統計量の計算とp値への変換

R カイ二乗検定 z検定

こちらでchisq.test()とprop.test()の話がある
それぞれの関数の中身を見るには

prop.test
chisq.test

とすればよい
統計量とp値計算がどうなっているかを確認したいとする
どちらの関数も、いろんな場合に対応しようとして、場合分けがたくさんあるので、面倒くさいが…
関数の値の返却の部分を確認することにする
- prop.testの場合

RVAL <- list(statistic = STATISTIC, parameter = PARAMETER, 
        p.value = as.numeric(PVAL), estimate = ESTIMATE, null.value = NVAL, 
        conf.int = CINT, alternative = alternative, method = METHOD, 
        data.name = DNAME)
    class(RVAL) <- "htest"
    return(RVAL)

- chisq.testの場合

structure(list(statistic = STATISTIC, parameter = PARAMETER, 
        p.value = PVAL, method = METHOD, data.name = DNAME, observed = x, 
        expected = E, residuals = (x - E)/sqrt(E), stdres = (x - 
            E)/sqrt(V)), class = "htest")

どちらの関数も、関数の中でSTATISTIC,PARAMETER,PVALを最終的に返しているので、それの出方をソースの上流にたどるのがよい
- prop.testの場合

STATISTIC <- sum((abs(x - E) - YATES)^2/E)
PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
        
z <- sign(ESTIMATE - p) * sqrt(STATISTIC)
        else z <- sign(DELTA) * sqrt(STATISTIC)
PVAL <- pnorm(z, lower.tail = (alternative == "less"))

- - と、なっていて、STATISTICを出してから、pchisq()を使ってp値化するか、STATISTICをzに変換してからpnorm()でp値化するかに場合分けしていることが見える
  - pchisq()で扱われるカイ二乗統計量とpnorm()で扱われる正規分布の統計量とは、sqrt(STATISTIC)とあるように、平方根←→二乗の関係であることも見える。また、正規分布は、２方向に裾を引くのでsign(...)と１次元の向きを気にしてあることも見える。カイ二乗分布は0以上と、片側にしか裾がない(正規分布を二つ折りにした(様なもの))ので、sign()でのような「向き」は気にできない(正規分布は、向きありの『距離』を測っているのに対して、カイ二乗分布は単純に『距離』を測っているだけ…)
    - YATESがイェーツの補正項。これに値を入れるかどうかがさらに上流で選択されていることも見える
- chisq.testの場合

PVAL <- pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
PVAL <- (1 + sum(ss >= almost.1 * STATISTIC))/(B + 1)

- …のようにSTATISTICからpchisq()を使ってPVALを計算している場合もあり、B(乱数使用・シミュレーションベースの計算のための変数)を使っているように、それぞれの違いがある

STATISTIC <- sum((abs(x - E) - YATES)^2/E)
STATISTIC <- sum((x - E)^2/E)

- 統計量についても２式ある。YATESの補正項を使う場合と使わない場合に分けてあることもわかる。YATES項を使う場合はさらに、その項に0を与える場合(結果として補正しない)と非0を与える場合とに分けてある。それぞれ理由がある(場合分け条件でわかる)
結局、prop.test()もchisq.test()も基本的には、観察値のセットと期待値のセットの差(向きなし〜abs()を使っているから)をセルごとの『距離』として、その二乗を足し合わせたい(これはユークリッド距離の定義)けれど、表のセルごとに重みづけを変えたいので、二乗を対応するセルの期待値で割った上で足し合わせている
- 期待値で割らないで足し合わせるとすると、 $\sum_i r_i^2$ を計算していることになって、これは、「ただの距離(の二乗」。この値が等しい２つの観察があったとすれば、それは、空間では同一半径の球面上の２点となる。
- 期待値で割ると、「方向によって、距離の重みが変わる」ので、カイ二乗値が等しい２つの観察表は、空間において、ある楕円(楕球)上の２点となる。→マハラノビス距離のWiki

STATISTIC <- sum((abs(x - E) - YATES)^2/E)
# もしくは
STATISTIC <- sum((x - E)^2/E)

さて。統計量という『距離』のようなものを、p値化しています。上ではpchisq()とかpnorm()とか、既知の分布関数を使ってp値にしています。
分布が知られていないときには、分布関数が役に立たなくなりますが、『距離』の方は役に立つので、分布関数が役に立たないような検定(SNPの２ｘ３表の遺伝形式検定(優性・劣性・相加その他)とかでは、距離からp値に変換するところで工夫をしたりすることになります。マルチプルテスティングも共通(２ｘ３表に遺伝形式モデルを持ち込んで検定することもマルチプルテスティングの一つ)。
「向き」
- 正規分布を考えたら、前後の２方向があったところが、カイ二乗分布にしたら、１方向になりました
- ２方向あるときに、「両側」とか「片側」とかの違いが登場しました
- 「方向」だけで、世界が分けられるのは、１次元空間だけです
- ２次元空間、多次元空間では、方向を考慮した位置を考えると、方向の数は無限にあります。
- したがって、「両側」とか「片側」とか、の２分はできなくなるわけですが、では、「方向」を考慮しなくなるのかと言えば、そんなことはやはりないわけです
- １次元のときに、「両側」「片側」と言ったのは、言い換えると、「全方向どちらに向かってもよい」「特定の方向への距離を考える」となります
- 後者は２以上次元でも使える表現です。「全方向」「特定の方向」。ただし「特定の方向」として、選べる方向の数が無限にあるので、「特定の方向を考慮しての検定は、無限に実施できる」と言うことになります
- さらに言えば、次元が上がると「方向」は向きが増えるだけでなく、「種類」も増やせます
- 「方向」は矢印が示すことで、矢印には２次元空間なら、垂線が、３次元空間なら垂直面が対応できます。このように、「方向」を気にするということは、２次元空間なら線を、３次元空間なら面を意識することに対応します
- じゃあ、２次元空間を、矢印と線とに考えること、３次元空間を、矢印と面とで考えること・・・というように扱っていて、１次元のときに「片側」で考えているのは、１次元空間を矢印と点とで考えていることになります
- n次元を矢印(１次元)とその他(n-1次元)とで考えることをしている、と、書き換えられるわけですが、そうしたら、n次元をk次元とn-k次元とに分けて考えていけないわけはなくなりますから、そうすることも可能になります
- マルチプルテスティングは、そういう意味合いもあります
さらに。ここまでは、矢印とそれに垂直なものとか、「まっすぐなもの」と「かくかくしたもの」で考えていたわけですが、それは「線形」なわけで、じゃあ、「線形ではない〜非線形」というわけで、こちらが向かっている先につながるわけでしょう
距離に関して言えば、DNA配列の異同を測る距離としては、５か所い点変異が入っていたときの距離は、５次元データだけれど、にするわけではなく、５にするなど、距離の取り方も変化させることがあるので：
- データ
- 「何かしらの決まり」で「距離」にする
- 「距離」をp値にする
…３段構えで、それぞれのところに決まりを作る必要があるわけですが、上の、prop.test(),chisq.test()の話しでは、「距離」にするところは、同じだったことになります