ベイズマックスへようこそ。
統計に関するあれやこれやを自由気ままに投稿します。
全記事一覧はこちらから。

2020年12月10日

p値で有意と言えない効果もベイズなら効果があると言える?──事後分布に基づく仮説評価について──

 この記事は,Open and Reproducible Science Advent Calendar 2020およびベイズ塾 Advent Calendar 2020の10日目の記事です。この記事では,ベイズ統計学と再現性に関わる次のようなケースについて考えてみます。

 研究者のAさんは,ある実験操作をしたときとしなかったときで変数Vの平均値に差が生じるかどうかを調べるために,60名の実験参加者を30名ずつ無作為に統制群と実験群に割り当てて実験データを取得しました。変数Vの平均値を計算してみたところ,統制群は-0.028,実験群は0.348であったので,実験操作によって変数Vの平均値が数値的には増加しているように見えます(図1)。

dataplot.png
図1. Aさんが収集したデータ。黒い横線は各群の平均値。

 ところが,この変数VについてWelchのt検定 (両側検定,α = .05) を行ったところp = .163という結果となり,有意差は認められませんでした。諦めきれなかったAさんはふと,ベイズで分析し直せば実験操作の効果があると主張できるのではないかと思い立ち,Welchのt検定に相当するモデルを使って平均値差μdiff (= μ実験群 - μ統制群) の事後分布を推定しました(図2)。

posterior.png
図2. 2群の平均値差 (μdiff) の事後分布。

 この事後分布を利用して,μdiffが0よりも大きい確率 (図2の赤い領域の面積) を計算したところ,p(μdiff > 0) = 91%という高い値が得られました。この結果をもってAさんは,「この実験操作は変数Vを高める効果があるということが高い確信を持って言える」と結論づけました。

 この例は筆者による創作でありフィクションですが,これに類する報告を学会や研究会,論文等でときどき見聞きします。このような推論には問題があるのですが,無自覚にこれを行ってしまっていると思われる実例としばしば遭遇してしまうので,この状況はまずいと思い,この推論の何が問題なのかについて自分なりに整理してまとめようと思いました。あくまでも私見なので,鵜呑みにするのではなく参考意見として受け止めて頂ければと思います。ちなみに,Aさんの実験データ (架空) の生成と分析を行うR・Stanコードはこの記事の末尾に掲載しています。

1. ベイズの枠組みにおける仮説評価の方法

 いわゆるベイズ的な仮説評価には複数の方法がありますが,研究者によって推奨する方法が異なっており,議論もたくさん行われています(この記事ではそれぞれの長短について深入りはしません)。仮説評価の方法は大雑把に,モデル比較に基づく方法とパラメータ推定に基づく方法の2種類に分けることができます。モデル比較に基づく方法にはベイズファクター (周辺尤度・自由エネルギー) やWAICによるモデル比較が代表的です。先ほどの例に出てきた事後分布を使った仮説の評価はパラメータ推定に基づく方法に含まれます (他にもROPEを使った方法や95%確信区間が0をまたぐかどうかを基準にする方法などがあります;ついでに,ROPEを利用したベイズファクターの計算方法も存在します)。Makowski et al. (2019) は,事後分布においてパラメータが中央値と同じ符号になる確率 (つまり,0を下回る確率と0を上回る確率のうち大きいほう) のことをprobability of direction (pd) と呼んでいますが,先ほどの例はこのpdを使った仮説評価を行っていることと等しいです。pdは[.5, 1] の範囲を取り得る指標で,値が大きいほど効果の存在 (effect existence) の強い証拠とみなすことができます。一方,値が小さかったとしても効果がないことの証拠とはみなせないという非対称性があります。Makowski et al. (2019) も言及しているように,(少なくとも単純なモデルにおいては) pdはp値と一対一に対応します (効果の事前分布に(-∞, +∞) の範囲の非正則な一様分布を仮定したとき,pd = 1 – p/2が成り立つ;ただし,この記事の分析では各群の平均の事前分布に一様分布を仮定しているためμdiffの事前分布は三角分布となり,この式とは一致しない)。よくよく考えてみたら,上述したpdの性質はp値の性質と何ら変わりませんよね。帰無仮説が正しい時 (真の効果の大きさが0のとき) にはサンプルの変動によってp値は (0, 1)の範囲の一様分布に従いますが,pdも同様に(.5, 1)の範囲の一様分布に従います (事前分布にもよりますが)。そう考えると,p = .163で有意じゃないのにpd = .91で効果があるように見えるのは不思議な気がしませんか?

2. 何が問題なのか

 Aさんのpdの使い方には,少なくとも2つの問題点があるように思われます(もっとあるかもしれませんが,ここでは2つに絞ります)。1つ目は,pdが示す確率の値を素朴に解釈してしまっている点でしょう。例えば,天気予報で降水確率が80%であることを知れば,多くの人は傘を持って出かけるでしょう。降水確率が90%なら尚のこと,間違いなく雨が降るだろうと確信することでしょう。僕だったら降水確率が50%でも傘を持って出かけます(何なら毎日折り畳み傘を持参していますけど)。これと同じように,pd = .91を根拠に「実験操作の効果がある」と判断するのは妥当でしょうか。例えば,100回同じ実験をしたら91回はμdiff > 0となる,といった具合に, pdが実験の再現率を表すと解釈できるのであればこの判断は正当化できるかもしれませんが,残念ながらそうではありません (実際,今回の設定においては100回同じ実験を繰り返した時にpd > .90となるのはたった20回だけです)。そもそも,pdは必ず.5以上の値を取りますしね。パラメータの事後分布とそこから計算された指標は,あくまでも今回のサンプルで条件付けられた結果に過ぎません。言い換えれば,別の実験参加者60人に対して同じ実験を行った場合には異なる事後分布が得られるということです (あるいは,今回の実験参加者のうち10名がたまたま何らかの事情で実験に参加できず,たまたま都合の付いた他の10名が実験に参加していた平行世界があったとしたら,異なる事後分布が得られていたでしょう)。たまたまサンプルが偏ってしまった場合には,真の効果の大きさが0であっても大きなpdが得られてしまうことは十分あり得ます (偶然p < .05になり得るのと同じように)。このように,pdはそもそも直感的に解釈するのが難しい指標であるように思われます (少なくとも僕にとっては)。少なくとも降水確率のように素朴に解釈してしまうのはまずい,というのが1点目の問題でした。

 2点目の問題は,基準を明確にせずに効果の有無について二値判断をしてしまっている点です(i.e., 「この実験操作は変数Vを高める」)。今回はpd = .91だったので素朴に確率を解釈すれば効果がありそうな感じがしますが,pd = .84だったらどうでしょうか。じゃあpd = .73だったら?二値判断をせずに,pdの値をそのまま量的に解釈するのであればこの点は問題にならないのですが,前述した通りpdはそれほど直感的に解釈しやすい指標とは言えず,量的な解釈は意外と難しいように思います。それに,効果の存在を示すことはなんやかんやで学術的にも社会的にも求められることが多く,二値判断せずに論文を書くことは現状かなり難しいと思います。たとえ留保付きだとしても,効果の有無についての言及はほとんどの場合避けられないし,意図せずとも暗黙に二値判断をしてしまうこともあるし,論文で明示しなかったとしても読者が勝手に効果の有無を判断してしまうこともあり得えます (余談ですが,ベイズファクターによるモデル比較の場合にも,「モデルXがモデルYよりも支持された」といった具合に離散的・質的な判断が行われることがほとんどだと思います)。それならばせめて,事前に「pdが〇〇以上であれば効果があったと判断する」という基準を決めておかないとまずい訳です。残念ながら,そういった基準を明示せずにpdの結果を見て「効果があった」と主張している報告を複数回観測したことがあるのですが…… (しかもpdの値が.80ぐらいのものもありました)。帰無仮説検定の問題点として二値判断である点が批判されることは多いですが,基準を明確にせずに事実上二値判断をするのに比べたら,α = .05と事前に決めた上でp値を使って二値判断をするほうが遥かに健全だと思います。pdの基準については特別な理由やこだわりがないのであれば両側検定のα = .05に相当するpd > .975を採用するのが無難だと思います。※もちろん,二値判断せずにpdを量的に解釈するという方針を採用するのであれば (それができるのであれば),このような基準を決める必要はありません。ここで批判の対象はあくまでもpdを恣意的に解釈して二値判断する行為に対するものであり,pdそのものを批判している訳ではありません。

3. 他の指標(95%確信区間・ベイズファクター)との比較

 他の指標を使ってAさんの仮説を評価したらどうなるでしょうか。試しに確信区間を計算してみたところ,μdiffの95%確信区間は[-0.177, 0.925]であり,0をまたいでいるようです。95%確信区間が0をまたがないことはpd > .975と同値なので,pd = .91のとき0をまたぐのは当たり前ですね。

 続いてベイズファクターではどうでしょうか。ベイズファクターは事前分布に依存しますが,今回は標準化効果量の事前分布をスケールパラメータr = 1/√2のコーシー分布としました (JASPのデフォルトの設定と同じ;Morey & Rouder, 2011およびRouder et al., 2009も参照)。帰無仮説 (効果ゼロ; μdiff = 0) の周辺尤度を分母とするベイズファクターを計算したところBF = 0.605となり,どちらかといえば帰無仮説に軍配が上がるが今回のデータからはどちらともいえない,という結果が得られました (Jeffreysの基準でいえばanecdotalな証拠)。なんにせよ,今回のAさんのデータから「効果がある」と結論づけるのは難しそうですね。
※ここでの分析方法の詳細が気になる方は後述のR & Stanコードを参照してください。

4. まとめ

 p値を使ったら有意にならないけどベイズなら効果があると主張できる,という幻想をぶち壊したくてこの記事を書きました。少なくとも僕は経験上,他の条件がほぼ同じであるにもかかわらずp値では有意にならないがベイズなら効果があると言えるケースに遭遇したことは一度もありません(もしそういう事例があったら教えてほしいです)。そういうことが起こった場合には,分析の設定か解釈を誤っていることを疑ったほうが良いと思います。例えば事前分布を不適切に設定した場合,ベイズファクターが誤ったモデルを選択してしまうといったことは起こり得ます (清水先生のこの記事が参考になります)。逆に,p値では有意なのにベイズファクターで評価するとむしろ帰無仮説に軍配が上がるといったようなことはときどき起こります。ベイズの利点は存在しない効果をあたかも存在するかのように主張できることではなくて,統計モデリングとの相性の良さや優れた仮説評価の方法を提供してくれるところにあると僕は考えています。今回紹介したpdも,適切に使用・解釈する分には豊かな情報を与えてくれる指標の1つだと思われます (議論はあるでしょうが)。

 ちなみに,今回使用したAさんのデータは,帰無仮説が正しい (真の効果がゼロ) という前提のもとで生成したデータです (ちょうどいいp値とpdになるようにseedの調整はしましたが)。つまり,本当は効果がないのにもかかわらず素朴にpdを解釈すると効果があるように見えてしまうことが起こりうるという例示でもありました。p値との対応を考えると,例えばpd > .90で効果ありとみなすという基準を設けたとしたら,単純なモデルにおいてこの基準はα = .20と等価になってしまいます。つまり,帰無仮説が正しかったとしても20%の確率でpd > 90となり「効果あり」と誤って判断してしまうことを意味します。分析の目的によってはこれで良いこともあるかもしれませんが,科学的な研究においてはこの二値判断はかなりリスキーですね。ベイズ的なアプローチには再現性の向上が期待されている側面もありますが,よく考えずになんとなーくベイズを使ってしまうと,p値を使う以上に再現性を低めてしまう危険性があります (査読で弾ければ良いのですが,査読者が必ずしも理解しているとは限らない)。ベイズに限らず,自分が依拠する研究法についてきちんと理解しておくことは非常に大切なことだと思います。僕も勉強し続けます。

5. RとStanのコード

 架空データの生成・分析に使用したRコードとStanコードは以下の通りです。

Rコード: run.R
library(tidyverse)
library(BayesFactor)
library(rstan)

set.seed(21)

# Making data -------------------------------------------------
N <- 60
Y1 <- rnorm(N/2,0,1)
Y2 <- rnorm(N/2,0,1)

# Visualization of data -------------------------------------------------
dg <- data.frame(Y1 = Y1, Y2 = Y2) %>%
  gather(key = "group", value = "value") %>%
  mutate(group = ifelse(group=="Y1", "統制群", "実験群")) %>%
  mutate(group = factor(group, levels = c("統制群", "実験群")))

g <- ggplot(data = dg, aes(x = group, y = value)) +
  theme_bw() +
  geom_jitter(width = 0.2, color = "99BB99", alpha = 0.8, size = 2) +
  stat_summary(fun = "mean", fun.min = function(x)mean(x), fun.max = function(x)mean(x), geom = "crossbar", color = "black", width = 0.3)
plot(g)
ggsave("dataplot.png", g, width = 3, height = 3, dpi = 600)

# Result of t-test -------------------------------------------------
t.test(Y1,Y2) # p = .1625

# MCMC sampling by Stan -------------------------------------------------
sm_t <- stan_model(file = "ttest.stan")

standata <- list(
  N1 = length(Y1),
  N2 = length(Y2),
  Y1 = Y1,
  Y2 = Y2
)

fit <- sampling(sm_t, data = standata, seed = 8823,
                chains = 4, iter = 25500, warmup = 500, thin = 1)

s <- summary(fit)$summary
e <- rstan::extract(fit) %>% as.data.frame()

# Summary of the posterior of mu_diff -------------------------------------------------
s["mu_diff",] %>% round(3)

# Result of probability of direction (pd) -------------------------------------------------
max(mean(e$mu_diff > 0), 1-mean(e$mu_diff > 0)) # pd = .91

# Visualization of posterior of mu_diff -------------------------------------------------
dg <- data.frame(
  mu_diff = density(e$mu_diff)$x,
  density = density(e$mu_diff)$y
)

g <- ggplot(dg, aes(x = mu_diff, y = density)) +
  theme_bw() +
  geom_area(fill = "#888888") +
  geom_area(data = dg[dg$mu_diff>0,], fill = "#DD8888") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  ylab("事後確率密度")

plot(g)

ggsave("posterior.png", g, width = 3, height = 3, dpi = 600)

# Result of Bayes factor -------------------------------------------------
ttestBF(x = Y1, y = Y2, rscale = 1/sqrt(2)) # BF = 0.605, in favor of the null hypothesis

Stanコード: ttest.stan
data{
  int N1;
  int N2;
  real Y1[N1];
  real Y2[N2];
}

parameters{
  vector[2] mu;
  vector<lower = 0>[2] sigma;
}

model{
  Y1 ~ normal(mu[1], sigma[1]);
  Y2 ~ normal(mu[2], sigma[2]);
}

generated quantities{
  real mu_diff = mu[2] - mu[1];
}

posted by mutopsy at 00:00 | 統計