症例数の限られた研究でベイズ法を使うと、95%信用区間が頻度論の信頼区間より狭く出ることがあります。区間が狭ければ結論は明確に見えやすく、報告する側にとっても扱いやすい結果になります。
ただし、そこで考えなければならないのは、その狭さがどこから生まれたのかという点です。データそのものが増えたわけではなく、事前分布という形で外部から加えた情報が効いています。その影響がどの程度なのかは計算できるため、今回は公開データを使って確かめました。
事前分布は情報を1つ足している
対数ハザード比を$\beta$とします。データだけから推定した時の標準誤差を$\hat{\sigma}$、事前分布として$\beta \sim N(0, \tau_{\mathrm{prior}}^{2})$を置いたとします。事後分布が正規分布で近似できる範囲では、精度（分散の逆数）が足し算になります。
$$\frac{1}{\tau_{\mathrm{post}}^{2}} \approx \frac{1}{\tau_{\mathrm{prior}}^{2}} + \frac{1}{\hat{\sigma}^{2}} \tag{1}$$
式(1)の$\tau_{\mathrm{post}}$は事後分布の標準偏差です。右辺を見ると、第1項は事前分布から持ち込まれる精度、第2項はデータが持つ精度に対応しています。事後分布の精度はこの2つの和になるため、事前分布を狭く設定するほど、事後分布もそれに応じて狭くなります。
ここで重要なのは、右辺の2項が同じ単位で足されていることです。データの精度はイベント数が増えるにつれて高くなりますが、事前分布の精度は研究者が決めた$\tau_{\mathrm{prior}}$だけで決まり、データを1例も見なくても最初から存在します。したがって、区間が狭いという事実だけを見ても、その狭さがデータによるものなのか、事前分布によるものなのかは判断できません。
一方で、両者の比を計算すれば、その内訳を分けて考えることができます。
平坦な事前分布ならCoxと一致する
ここでは、ドイツ乳がん研究グループの試験データを使います。survivalパッケージのgbsgデータセットにはリンパ節転移陽性の乳がん686例が収録されており、無再発生存のイベントは299件です。Schumacherらが1994年にJournal of Clinical Oncologyで報告した2×2試験で、ホルモン療法の有無が割り付けられています。
ハザードが時間によらず一定という指数モデルを置き、ホルモン療法群のハザードを$\lambda_0 e^{\beta}$、対照群を$\lambda_0$とします。$\lambda_0$を積分で消すと、$\beta$の事後分布は閉じた形で書けます。
$$p(\beta \mid \mathcal{D}) \propto \frac{e^{\beta D_1}}{\left(T_0 + e^{\beta} T_1\right)^{D_0 + D_1}}\, p(\beta) \tag{2}$$
式(2)の$D_0$と$D_1$は各群のイベント数、$T_0$と$T_1$は各群の観察年数の合計を表します。$p(\beta)$が事前分布で、それ以外の部分が尤度に相当します。このモデルでは積分や乱数生成を行う必要がなく、$\beta$を細かい格子に区切って式を評価するだけで事後分布を求められます。
平坦な事前分布を置いた場合、事後中央値は0.699、95%信用区間は0.545から0.890でした。一方、同じデータにCox比例ハザードモデルを当てはめると推定値は0.695、信頼区間は0.544から0.888です。Cox比例ハザードモデルはハザードの形を指定しませんが、この試験ではハザードがおおむね一定であり、指数モデルとの差は3桁目に現れる程度でした。
この結果から分かるのは、ベイズ法そのものが特別な推定値を生み出しているわけではないということです。平坦な事前分布を使えば、事後分布はほぼ尤度そのものになります。それでは、ベイズ解析で時に見られる狭い信用区間は、どこから生じるのでしょうか。
46イベントでは幅の3割が事前分布から来る
信用区間が狭ければ推定の精度が高い、と考えたくなります。しかし、ベイズ解析では必ずしもそうとは限りません。事前分布を狭く設定するだけでも、データを1例追加することなく区間を縮められるからです。
そこで、686例のうち120例を無作為に取り出し、イベント46件のデータを作りました。平坦な事前分布では信用区間が0.375から1.384となり、幅は1.009です。これに対して、対数ハザード比の標準偏差0.25の正規分布を事前分布として置くと、区間は0.609から1.315、幅は0.706まで縮まりました。データは同じままですが、区間は約3割狭くなっています。
この内訳は式(1)から計算できます。今回のデータでは標準誤差が0.333なので、データの精度は9.05です。一方、標準偏差0.25の事前分布が持つ精度は16であり、事後分布の精度のうち64%を事前分布が占める計算になります。さらに標準偏差0.10まで狭くすると、その割合は92%に達し、区間の幅は0.366まで縮みました。平坦な事前分布の場合と比べると、およそ3分の1です。
これに対してイベント299件の全データでは、同じ標準偏差0.25の事前分布を置いても、区間の幅は0.344から0.327へと5%縮むだけで、事前分布の寄与も20%にとどまります。つまり、同じ強さの事前分布を使っていても、データが少ない研究ほどその影響は大きくなります。
狭くなった区間は中心もズレている
もう1つ注意したいのは、区間の幅だけに注目すると、事前分布によって推定値の中心そのものが動いていることを見落としやすい点です。
イベント46件の場合、平坦な事前分布で0.743だった中央値は、標準偏差0.25では0.898、標準偏差0.10では0.975まで動きました。区間が約3分の1に縮むのと同時に、推定されたハザード比そのものもほぼ1に近づいています。イベント299件の全データでさえ、標準偏差0.10の事前分布を置けば、中央値は0.699から0.867まで変化します。
ハザード比1を中心とする事前分布を置くことは、「効果はない」という情報を解析前に与えることでもあります。その結果として狭い区間が得られたとしても、必ずしも効果を精密に推定できたことを意味しません。単に、解析結果があらかじめ置いた仮定へ強く引き寄せられている場合もあります。Spiegelhalterらが1994年に懐疑的事前分布として整理した方法は、むしろこの性質を意図的に利用するもので、「効果はない」という立場を事前に置いたうえでも効果が残るかを確かめるために使われます。
もちろん、逆方向にも同じことが起こります。先行研究で得られた効果を事前分布に組み込めば、区間は狭くなり、推定値の中心も先行研究の結果へ近づきます。もし先行研究に偏りが含まれていれば、その偏りまで新しい解析に持ち込むことになります。
事前分布の強さをイベント数に換算する
「事前分布の標準偏差を0.25にした」とだけ書かれていても、それがどの程度強い仮定なのかを直感的に理解するのは難しいでしょう。そこで、式(1)の第2項がイベント数に比例することを使うと、事前分布を何イベント分の情報かに言い換えられます。
このデータでは、イベント46件、131件、299件の場合のデータの精度が、それぞれ9.05、27.9、64.2でした。1イベントあたりに直すと0.208です。一方、事前分布の精度は、標準偏差0.25なら16、0.10なら100になります。この値を1イベントあたりの精度で割ると、標準偏差0.25の事前分布は約77イベント分、標準偏差0.10なら約480イベント分の情報に相当します。同じ計算をすると、標準偏差0.50では約19イベント分、1.00では約5イベント分です。
私なら、事前分布の設定を報告する際には、この換算値も併記します。「標準偏差0.10」と書くだけでは、その強さを読者が判断するのは簡単ではありません。しかし、「外部から約480イベント分に相当する情報を加えた」と表現すれば、研究そのもののイベント数と並べて重みを比較できます。46イベントのデータに480イベント分の仮定を加えた解析なら、得られた結果が事前分布の要約にかなり近くなることも理解しやすくなります。
ただし、この換算は正規近似に基づくため、あくまで目安として扱う必要があります。事後分布が大きく歪んでいる場合や、事前分布がハザード比1から離れた位置に置かれている場合には、精度の単純な足し算だけでは十分に説明できません。それでも、標準偏差の数字だけを示すより、事前分布が持つ情報量を読者に伝えやすい方法です。
指数モデルなら事後分布は式で書ける
ここまで示した数値を再現する手順も載せておきます。survivalパッケージのgbsgデータセットを使い、式(2)を格子上で直接評価しています。マルコフ連鎖モンテカルロ法は使わないため、収束診断も必要ありません。
library(survival)
yrs <- gbsg[["rfstime"]] / 365.25
ev <- gbsg[["status"]]; grp <- gbsg[["hormon"]]
beta <- seq(-2, 1.2, length.out = 4001)
cred <- function(rows, sigma) { # rows は使う症例の添字
d0 <- sum(ev[rows][grp[rows] == 0]); t0 <- sum(yrs[rows][grp[rows] == 0])
d1 <- sum(ev[rows][grp[rows] == 1]); t1 <- sum(yrs[rows][grp[rows] == 1])
lp <- beta * d1 - (d0 + d1) * log(t0 + exp(beta) * t1) +
dnorm(beta, 0, sigma, log = TRUE) # 対数事後分布
cdf <- cumsum(exp(lp - max(lp))); cdf <- cdf / max(cdf)
round(exp(approx(cdf, beta, c(0.025, 0.5, 0.975))[["y"]]), 3)
}この関数を、全データと120例を無作為抽出した部分集合の両方に当てはめます。
set.seed(42)
sub <- sample(nrow(gbsg), 120); all <- seq_len(nrow(gbsg))
out <- rbind("299 events, flat" = cred(all, 10), "299 events, SD 0.25" = cred(all, 0.25),
" 46 events, flat" = cred(sub, 10), " 46 events, SD 0.25" = cred(sub, 0.25))
colnames(out) <- c("lower", "median", "upper"); out#> lower median upper
#> 299 events, flat 0.545 0.699 0.890
#> 299 events, SD 0.25 0.603 0.750 0.930
#> 46 events, flat 0.375 0.743 1.384
#> 46 events, SD 0.25 0.609 0.898 1.315この出力は、行を比較して読むと意味が分かりやすくなります。同じ標準偏差0.25の事前分布でも、299イベントでは中央値が0.699から0.750へ動く程度なのに対し、46イベントでは0.743から0.898まで変化しています。事前分布そのものの強さは同じなので、この差を生んでいるのはデータ側の情報量です。
一方、平坦な事前分布の行は少し注意して読む必要があります。ここでは標準偏差10の正規分布を「平坦」と呼んでいますが、厳密な意味での一様分布ではありません。ただし、この範囲では事後分布への影響が3桁目より小さいことを確認したうえで、比較の基準として使っています。また、パラメトリック生存モデルで扱ったように、指数分布を置くこと自体も1つの仮定です。この試験ではCoxモデルの推定値とほぼ一致しましたが、ハザードが時間とともに大きく変化するデータであれば、事前分布以前にこのモデルの前提から確認し直す必要があります。
事前分布を書かない論文は読めない
Kruschkeが2021年にNature Human Behaviourで公表したBayesian Analysis Reporting Guidelinesでは、事前分布そのものだけでなく、それを選んだ理由と、事前分布を変えた時に結果がどう動くかを報告するよう求めています。事後分布の要約しか示されていなければ、信用区間の狭さがデータから得られた情報によるものなのか、それとも事前に置いた仮定によるものなのかを読者は判断できません。
では、感度分析では事前分布をどこまで動かせばよいのでしょうか。「事前分布を変えても結果は変わらなかった」と報告していても、試した範囲がごく狭ければ、十分な感度分析とは言えません。図2のように、実質的に無情報とみなせる設定から、その分野で正当化できる最も強い設定までを並べて示すと、結果がどの程度事前分布に依存しているのかを把握しやすくなります。これは感度分析で扱った考え方と同じで、どの範囲まで動かすのかという基準も含め、解析前に決めておくことが重要です。
こうした報告は手間が増えるように感じるかもしれませんが、事前分布を明示できること自体はベイズ解析の利点でもあります。どの仮定を入れ、それが結果にどの程度影響したのかを数値として開示できる解析の形式は、頻度論ではそれほど一般的ではありません。
私が事前分布を決める順序
私が事前分布を考える時、最初に確認するのはイベント数です。イベントが200件を超える研究なら、標準偏差0.5程度の弱情報事前分布を置いても事後分布はほとんど変わらないため、事前分布の細かな設定に多くの時間を使う必要はありません。反対に、イベント数が50件を下回るような研究では、事前分布の選び方が結論を大きく左右します。そのため、この段階では解析計画を作る時点で設定を決め、その根拠も記載しておくようにしています。
次に確認するのは、事前分布を何に基づいて設定するかです。先行研究のメタ解析があるなら、その推定値を利用する方法がありますが、そのためには先行研究と自分の研究で患者背景が十分に近いことが前提になります。外部情報として使える根拠が乏しいのであれば、まずハザード比1を中心とする弱情報事前分布を置き、感度分析で幅を変えて結果の安定性を確認します。
私がベイズ法を勧めるのは、イベント数が少なく、頻度論による区間推定では幅が広すぎて十分な判断が難しい一方で、外部から取り込める情報の出どころが明確な場合です。逆に、外部情報の根拠を説明できないまま、信用区間を狭くすることだけを目的として事前分布を使う方法は適切ではありません。区間の狭さは解析手法そのものの性能を示しているのではなく、どれだけ強い仮定を解析に入れたかを反映しているからです。予測を目的とした解析で手法を比較する場合については、機械学習的生存時間解析で扱っています。
実際の研究でどのような事前分布を置くべきか、また感度分析をどこまで示すべきかは、イベント数だけで決められるものではありません。先行研究の質や、自分の研究との類似性、さらに投稿先が求める報告内容まで含めて判断する必要があります。Dr.データサイエンスでは、こうしたベイズ解析の設計や報告方法についてのご相談を承っています。
参考文献
Spiegelhalter DJ, Freedman LS, Parmar MKB. Bayesian approaches to randomized trials. Journal of the Royal Statistical Society Series A. 1994;157(3):357-387.
Schumacher M, Bastert G, Bojar H, et al. Randomized 2×2 trial evaluating hormonal treatment and the duration of chemotherapy in node-positive breast cancer patients. Journal of Clinical Oncology. 1994;12(10):2086-2093.
Ibrahim JG, Chen MH, Sinha D. Bayesian Survival Analysis. Springer; 2001.
Kruschke JK. Bayesian analysis reporting guidelines. Nature Human Behaviour. 2021;5(10):1282-1291.
Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. Bayesian Data Analysis. 3rd ed. CRC Press; 2013.
