Matching Adjusted Indirect Comparison (MAIC)という比較の手法があることを知った. 直接比較されていない薬剤の効果を間接的に比較したい場合, よく知られているのはネットワークメタアナリシス(Network Meta-analysis : NMA)だと思われる. NMAは, 関心のある薬剤と関連する比較研究を集めて, 公表されている集団レベルのデータ(Aggregate Data : AgD)を繋ぎ合わせて間接比較する. ただし, 普通のNMAでは, 研究間に効果修飾因子の異質性がある場合, バイアスを生んでしまう. もしも, 関心のある薬剤の研究で, 個人レベルのデータ(Individual Patient Data : IPD)があるならば, その問題を解決手法がある. その間接比較は, Population-adjusted indirect comparisonsと呼ばれ, その中には, Matching Adjusted Indirect Comparison (MAIC)というものと, Simulated Treatment Comparison (STC)というものがある. 今回は前者のMAICについて学ぶ必要があったので, MAICが具体的にどういう処理をしているのか、それは何を意図しているのかを理解するため, NICE(The National Institute for Health and Care Excellence)のTSD(Technical Support Documents)の18のAppendixのシミュレーションと, MAICの元論文のSignorovitch et al. (2010)を参考に学習していく.
薬剤Aと薬剤BのRCT(AB試験)と薬剤Aと薬剤CのRCT(AC試験)の結果があって, BとCの効果の違いを求めたい場合, Aを軸にBとCを間接的に比較することが考えられる. このときAはBとCを繋いでおり, アンカー(anchor)であると言う. 「AとBの効果の差」と「AとCの効果の差」で単純に差分を取ると, 効果修飾因子がAB試験とAC試験で分布が異なる場合, バイアスが生じてしまう(効果修飾については昔作成したスライドのフェアな比較を崩すもの ～交絡と効果修飾～ / Confounding EffectModification - Speaker Deckが参考になるかもしれない). このとき, AB試験の個々の患者データがある場合, AC試験側の効果修飾因子の分布になるように, AB試験の個々のデータをうまいこと重み付けすることで, AC試験側の集団を疑似的に再現して, 比較することができる. これがMAICの発想である.
このような状況は実際に起きうる話である. 薬剤Bと薬剤Cの効果の差を示したい薬剤Bのメーカーが, BとCのRCTは実施していないけれど, 薬剤Aと薬剤Bを比較したAB試験については自社のデータとして持っていて, AC試験については論文の公表値で入手できる状況で, AB試験とAC試験では効果修飾因子となる患者背景が異なると考えられる場合, MAIC(もしくはSTC)の使いどころというわけである. ちなみに, AB試験とAC試験でAがアンカーになる場合を例に説明したが, アンカーがないケースでもMAICはできる. ただし, 必要な仮定が厳しくなる. アンカーがないケースは今回は考えないことにする.
TSD18のAppendix DはAnchoredの研究のケースを想定したシミュレーションをしてMAICをしている. このシミュレーションモデルについて説明する. AB試験とAC試験の比較をするものとし. AB試験は薬剤Aと薬剤BのRCT試験, AC試験は薬剤Aと薬剤CのRCT試験に相当するものとする. AB試験はIPD(個人レベルのデータ)があり, AC試験はAgD(集団レベルのデータ)のみがある. Aをアンカーとして、BとCの効果を比較するのが今回の目的となっている.
AB試験とAC試験でNは, 500と300, 年齢の範囲は, 45～75と45～55(IPDがあるAB試験がACの範囲を包括していることに注意), 女性比率は, 64%と80%としている.
個人iの治療tのときのアウトカムは以下のロジットモデルで生成されると考える.
ただし, Aの効果を基準としたときのBとCの効果量は,
と設定する. 今回のシミュレーションモデルではBよりCの方が効果があるという設定になっている. 以上のモデルでは年齢がeffect modifier(効果修飾因子)となっているため, NMAではバイアスが生じてしまう.年齢を調整する必要があるので, MAICで調整を行うというのが手法選択の理由である.
シミュレーションデータの確認
Appendix DのRコード(https://www.sheffield.ac.uk/nice-dsu/tsds/population-adjusted)のデータ生成の部分を実行したものとして, データを確認する.
AB試験のデータを6行(こちらはIPDを持っている設定のため, アクセス可能なデータ).
AB.IPD %>% head()
## # A tibble: 6 x 5
## ID age gender trt y
## <chr> <int> <fct> <chr> <int>
## 1 001 45 Male A 1
## 2 002 71 Female A 1
## 3 003 58 Female A 1
## 4 004 48 Female A 1
## 5 005 69 Male A 1
## 6 006 48 Female A 1
AC試験のデータを6行(こちらはAgDのみの設定なので, アクセス不可能なデータで, シミュレーションのため発生している. 推定の際には個々のデータは使えない).
AC.IPD %>% head()
## # A tibble: 6 x 5
## ID age gender trt y
## <chr> <int> <fct> <chr> <int>
## 1 001 46 Female A 0
## 2 002 52 Female A 1
## 3 003 51 Female A 1
## 4 004 55 Female A 1
## 5 005 53 Male A 1
## 6 006 51 Female A 1
AB試験のIPDデータの要約統計量は
AB.IPD %>% group_by(trt) %>%
summarise(n(), mean(age), sd(age), `n(male)`=sum(gender=="Male"),
`%(male)`=mean(gender=="Male"),sum(y),mean(y))
## # A tibble: 2 x 8
## trt `n()` `mean(age)` `sd(age)` `n(male)` `%(male)` `sum(y)` `mean(y)`
## <chr> <int> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 A 250 60.5 9.13 88 0.352 216 0.864
## 2 B 250 60.1 9.07 91 0.364 39 0.156
AB試験のIPDデータの要約統計量は
AC.IPD %>% group_by(trt) %>%
summarise(n(), mean(age), sd(age), `n(male)`=sum(gender=="Male"),
`%(male)`=mean(gender=="Male"),sum(y),mean(y))
## # A tibble: 2 x 8
## trt `n()` `mean(age)` `sd(age)` `n(male)` `%(male)` `sum(y)` `mean(y)`
## <chr> <int> <dbl> <dbl> <int> <dbl> <int> <dbl>
## 1 A 150 50.1 3.29 36 0.24 120 0.8
## 2 C 150 49.7 2.96 30 0.2 17 0.113
各試験においては, ランダム化しているので共変量は群間でうまくバランスしているが, AB試験とAC試験では平均年齢と男性割合が異なっていることが分かる(要約統計量はどちらも確認できる. ただし, 論文のTable1にある変数に限られる).
今回, AC試験について使用する集計値データは以下のようになる.
AC.AgD
## age.mean age.sd N.male prop.male y.A.sum y.A.bar N.A y.C.sum y.C.bar
## 1 49.9 3.128785 66 0.22 120 0.8 150 17 0.1133333
## N.C
## 1 150
MAICの重みの推定
調整したい共変量(ここではage)をバランスさせるような重みを考える.
以下のように, 共変量を条件付けたときの各試験に組み込まれる確率を考えて, その確率をロジットモデルで仮定すると,
となる.
となっているのは, 計算を簡略化するために, IPDのある試験の共変量から,集計値のみの試験の共変量の平均値を引いたものを出している(両方の分布を平行移動させているイメージ).
上のウェイトで加重平均を出したものが0に(平行移動させた後の共変量の平均, つまり0に)一致するようにモーメントが釣り合うような条件式は
となる。
上の式から
が導ける. この式を満たすようなパラメータα1を求めたい.
この解は,
を最小化するパラメータα1に等しく, 一意であり一致推定量となることが示されている. これは、
の右辺が
の一階条件であり, 二階微分した
は正定値であり,
を最小化するαが上の方程式の解に等しく, また一意であるため.
要点をまとめると
モーメントが釣り合うようなウェイトのパラメータαを求める.
モーメント条件の方程式が出て, 平行移動しても同じなので, 平行移動して
「αの関数=0」
の形の方程式にして, それの解wを求める.
expを足し合わせる関数になっていて, 数値計算の都合上と思われる(たぶん).
方程式の解を求めるのではなく, その方程式の右辺の関数は
の関数の一階条件となっているうえ, その関数の二階微分したもの(ヘッセ行列)が正定値行列になっており, その関数を一意に最小にするようなαは, 一階条件の解, すなわち上の方程式の解となっている. 方程式の解を求めるのではなく, 最小化の形に持っていく等価な問題に置き換えているということである.
したがって, これを最小化する.
objfn <- function(a1, X){
sum(exp(X %*% a1))
}
# gradient関数を求めているのは最適化を高速化させるため.
gradfn <- function(a1, X){
colSums(sweep(X,1,exp(X%*%a1),"*"))
}
共変量を平行移動したもの.
X.EM.0 <- sweep(with(AB.IPD, cbind(age, age^2)),2,
with(AC.AgD, c(age.mean, age.mean^2+age.sd^2)),"-")
最適化の実行.
opt1 <- optim(par=c(0,0), fn=objfn, gr=gradfn, X=X.EM.0, method = "BFGS")
opt1
## $par
## [1] 3.32289873 -0.03404291
##
## $value
## [1] 192.6877
##
## $counts
## function gradient
## 48 11
##
## $convergence
## [1] 0
##
## $message
## NULL
# gradientなくても使える。速度は多少遅くなるが.
# optim(par=c(0,0), fn=objfn, X=X.EM.0, method = "BFGS")
推定されたパラメータからウェイトを求める. ここでは, ウェイトが何人分にあたるのかを理解しやすいようにリスケールしている.
a1 <- opt1$par
wt <- exp(X.EM.0 %*% a1)
wt.rs <- wt/sum(wt) * N_AB
ウェイトの要約統計量とヒストグラム.
## V1
## Min. :0.000000
## 1st Qu.:0.000011
## Median :0.052907
## Mean :1.000000
## 3rd Qu.:2.071826
## Max. :3.767246
近似的な有効サンプルは
sum(wt)^2/sum(wt^2)
## [1] 173.845
となっている.
重み付けしたABの年齢分布(平均と標準偏差)とACの年齢分布(平均と標準偏差)がバランスしていることが確認できる(上がウェイト付けしたAB試験の年齢の平均と標準偏差, 下がAC試験の年齢の平均と標準偏差).
## # A tibble: 1 x 2
## age.mean age.sd
## <dbl> <dbl>
## 1 49.9 3.14
## age.mean age.sd
## 1 49.9 3.128785
ウェイト付けデータによる推定結果
AB試験のデータにウェイトをつけてロジスティック回帰. AC試験の集団と年齢の分布がバランスした分布で行っていることになる.
fit1 <- AB.IPD %>% mutate(y0 = 1-y, wt=wt) %>%
glm(cbind(y,y0) ~ trt, data = ., family = binomial, weights = wt)
fit1
##
## Call: glm(formula = cbind(y, y0) ~ trt, family = binomial, data = .,
## weights = wt)
##
## Coefficients:
## (Intercept) trtB
## 1.331 -2.687
##
## Degrees of Freedom: 499 Total (i.e. Null); 498 Residual
## Null Deviance: 267
## Residual Deviance: 196.5 AIC: 168.2
サンドウィッチ分散を出す.
V.sw <- vcovHC(fit1)
d.AB.MAIC <- coef(fit1)["trtB"] #Bの効果(係数が対数オッズ比)
var.d.AB.MAIC <- V.sw["trtB","trtB"] #Bの効果の分散
# Estimated log OR of C vs. A from the AC trial
d.AC <- with(AC.AgD, log(y.C.sum * (N.A - y.A.sum) / (y.A.sum * (N.C - y.C.sum))))
var.d.AC <- with(AC.AgD, 1/y.A.sum + 1/(N.A - y.A.sum) + 1/y.C.sum + 1/(N.C - y.C.sum))
ACの対数オッズ比からABの対数オッズ比を引いて, BCの効果を求める.
# Indirect comparison
print(d.BC.MAIC <- d.AC - d.AB.MAIC)
## trtB
## -0.756381
その分散を求める.
print(var.d.BC.MAIC <- var.d.AC + var.d.AB.MAIC)
## [1] 0.2552748
真の値(βの差分で-0.4)とMAICで求めた値と未調整で求めた値(95%信頼区間付き)を図示すると 
となり, 未調整で推定した場合よりはMAICをした方がいい！！とは今回の結果からは必ずしも言えない. サンプルサイズを増やしてみると良い感じというのが実感できた(ここには結果を載せていない). また, シミュレーションを1回でなく, たくさん回して真の値のカバー率とかを確認してみるのが良さそう(こっちはコードを書いていない).
参考文献
以下の文献が参考になる.
・NICEのTSD18
https://www.sheffield.ac.uk/nice-dsu/tsds/population-adjusted
・Signorovitch et al. (2010)
"Comparative effectiveness without head-to-head trials: a method for matching-adjusted indirect comparisons applied to psoriasis treatment with adalimumab or etanercept"
https://pubmed.ncbi.nlm.nih.gov/20831302/
・製薬協の資料
間接比較において母集団調整法は有用か？
Matching Adjusted Indirect Comparison及びSimulated Treatment Comparison
https://www.jpma.or.jp/information/evaluation/results/allotment/indirect_comparison.html
【学習予定】
MAICの理解がまだ浅いので, もう少し深めていきたい. Population-adjusted indirect comparisonsには他にも手法があるので, そちらも学んでいきたいところ.
