確率変数としてのp-valueの感覚をRで身につける

R statistics

R Advent Calendar 2012の7日目です．もう一週間経ったのですね．


The American Statistician, August 2008, Vol. 62, No. 3 P-Values are Random Variables (PDF) で「学生がp-valueについて理解するためにはp-valueが確率変数であることを強調した方がいいんでない」という提案がありました．
そこで，何の脈絡もありませんが，p-valueは確率変数であり，帰無分布からのp-valueは一様分布に従うことをRで簡単に確かめたいと思います．

帰無分布からのp-valueは一様分布に従う

まず，ある統計検定量Tが帰無仮説H0の下で確率密度関数 f(T) に従うとします．
p-valueの定義は，帰無仮説H0の下である検定統計量Tがt以上の値を取る確率です．
p=P(T \gt = t|H_{0}) (1)
p-valueは以下のようにも書けます．
p=1-P(T \lt t|H_{0}) (2)
ここで，累積分布関数 F(T) を考えると，次のように書けます．
p=1-F(t) (3)
F(T)は単調増加なので，
P(T \lt t|H_{0})=P(F(T) \lt F(t)) (4)
(2), (3), (4) より
P(F(T) \lt F(t)) = F(t) (5)
(5) より，F(t)は一様分布となり，1-F(t)も一様分布となるので，p-valueも一様分布となります．

Rで実際に確かめる

（確率に関連する関数について詳しくはこちらをご覧ください scratch-R: probability plots）


まず，ある統計量が平均1，分散1の正規分布に従うとします（帰無仮説）．そこから10000個のデータをサンプリングします．これはある統計検定（棄却されなそうな）を10000回やっているものと考えて下さい．

m <- 1  #平均
v <- 1  #分散
n <- 10000  #データ数
data <- rnorm(n, m, v)  #乱数生成


次に，この統計量が平均1，分散1の正規分布に従うときのp-valueの分布を描くと，一様分布になります．

pvalues <- 1 - pnorm(data, m, v)
hist(pvalues)

では，「平均0.5，分散1の正規分布」もしくは「平均2，分散1の正規分布」に従うときのp-valueの分布はどうなるでしょうか．

m.fake1 <- 0.5
pvalues.fake1 <- 1 - pnorm(data, m.fake1, v)
hist(pvalues.fake1)


m.fake2 <- 2
pvalues.fake2 <- 1 - pnorm(data, m.fake2, v)
hist(pvalues.fake2)

偏った分布になります．これは，さまざまな対立仮説ではp-valueの値のとりやすさが変わるということを直感的に表しています．

参考