ベイズ解析を使った臨床研究の原稿で、結果の節に事後平均と95%信用区間だけが並び、事前分布をどう決めたのかも、計算が収束したのかも書かれていない、ということがあります。その区間がどの事前分布とどの計算手続きから出たのかが欠けたままでは数値の意味を評価できないため、ここでは2群ランダム化試験を模したランダムに生成したデータを使い、事前分布の指定、収束診断、事後量をどの数値でそろえるかを扱います。
報告の起点になる推定対象と記号
解析の対象は、対照群250例、治療群250例からなる500例のランダムに生成したデータセットで、アウトカムは28日死亡としました。対照群では82例、治療群では55例に死亡が観察され、リスクはそれぞれ0.328と0.220、リスク差は-0.108であり、ここから推定したいのは治療群と対照群の対数オッズ比です。
この対数オッズ比を$\theta$、対照群の対数オッズを$\alpha$と置くと、観測データ$\mathbf{y}$のもとでの事後分布は次のように書けます。
$$p(\theta, \alpha \mid \mathbf{y}) \propto L(\theta, \alpha; \mathbf{y})\, p(\theta)\, p(\alpha) \tag{1}$$
式(1)は、データから計算される尤度$L(\theta, \alpha; \mathbf{y})$に、$\theta$と$\alpha$それぞれの事前分布を掛けたものが事後分布に比例することを表しています。$p(\theta)$が治療効果についての事前分布、$p(\alpha)$が対照群のリスク水準についての事前分布です。
この式の形から、報告に必要な要素は左辺と右辺に分かれます。事前分布の記載が抜けた原稿は、式(1)の右辺の半分を伏せたまま左辺だけを見せていることになり、同じ区間が別のデータから出たものかどうかを読み手が判断できません。
なお、事前分布の幅が結果をどう動かすかという論点は、生存時間アウトカムでも同じ形で現れます。区間の狭さがデータの情報量から来ているのか事前分布の仮定から来ているのかという区別については、ベイズ生存時間解析：狭くなった区間は情報か仮定かで扱っています。
事前分布の幅と由来の書き方
事前分布を「無情報事前分布を用いた」とだけ書いた原稿は珍しくありませんが、その一語からはどの範囲の効果をあらかじめ許容したのかが読み取れません。今回の解析で$\theta$に置いた弱情報事前分布Normal(0, 1)は、オッズ比の尺度に直すと中央値1、95%の範囲が0.141から7.099に広がる設定であり、この幅は治療によって死亡オッズが7倍になる可能性まで事前に認めていることを意味します。
そのため、事前分布の記載は分布名とパラメータだけで終えず、オッズ比やハザード比の尺度へ引き直した中央値と95%区間まで書くことを推奨します。同じNormalでも標準偏差を0.354に絞った懐疑的事前分布は、オッズ比の95%範囲が0.500から2.001となり、2倍を超える効果をあらかじめ強く疑う設定になるため、読み手が判断できるのは分布名ではなくこの区間の広さです。
この標準偏差0.354という値自体も、オッズ比が2倍を超える効果は考えにくいという条件から逆算して決めたものです。log 2は0.693で、これを正規分布の97.5%点1.96で割ると0.354になり、事前分布の95%上限がちょうどオッズ比2に来るので、読み手は設定の根拠を一行で追えます。
由来についても同じことが言えます。今回は比較のために治療有利側へ中心を置いた楽観的事前分布Normal(log 0.7, 0.354)も用意しましたが、この設定ではオッズ比が1を下回る事前確率が0.843となり、データを見る前から治療有利に84%の確信を置いた状態から解析が始まります。
この由来の書き方について、Sungら（2005）は、臨床研究のベイズ解析を報告する際に含めるべき項目を検討し、事前分布の指定とその根拠を中核に置いています。またZampieriら（2021）は、集中治療領域の試験を再解析した総説で、事前分布を解析前に定めないと結果に合わせて選んだという疑いを招くと指摘し、懐疑的な設定と楽観的な設定を並べる実務を示しているため、データを見た後に事前分布を選び直してその経緯を書かないまま一つの結果だけを載せる進め方は推奨できません。
事前分布の違いが事後量へ与える影響
三つの事前分布で同じデータをマルコフ連鎖モンテカルロ法にかけると、事後中央値のオッズ比は弱情報事前分布で0.584、懐疑的事前分布で0.659、楽観的事前分布で0.606となりました。95%信用区間はそれぞれ0.388から0.869、0.469から0.925、0.422から0.855で、点推定の差に比べて区間の重なりは大きく出ています。
密度を0（オッズ比1）の周りに集めた事前分布は推定値を効果なし側へ引き戻すため、標準偏差を0.354に絞った懐疑的事前分布では、データが示す0.578という観測オッズ比が0.659まで戻り、三つの中で最も1に近くなりました。楽観的事前分布は中心をオッズ比0.7に置いているので、0.578を0.7の方向へ引き上げ、事後中央値は0.606と弱情報事前分布の0.584より1に近づきます。
これに対して、事前分布の幅が尤度の情報量に比べて十分に広いとき、事後分布はほぼ尤度の形をなぞります。弱情報事前分布のもとでの事後中央値0.584が最尤推定のオッズ比0.578とほとんど変わらないのはこのためで、95%信用区間0.388から0.869も最尤推定の95%信頼区間0.388から0.861とほぼ重なりますが、一致しているのは区間の位置であって、区間が意味する内容ではありません。
収束診断で確認する値と判定の目安
事後量はマルコフ連鎖モンテカルロ法の出力から計算するので、連鎖が目的の事後分布を十分に探索できていなければ、どれだけ小数点以下を並べても意味を持ちません。今回は過分散な初期値から4本の連鎖を8,000回ずつ走らせて前半4,000回を捨て、16,000個の標本を残しましたが、GelmanとRubin（1992）の収束診断は、連鎖間の分散と連鎖内の分散を比べる考え方です。
その比較を一つの値にまとめたものが$\hat{R}$です。
$$\hat{R} = \sqrt{\frac{\widehat{\mathrm{var}}^{+}}{W}}, \qquad \widehat{\mathrm{var}}^{+} = \frac{n-1}{n}W + \frac{1}{n}B \tag{2}$$
式(2)の$W$は連鎖内分散の平均、$B$は連鎖間分散、$n$は各連鎖の長さで、これは分割や順位正規化を加える前の基本形です。連鎖がまだ別々の領域にとどまっていると$B$が大きくなり、$\hat{R}$は1から離れます。
この診断について、Vehtariら（2021）は順位正規化と分割による改良を示し、判定の目安として1.01を挙げています。今回報告する$\hat{R}$もこの改良版で計算しており、三つの設定でそれぞれ1.0006、1.0032、1.0023とどれも目安を下回りましたが、目安を下回ったことは連鎖が事後分布全体を探索し終えた保証ではなく、探索していない領域があれば、その領域のことは連鎖間の比較には現れません。
そこで、$\hat{R}$と合わせて有効サンプルサイズとモンテカルロ標準誤差も報告します。今回の有効サンプルサイズは対数オッズ比について1,576、1,723、1,555で、16,000個の標本から実質1割前後の独立な情報が得られている状態であり、対応するモンテカルロ標準誤差は0.005、0.004、0.004、事後標準偏差は0.20前後なので、報告する小数第2位までは計算の揺らぎで動かず、小数第3位はこの誤差の範囲で動きます。
逆に、提案分布の刻み幅を100分の1に縮めて同じモデルを走らせると、採択率は0.956まで上がるのに$\hat{R}$は3.00、有効サンプルサイズは260まで落ちます。刻み幅が小さすぎる連鎖は現在位置のすぐ近くばかりを提案するため、ほとんどの提案が受理されながら分布全体を回れていません。
事後量の選び方と閾値の決め方
診断を通過した標本から何を取り出すかは、研究目的に対応する量を先に決めておくことで定まります。事後中央値と95%信用区間に加えて、ベイズ解析で報告できるのは特定の閾値$c$を下回る事後確率です。
$$\Pr(\mathrm{OR} < c \mid \mathbf{y}) = \int_{-\infty}^{\log c} p(\theta \mid \mathbf{y})\, d\theta \tag{3}$$
式(3)は、事後分布のうち対数オッズ比が$\log c$より小さい部分の面積が、オッズ比が$c$を下回る事後確率になることを表しています。$c=1$と置けば、治療がわずかでも有利である確率になります。
実際に、今回のデータでオッズ比が1を下回る事後確率を求めると、弱情報事前分布で0.997、懐疑的事前分布で0.991、楽観的事前分布で0.998となりました。ただし、この確率は効果の方向についての確率にとどまり、臨床的に意味のある大きさかどうかは別に評価する必要があるので、閾値を0.9へ動かすと同じ順で0.984、0.963、0.988となり、懐疑的事前分布のもとでは2ポイント以上下がります。
この差をどう扱うかは、閾値をいつ決めたかで変わります。解析計画の段階で0.9と定めていれば三つの事前分布で0.96を超えるという結果は解釈しやすくなりますが、結果を見てから0.9を選んだ場合その数値は報告の見せ方の選択になるため、閾値を先に決めて根拠まで書くことを推奨します。
また、効果がほぼないと見なす範囲を置いて、その範囲に事後分布が入る確率を報告する方法もあります。オッズ比0.9から1.111を実質的に差がない範囲とすると、そこに入る事後確率は弱情報事前分布で0.015、懐疑的事前分布で0.037、楽観的事前分布で0.012であり、差がある方向の確率だけでなく差がない範囲の確率まで並べると、読み手が自分の基準で判断できます。
確率だけでなく点推定として何を載せるかにも、同じ種類の注意が要ります。弱情報事前分布のもとでの対数オッズ比の事後平均は-0.541で、これを指数変換すると0.582になり、オッズ比の尺度で直接求めた事後中央値0.584とはわずかに離れるため、対数尺度で推定して指数変換する量では、尺度を変えても位置が保たれる中央値を報告した方が読み取りやすくなります。
報告ガイドラインと公表研究での記載
ここまで挙げてきた項目は、報告ガイドラインの中でも整理されています。Kruschke（2021）のBayesian Analysis Reporting Guidelinesは、事前分布の明記、収束診断とモンテカルロ標準誤差の報告、事後分布の要約という並びで記載事項を定め、診断値を本文か補遺のどちらかに必ず置くことを求めており、Sungら（2005）も臨床研究に絞った形で同様の項目を導出しています。
この記載事項に沿った公表研究の例としては、Goligherら（2018）による重症急性呼吸窮迫症候群の体外式膜型人工肺に関する事後ベイズ解析があります。この解析では複数の事前分布のもとで死亡が減る事後確率を並べて示しており、事前分布の選択に対して結果がどれだけ動くかが読み取れる形になっています。
もっとも、どこまでの設定を並べるかは判断が要ります。事前分布を無数に試して都合のよい結果だけを載せるのは避けるべきですが、一つの設定だけでは結論をどれだけ支えているかが分からないので、解析計画の段階で弱情報、懐疑的、楽観的という程度の設定を決めておき、全ての結果を表として出す形が扱いやすくなります。
Rで再現する事前分布と事後量
ここからは、事前分布と事後量をbase Rとstatsだけで再現します。収束診断の数値はランダムウォークMetropolis法の実装から出したものでこの節のコードには含めず、代わりに同じ事後分布を格子求積で確かめる形にしました。
set.seed(42)
n_per_arm <- 250
treatment <- rep(c(0, 1), each = n_per_arm)
risk_true <- ifelse(treatment == 1, 0.24, 0.30)
death <- rbinom(n = length(treatment), size = 1, prob = risk_true)
fit_ml <- glm(death ~ treatment, family = binomial(link = "logit"))
ml_coef <- summary(fit_ml)$coefficients["treatment", ]
cat(sprintf("OR=%.3f 95%%CI=%.3f-%.3f p=%.4f\n", exp(ml_coef[1]),
exp(ml_coef[1] - qnorm(0.975) * ml_coef[2]),
exp(ml_coef[1] + qnorm(0.975) * ml_coef[2]), ml_coef[4]))#> OR=0.578 95%CI=0.388-0.861 p=0.0071次に、式(1)の右辺にあたる対数事後密度を関数として書き、$\alpha$と$\theta$の値を格子状に並べました。二値アウトカムの2群比較では未知数が2個しかないので、格子の上で密度を計算して足し上げる求積でも事後分布を十分な精度で求められます。
events_control <- sum(death[treatment == 0])
events_treated <- sum(death[treatment == 1])
ml_log_or <- unname(coef(fit_ml)["treatment"])
ml_se <- unname(sqrt(vcov(fit_ml)[2, 2]))
alpha_grid <- seq(coef(fit_ml)[1] - 8 * sqrt(vcov(fit_ml)[1, 1]),
coef(fit_ml)[1] + 8 * sqrt(vcov(fit_ml)[1, 1]), length.out = 801)
theta_grid <- seq(ml_log_or - 8 * ml_se, ml_log_or + 8 * ml_se, length.out = 801)
log_posterior <- function(alpha, theta, prior_mean, prior_sd) {
eta_control <- alpha
eta_treated <- alpha + theta
events_control * eta_control - n_per_arm * log1p(exp(eta_control)) +
events_treated * eta_treated - n_per_arm * log1p(exp(eta_treated)) +
dnorm(alpha, mean = 0, sd = 10, log = TRUE) +
dnorm(theta, mean = prior_mean, sd = prior_sd, log = TRUE)
}この関数を使って、$\alpha$方向に足し上げた$\theta$の周辺事後分布から、中央値、95%信用区間、オッズ比が1を下回る事後確率を取り出します。事前分布の平均と標準偏差だけを変えて三つの設定を並べました。
posterior_summary <- function(prior_mean, prior_sd) {
log_dens <- outer(alpha_grid, theta_grid,
function(a, t) log_posterior(a, t, prior_mean, prior_sd))
weight <- colSums(exp(log_dens - max(log_dens)))
weight <- weight / sum(weight)
cdf <- cumsum(weight)
quantile_or <- function(p) exp(approx(x = cdf, y = theta_grid, xout = p)$y)
c(quantile_or(0.5), quantile_or(0.025), quantile_or(0.975), sum(weight[theta_grid < 0]))
}
for (prior in list(c(0, 1), c(0, 0.354), c(log(0.7), 0.354))) {
summary_or <- posterior_summary(prior_mean = prior[1], prior_sd = prior[2])
cat(sprintf("median OR=%.3f 95%%CrI=%.3f-%.3f P(OR<1)=%.3f\n", summary_or[1],
summary_or[2], summary_or[3], summary_or[4]))
}#> median OR=0.588 95%CrI=0.396-0.868 P(OR<1)=0.996
#> median OR=0.660 95%CrI=0.467-0.930 P(OR<1)=0.991
#> median OR=0.604 95%CrI=0.427-0.852 P(OR<1)=0.998この出力は上から弱情報、懐疑的、楽観的の順に並んでおり、懐疑的事前分布の行だけ中央値が1寄りで区間の上端も0.930まで伸びています。この求積で得た値は、先に示したマルコフ連鎖モンテカルロ法による事後中央値0.584、0.659、0.606と小数第2位までおおむね一致しており、標本の数が足りていることの確認になります。
この二つの計算経路の値は小数第3位で食い違い、そこが読み違えやすい点ですが、差の由来は事後分布そのものの不確実性ではなく計算の仕方なので、モンテカルロ標準誤差0.005の範囲に収まっているかどうかで判断し、報告する桁数もその水準に合わせます。なお、図1と図3はここで作った周辺事後分布から描き、図2はマルコフ連鎖モンテカルロ法の標本をそのまま描いたものです。
研究計画に持ち帰る判断
ベイズ解析の報告で読み手が確かめたいのは、事前分布がどの範囲の効果を許容していたか、計算が事後分布を探索し終えているか、そして報告された確率がどの閾値に対するものかという点に絞られます。今回のデータでは、事前分布を懐疑的な設定に変えると事後中央値のオッズ比が0.584から0.659へ動き、オッズ比が0.9を下回る事後確率は0.984から0.963へ下がりました。
この動き幅を解析計画の段階で見込んでおくかどうかが、実際の研究での分かれ目になります。事前分布と閾値を結果が出る前に決めて表の形を固め、$\hat{R}$と有効サンプルサイズを補遺に載せる準備をしておけば、事後確率の数値をめぐる説明は後から作らずに済みます。
実際のデータでは、先行研究から引ける効果の大きさ、アウトカムのイベント数、事前に定めた解析計画の記載まで含めて判断することになり、Dr.データサイエンスでは、こうした解析方針のご相談を承っています。
参考文献
- Sung L, Hayden J, Greenberg ML, Koren G, Feldman BM, Tomlinson GA. Seven items were identified for inclusion when reporting a Bayesian analysis of a clinical study. Journal of Clinical Epidemiology. 2005;58(3):261-268.
- Kruschke JK. Bayesian Analysis Reporting Guidelines. Nature Human Behaviour. 2021;5:1282-1291.
- Vehtari A, Gelman A, Simpson D, Carpenter B, Bürkner PC. Rank-Normalization, Folding, and Localization: An Improved R-hat for Assessing Convergence of MCMC (with Discussion). Bayesian Analysis. 2021.
- Gelman A, Rubin DB. Inference from Iterative Simulation Using Multiple Sequences. Statistical Science. 1992;7(4):457-472.
- Zampieri FG, Casey JD, Shankar-Hari M, Harrell FE Jr, Harhay MO. Using Bayesian Methods to Augment the Interpretation of Critical Care Trials. An Overview of Theory and Example Reanalysis of the Alveolar Recruitment for Acute Respiratory Distress Syndrome Trial. American Journal of Respiratory and Critical Care Medicine. 2021;203(5):543-552.
- Goligher EC, Tomlinson G, Hajage D, Wijeysundera DN, Fan E, Jüni P, Brodie D, Slutsky AS, Combes A. Extracorporeal Membrane Oxygenation for Severe Acute Respiratory Distress Syndrome and Posterior Probability of Mortality Benefit in a Post Hoc Bayesian Analysis of a Randomized Clinical Trial. JAMA. 2018;320(21):2251-2259.
