Stanでベイジアンメタ分析（一対比較の固定効果モデル）

はじめに

『Network Meta-Analysis for Decision-Making (Statistics in Practice)』は，ベイズ統計学ネットワークメタ分析を学ぶ上で良い本ですが，紹介されているコードがWinBUGSだったりします（これ以外の書籍もWinBUGSが多い）。Mac & Stanユーザーとしては，WinBUGSはきついので，第2章で紹介されている一対比較(つまりネットワークメタ分析ではない)データの固定効果モデルをStanコードに書き直してみました。

Network Meta-Analysis for Decision-Making (Statistics in Practice)

使用するパッケージ

以下のパッケージを使います。これら一式がインストールされたDockerfileも公開しているので，こちらの記事も参照ください。

library(rstan)
library(tidybayes)
library(tidyverse)
library(bayesplot)
library(loo)

データ

使用するのは，『Network Meta-Analysis for Decision-Making』の２章で紹介されている血栓溶解薬のデータです(Caldwell et la., 2005のデータ)。データセット全体の中からACC t-PA(治療0,今回はこれをコントロールにする)とPTCA(治療1)を比較した11の試験を使います。つまり，一対比較(Pairwise comparison)データです。

studyは研究のID，treatmentは治療の種類(0=ACC t-PA，1=PTCA)，deadは死者数，sampleSizeはその治療に参加した患者数です。studyNameは第1著者の姓か研究プロジェクト名， studyYearは論文の出版年です。

Pairwise data for meta-analysis

以下のような感じのデータです。

# A tibble: 6 x 7
  study treatment  dead sampleSize treatmentName studyName studyYear
  <dbl>     <dbl> <dbl>      <dbl> <chr>         <chr>         <dbl>
1     1         0     3         55 Acc_tPA       Ribichini      1996
2     2         0    10         94 Acc_tPA       Garcia         1997
3     3         0    40        573 Acc_tPA       GUSTO-2        1997
4     4         0     5         75 Acc_tPA       Vermeer        1999
5     5         0     5         69 Acc_tPA       Schomig        2000
6     6         0     2         61 Acc_tPA       LeMay          2001

固定効果モデルのStanコード

まず，data{}ブロックにおいて，使用するデータの定義をしています。教科書は行列形式でdeadやsampleSizeを読み込む形式ですが，少し今後の拡張を考えると面倒です(WinBUGSでやりやすいこととStanでやりやすいことは微妙に違ったりします)。まずlong型のデータセットにしてから（上記のデータはすでにそうなっています），各列をStanに読み込ませます。

parameters{}ブロックでは，推定するパラメータとして，mu(各研究におけるベースライン=各研究での治療0の効果)とd01(治療0と治療1の差=治療0よりも治療1が優れるor劣る効果の大きさ)を準備しています。

model{}ブロックでは，死者数が二項分布に従うとして，死者数が，binomial_logit(試験の参加人数，死亡確率を構成する式)から生成されます。なお，死亡確率を構成する式は，治療が0の場合，はその試験のmuのみで，治療1の場合は，その試験のmuにd01を足したものになります。d01とmuの事前分布としては，幅のひろーい正規分布としました。

generated quantities{}ブロックでは，d01から治療0に対する治療1のオッズ比や有害の確率を計算したり，モデル比較用の対数尤度(log_lik)も計算しています。

Stan code of pairwise data using fixed effect mode ...

パラメータ推定

Stanコードが書けましたので，早速，コンパイル＆サンプリングをします。

ld = length(study)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
options(max.print = 99999)
fit_fixed <- stan("netmeta_pairwise_fixed_effect.stan",data=list(ld = ld, nt = 2, ns = 11, study = study, treatment = treatment, dead = dead, sampleSize = sampleSize), chains = 4, iter = 5500, warmup = 500, thin = 1)

推定結果

結果を簡単に確認します。

print(fit_fixed,digit = 3)

見にくいので，一部の結果のみを示します。若干ズレはありますが，教科書とほぼ同じ推定値になりました（関心のあるパラメータのみ掲載）。Rhatやn_effからもサンプリングも問題なさそうです。

	mean	se_mea	sd	2.5%	97.5%	n_eff	Rhat
d01	-0.2344	0.0009	0.1194	-0.4689	-0.0026	17330	0.9999
OR01	0.7967	0.0007	0.0954	0.6257	0.9974	17086	0.9999
Prob_harm	0.0236	0.0012	0.1520	0.0000	0.0000	17264	0.9999

なお，収束判定は以下のように可視化してもできます（コードのみで図は割愛します）。R hat，トレースプロット，自己相関，有効サンプルサイズの順番です。

stan_rhat(fit_fixed, pars = c("d01", "mu"))
stan_trace(fit_fixed, pars = c("d01", "mu"),inc_warmup=T)
stan_ac(fit_fixed, pars = c("d01", "mu"))
stan_ess(fit_fixed, pars = c("d01", "mu"))

治療0に対する治療1のオッズ比をプロットしてみます。ACC t-PAと比べて，PTCAが死亡率を下げることが分かりますね。

fit_fixed %>% 
  spread_draws(OR01) %>% 
  ggplot(aes(x = OR01)) +
  stat_density(fill = "gray75") +
  stat_pointintervalh(aes(y = -0.05), point_interval = mean_qi, .width = .95) +
  annotate("text", label = "mean, 95% quantile intervals", x = 1.1, y = -0.05, hjust = 0, vjust = 0.3) +
  xlim(0.5,1.45)+
  labs(x="Odds ratio")

単体では意味がないですが，モデル比較をすることもあるかと思い，Stanコードでは対数尤度も計算をしています。以下のコードでWAICも算出できます。

log_like <- extract_log_lik(fit_fixed)
waic(log_like)

これで，メタ分析の基本形である固定効果モデルをStanで実行できました！