廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

「研究仮説が正しい確率」(豊田、2017)のシミュレーション

今日の川柳

「研究仮説が正しい確率」(豊田, 2017)の分布 - 廿TT ではやや甘の結論だったので、このエントリでは「研究仮説が正しい確率」という言葉は使うべきできない、とはっきり言い切ることを目指します。

「p 値を使って学術論文を書くのは止めよう」(https://www.jstage.jst.go.jp/article/sjpr/60/4/60_379/_pdf)を読むと、文脈上、豊田氏の主張は 「事前分布として幅の広い一様分布を選んで、事後分布からサンプリングを行い、事後分布からのサンプルにおいて『研究仮説』が正しければ1、さもなくば0とおく。その期待値をとれば研究仮説が正しい確率になっている」 と見受けられます。

「『研究仮説が正しい確率が 97.6%であり,正しくない確率は 2.4% である』という解釈は,きわめて自然であり誤解の生じようがない。」(https://www.jstage.jst.go.jp/article/sjpr/60/4/60_379/_pdf)とされていますが、具体的な問題を考えると、「研究仮説が正しい確率」は解釈のしかたがかなり難しい数字になってしまいます。

仮にデータ x=8 が試行回数パラメータ n=10、成功確率パラメータ p=0.49 の二項分布から生成されたことを知っているとします。尤度をデータを生成した確率分布と同じく二項分布として、事前分布に区間 0 から 1 の一様分布を選んだ場合、事後分布はパラメータ a=x+1、b=n-x+1 のベータ分布になります。

この場合に p>0.5 という「研究仮説が正しい確率」を計算すると、約0.97になります。

R による計算例はこちら。

set.seed(4444)
n <- 10
p <- 0.5
x <- rbinom(1,n,p)
print(x)
psmp <- rbeta(5000,x+1,n-x+1)
print(mean(psmp>0.5))
#[1] 0.968

ちなみに研究仮説が p>0.5 という単純な形なので、サンプリングをしなくても確率を計算できます。

print(pbeta(0.5,x+1,n-x+1,lower.tail = FALSE))
#[1] 0.9672852

0.97という一見高い確率で研究仮説が支持されたとしても、この仮説が間違っていることをシミュレーションを設定したぼくらはすでに知っています。

「『研究仮説が正しい確率が 97.6%であり,正しくない確率は 2.4% である』という解釈は,きわめて自然であり誤解の生じようがない。」と、豊田先生は述べていますが、誤解は生じるでしょう。

97%の確率で研究仮説が支持される、と言われたら、実運用上は研究仮説に統計的なお墨付きを得たと判断してしまう場合が多いと思います。

もうちょっとシミュレーションを続けます。

試行回数が10というのはちょっと少ない気がするので、1000に増やし、シミュレーションも10000回繰り返します。

set.seed(1)
g <- numeric(10000)
n <- 1000
for(i in 1:10000){
  x <- rbinom(1,n,p)
  ptilde <- rbeta(5000,x+1,n-x+1)
  g[i] <- mean(ptilde>0.5)
}
hist(g, main="probability of p > 0.5")

これが研究仮説が正しい確率の分布です。

f:id:abrahamcow:20200131064446p:plain

低い値が出ることがやや多めですが、高い確率が出ることも十分あります。

そうなってくると知りたいのは「研究仮説が正しい確率」がどういう性質を満たしていれば、「研究仮説が正しい確率」と呼べるのか、ということです。

それについて説明している人はぼくの知る限り一人もいません。

そのため、「研究仮説が正しい確率」を使うとしたら、こう計算したらとにかく研究仮説が正しい確率なんだと思い込む以外にありません。

それはデータをもとに仮説を検証するという本来の目的から、かけ離れた行為になってしまいます。

「研究仮説が正しい確率」という言葉を使わないことを強くおすすめします。

おまけ。もう一つシミュレーション。データを生成した二項分布のpが0.5のとき、p>0.5という「研究仮説が正しい確率」の分布。

f:id:abrahamcow:20200131070018p:plain

p=0.5 なので p> 0.5という仮説は間違っていますが、ほぼ一様分布。

(今の段階で「研究仮説が正しい確率」に解釈を与えるとしたら、「サンプルサイズが十分大きくてモデルが単純なとき、p値と似たような値を与える指標」ということになりそうな気がしますが、これについてはまだ自信がない。)