LoginSignup
1
0

More than 1 year has passed since last update.

@kenkenvw

p-hacking問題をシミュレーションしてみる

Last updated at Posted at 2023-04-01

１．はじめに

　twitterで、「p値を下げるために後からサンプルサイズを増やすな」って正直意味がわからないっていう内容のツィートを見かけました。
　確かに最終のサンプルサイズによるp値は同じだし、データも全く同じものを使っているのにp-hackingだ！不正であると言われても今ひとつ、腑に落ちないところがあります。そこで、ちょっとシミュレーションしてみたいと思います。

２．DATAの設定

帰無仮説として２グループの平均値のt検定を試みます。

同じ平均と分散を持つグループからのサンプリングのため、p値は一様分布し、有意差を計算するためにα水準と同じ割合にならないといけません。

（α水準５％なら、同じグループであっても２０回に１回は有意差がでる割合になる）

試しに、同じグループからのp値の分布を見てみます。平均は50、標準偏差５に従います。サンプルサイズは100です。10000回検定を繰り返します。

rm(list = ls())
set.seed(12345)
record = NULL

for (i in 1:10000){
group1 = rnorm(100,50,5)
group2 = rnorm(100,50,5)

record = append(record,t.test(group1,group2)$p.value[1])　#record変数にp値を記録する

}

hist(record)

image.png
ほぼ理論どおり、一様分布しました。5％有意水準なら、ほぼ棄却割合は５％になります。

sum(record <= 0.05) / 10000

>[1] 0.0528

有意差のでる割合はα水準どおりほぼ、5%となりました。

３．p-hackingのシミュレーション

p-hacckingのグループは計画のサンプルサイズは100ですが、サンプルサイズが50になった段階で一度、検定をします。ここではα水準を5%とします。そして有意差がなければ、計画通りの検定をします。一度、検定をするかどうかだけの違いとなります。両者ともにデータは全く同じものです。

set.seed(12345)

record1=NULL
record2=NULL

for( i in 1:10000){
group1 = rnorm(100,50,5)
group2 = rnorm(100,50,5)

#通常のグループ。サンプルサイズ100で検定（計画どおり）。
if(t.test(group1,group2)$p.value[1] <= 0.05) {
record1 = append(record1,1) 
}else{record1=append(record1,0)}


#p-hacikibのグループ。サンプルサイズ50になった段階で一度検定する。
hack1=group1[1:50]
hack2=group2[1:50]

#p-hacikibのグループは、有意差がなければ、計画どおり検定を実施する
if(t.test(hack1,hack2)$p.value[1] <= 0.05) {
record2 = append(record2,1) 
}else{
  record2 = append(record2,0)
  hack1=group1[1:100]
  hack2=group2[1:100]
 
    if(t.test(hack1,hack2)$p.value[1] <=0.05) {
    record2 = append(record2,1)
    }else{record2 = append(record2,0)}
}
}

有意差がでた回数を比較してみましょう。10000回の試行なので、500回程度が有意差がでる期待値となります。

sum(record1)  #通常の検定グループ
sum(record2)　#p-hackingグループ

>[1] 528
>[1] 857

通常のグループは、ほぼ期待値どおりの回数だったのに対し、p-hackingグループは期待値よりも大きくなり有意差がでやすくなっています。p-hackingグループは、同じデータで途中で検定を1回実施しただけですが、有意差がでやすくなっていることがわかります。どうしてこうなったのでしょうか？
p-hackingグループは試行回数は同じですが、検定の回数が多くなっていることが原因のようです。

length(record1)　#通常のグループ
length(record2)　#p-hackingグループ

>[1] 10000
>[1] 19486

では、今度はp-hackingグループは、もっと有意差を出やすくするために、同じデータを用いて20回ずつ小出しに検定を実施し、検定回数をもっと増やしてみます。

set.seed(12345)
record1=NULL
record2=NULL
end=0
for( i in 1:10000){
group1 = rnorm(100,50,5)
group2 = rnorm(100,50,5)

#通常のグループはサンプルサイズ100で検定（計画どおり）。
if(t.test(group1,group2)$p.value[1] <= 0.05) {
record1 = append(record1,1) 
}else{record1=append(record1,0)}

#p-hackingグループはサンプルサイズ20ずつ追加し、検定を繰り返す。
hack1=group1[1:20]
hack2=group2[1:20]

if(t.test(hack1,hack2)$p.value[1] <= 0.05) {
  record2 = append(record2,1) 
  }else{
  record2 = append(record2,0)
  hack1=group1[1:40]
  hack2=group2[1:40]
  
  if(t.test(hack1,hack2)$p.value[1] <=0.05) {
  record2 = append(record2,1)
    }else{record2 = append(record2,0)
     hack1=group1[1:60]
     hack2=group2[1:60]
     if(t.test(hack1,hack2)$p.value[1] <=0.05) {
     record2 = append(record2,1)
     }else{record2 = append(record2,0)
        hack1=group1[1:80]
        hack2=group2[1:80]
        if(t.test(hack1,hack2)$p.value[1] <=0.05) {
        record2 = append(record2,1)
        }else{record2 = append(record2,0)
        hack1=group1[1:100]
        hack2=group2[1:100]
        if(t.test(hack1,hack2)$p.value[1] <=0.05) {
        record2 = append(record2,1)
        }else{record2 = append(record2,0)
        }
      } 
     }
    }
  }

}

それでは、有意差がでた回数を比較してみましょう。

sum(record1)　#通常のグループ
sum(record2)　#p-hackingグループ

>[1] 528
>[1] 1441

p-hackingグループは、ねらいどおり同じデータを用いて、より有意差をでやすくすることができました。
このように検定回数を増やせば増やすほどに、有意差がでやすくなります。

今度はt検定で有意差がでなければ、ノンパラであるwilcox検定をしてみるという、検定方法をいろいろと試し、有意差が出た検定法だけを選ぶというp-hackingを試みます。

set.seed(12345)

record1=NULL
record2=NULL

for( i in 1:10000){
group1 = rnorm(100,50,5)
group2 = rnorm(100,50,5)

#通常のグループ。サンプルサイズ100で検定（計画どおり）。
if(t.test(group1,group2)$p.value[1] <= 0.05) {
record1 = append(record1,1) 
}else{record1=append(record1,0)}


#p-hacikibのグループは、t検定またはwilcox検定を実施し、どちらか有意差がでた場合だけを選ぶ
if(t.test(group1,group2)$p.value[1] <= 0.05 | wilcox.test(group1,group2)$p.value[1] <= 0.05) {
record2 = append(record2,1) 
}else{record2 = append(record2,0)}
  
}

それでは、比較してみましょう。

sum(record1)　#通常のグループ
sum(record2)　#p-hackingグループ

>[1] 528
>[1] 624

p-hackingグループの方が若干ですが有意差が出やすくなりました。

最後に、p-hackingグループは両者の方法を組み合わせてみます。

set.seed(12345)

record1=NULL
record2=NULL

for( i in 1:10000){
group1 = rnorm(100,50,5)
group2 = rnorm(100,50,5)

#通常のグループ。サンプルサイズ100で検定（計画どおり）。
if(t.test(group1,group2)$p.value[1] <= 0.05) {
record1 = append(record1,1) 
}else{record1=append(record1,0)}


#p-hacikibのグループ。サンプルサイズ50になった段階で一度検定する。
hack1=group1[1:50]
hack2=group2[1:50]

#p-hacikibのグループは、2通りの検定を実施し、有意差がなければ、計画どおり検定を実施する
if(t.test(hack1,hack2)$p.value[1] <= 0.05 | wilcox.test(group1,group2)$p.value[1] <= 0.05) {
record2 = append(record2,1) 
}else{
  record2 = append(record2,0)
  hack1=group1[1:100]
  hack2=group2[1:100]
 
    if(t.test(hack1,hack2)$p.value[1] <=0.05) {
    record2 = append(record2,1)
    }else{record2 = append(record2,0)}
}
}

sum(record1)　#通常のグループ
sum(record2)　#p-hackingグループ

[1] 528
[1] 939

p-hackingグループは、途中検定の追加だけのときが 857だったのに対し、さらに有意差の出やすさを向上させることができました。
これで、検定回数も1回だけとし、検定方法も最初からどちらかの方法だけで実施したと報告すれば、誰にもバレないで有意差を出しやすくすることができます。さらに検定方法も、いろいろな方法を組み合わせれば組み合わせるほど有意差が出やすくなることが期待されます。

このようにp-hackingの実態は、検定の繰り返しです。特に有意差が出ない場合は、つい違う検定方法を試してみたりしますが、これもりっぱなp-hackingになりますね。

4.終わりに

p-hacking問題は、意外とそのことを知らないでやってしまっている人が多いような気がします。
検定をするときは、気をつけましょう。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0