- パッケージとデータの準備
- マッチング前のデータで共変量のバランスを確認する
- 傾向スコアでマッチさせたコホートを作成する
- マッチさせたコホートで共変量のバランスを確認する
- マッチングの様子を図示する
- 治療効果を推定する
- おわりに
- 参考資料
傾向スコア（propensity score, PS）を使った交絡調整方法には、
- マッチング
- 層別化
- PSを共変量とした回帰モデル
- 逆確率重み付け法
がある。
ここでは「1. マッチング（propensity score matching, PSM）」について説明する。
パッケージとデータの準備
PSAgraphicsパッケージのlindnerデータセットを使う。
lindnerデータは、Lindnerセンターで経皮的冠動脈インターベンション（percutaneous cardiac intervention, PCI）を受けた996例において、abciximabという抗血小板薬を追加したときの効果をみた研究のデータ。*1
下のコードでPCI後6か月時点の予後をdeathとして追加した。
library(tidyverse) library(PSAgraphics) data("lindner") d <- lindner %>% mutate(death = if_else(lifepres==0,1,0))
> head(d) lifepres cardbill abcix stent height female diabetic acutemi ejecfrac ves1proc death 1 0.0 14301 1 0 163 1 1 0 56 1 1 2 11.6 3563 1 0 168 0 0 0 56 1 0 3 11.6 4694 1 0 188 0 0 0 50 1 0 4 11.6 7366 1 0 175 0 1 0 50 1 0 5 11.6 8247 1 0 168 1 0 0 55 1 0 6 11.6 8319 1 0 178 0 0 0 50 1 0
このうち、使用する変数は以下のとおり。
death：PCI後6か月時点の予後（0=生存, 1=死亡）abcix：abciximab使用有無（0=PCI単独, 1=PCI+abciximab）stent：冠動脈ステント留置（0=なし, 1=あり）height：身長（cm）female：性別（0=男性, 1=女性）diabetic：糖尿病有無（0=なし, 1=あり）acutemi：7日前以内の急性心筋梗塞有無（0=なし, 1=あり）ejecfrac：左室収縮率（%）ves1proc：病変血管数
マッチング前のデータで共変量のバランスを確認する
tableoneパッケージを使って要約する方法を紹介したが、ここではMatchItパッケージを使う。「治療群を示す変数 ~ 共変量」の要領で指定する。
ポイントはmethod=NULLを指定すること。本来はマッチングの方法を指定する引数だが、NULLを指定するとマッチングが行われないので、マッチング前のデータに関してバランス評価を行うことができる。
library(MatchIt) unmatched <- matchit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data=d, method=NULL)
> summary(unmatched) Call: matchit(formula = abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, method = NULL) Summary of Balance for All Data: Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max distance 0.7265 0.6406 0.6609 1.1161 0.1714 0.2760 stent 0.7049 0.5839 0.2652 . 0.1210 0.1210 height 171.4427 171.4463 -0.0003 1.0201 0.0079 0.0250 female 0.3309 0.3859 -0.1168 . 0.0550 0.0550 diabetic 0.2049 0.2685 -0.1575 . 0.0636 0.0636 acutemi 0.1791 0.0604 0.3095 . 0.1187 0.1187 ejecfrac 50.4026 52.2886 -0.1810 1.0238 0.0356 0.1138 ves1proc 1.4628 1.2047 0.3654 2.1614 0.0433 0.1884 Sample Sizes: Control Treated All 298 698 Matched 298 698 Unmatched 0 0 Discarded 0 0
均衡化の指標については後ほど。
傾向スコアでマッチさせたコホートを作成する
MatchItパッケージのmatchit( )を使って治療群の各症例に対して対照群から近い症例を選んでペアを作る。近さの基準となる距離として傾向スコアが用いられることがほとんど。 マッチングの方法には色々あるが、Austin先生*2も「復元抽出なしの1:1最近傍マッチング」が一番いいんじゃない？って仰っているようなので、ひとまずこれを使う。
「治療群を示す変数 ~ 共変量」の他には、以下の引数を指定する。
method：マッチングの方法nearest：最近傍マッチング（nearest matching別名）。別名：貪欲マッチング（greedy matching）
distance：前述のとおり、症例間の距離としてPSが用いられることがほとんどなので、この引数でPSの推定方法を指定することができる。"glm"：デフォルトはdistance="glm"でlink="logit"`が指定されているとロジスティック回帰モデルを使ってPSを推定する。- この他のPS推定法として、一般化加法モデル（generalized additive model, GAM）, 一般化ブーストモデル（generalized boosted model, GBM）, Lasso/ridge/elastic net, 決定木, ランダムフォレストなどが使える。
"mahalanobis"：マハラノビス距離
estimand：効果を推定したい対象集団。デフォルトは"ATT"（Average Treatment effect on the Treated）。replace：対照を復元抽出するかどうか。TRUEとすると、1つの対照症例が複数の治療群症例にマッチさせることが可能。ratio：治療群患者1例に対してマッチさせる対照症例数
matched_1 <- matchit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data=d, method="nearest", distance = "glm", link = "logit", estimand = "ATT", replace=FALSE, ratio=1, seed=1234)
症例間の距離が遠すぎる人同士がペアにならないように、ある値以上近い場合にだけペアにすることがある。これは引数caliperで指定できる。例えばcaliper=0.2とすると、症例間の距離の標準偏差（SD）の0.2倍より小さい場合にだけペアにする。
論文ではときどき"Nearest matching with caliper width of 0.2 SD of the logit PS"のように、PSではなくPSのlogitをcaliperの基準にしているのを見かけるが、この場合は次のコードのようにdistanceを指定する必要がある（単にcaliper=0.2とするとPSを基準としたキャリパーになってしまうから）。
まずはglm( )でPS, logit PSを推定してデータセットに加えておく。
ps_model <- glm(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, family = binomial(link="logit")) d <- d %>% mutate(ps_logit = predict(ps_model, type = "link"), ps = predict(ps_model, type = "response"))
上記で作成したps_logitをdistanceに指定する（推定したPS, logit PSを使用するのでコード実行必要）。
std.caliper=TRUEとした場合は「キャリパーをdistanceのSDの何倍にするか」をcaliperで指定する。
std.caliper=FALSEとした場合は、キャリパーを絶対値で指定する。
matched_2 <- matchit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data=d, method="nearest", distance = d$ps_logit, caliper=0.2, std.caliper = TRUE, ratio=1, replace=F, seed=1234)
マッチングの概要を確認すると、次のようにキャリパーが0.156になっている。
> matched_2 A matchit object - method: 1:1 nearest neighbor matching without replacement - distance: User-defined [caliper] - caliper: <distance> (0.156) - number of obs.: 996 (original), 584 (matched) - target estimand: ATT - covariates: stent, height, female, diabetic, acutemi, ejecfrac, ves1proc
> sd(d$ps_logit)*0.2 [1] 0.1556522
logit PSの0.2倍がキャリパーの設定値と一致していることが確認できた。
マッチさせたコホートで共変量のバランスを確認する
マッチングの様子はsummary( )で確認できる。
> summary(matched_1) Call: matchit(formula = abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, method = "nearest", distance = "glm", link = "logit", estimand = "ATT", replace = F, ratio = 1, seed = 1234) Summary of Balance for All Data: Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max distance 0.7265 0.6406 0.6609 1.1161 0.1714 0.2760 stent 0.7049 0.5839 0.2652 . 0.1210 0.1210 height 171.4427 171.4463 -0.0003 1.0201 0.0079 0.0250 female 0.3309 0.3859 -0.1168 . 0.0550 0.0550 diabetic 0.2049 0.2685 -0.1575 . 0.0636 0.0636 acutemi 0.1791 0.0604 0.3095 . 0.1187 0.1187 ejecfrac 50.4026 52.2886 -0.1810 1.0238 0.0356 0.1138 ves1proc 1.4628 1.2047 0.3654 2.1614 0.0433 0.1884 Summary of Balance for Matched Data: Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max Std. Pair Dist. distance 0.8517 0.6406 1.6242 0.2455 0.4407 0.8289 1.6242 stent 0.7852 0.5839 0.4414 . 0.2013 0.2013 0.9123 height 171.4295 171.4463 -0.0016 1.1500 0.0109 0.0403 1.0803 female 0.2919 0.3859 -0.1997 . 0.0940 0.0940 0.8986 diabetic 0.1644 0.2685 -0.2577 . 0.1040 0.1040 0.7400 acutemi 0.4128 0.0604 0.9190 . 0.3523 0.3523 1.0065 ejecfrac 47.7953 52.2886 -0.4313 1.2340 0.0800 0.2315 1.0980 ves1proc 1.9228 1.2047 1.0170 2.7619 0.1197 0.5134 1.1405 Percent Balance Improvement: Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max distance -145.7 -1179.0 -157.1 -200.3 stent -66.4 . -66.4 -66.4 height -364.1 -604.0 -37.5 -61.3 female -71.0 . -71.0 -71.0 diabetic -63.6 . -63.6 -63.6 acutemi -196.9 . -196.9 -196.9 ejecfrac -138.2 -795.3 -124.7 -103.4 ves1proc -178.3 -31.8 -176.3 -172.5 Sample Sizes: Control Treated All 298 698 Matched 298 298 Unmatched 0 400 Discarded 0 0
以下の4つのパートで構成されている。
- Summary of Balance for All Data：マッチング前のデータのバランス
- Summary of Balance for Matched Data：マッチング後のデータのバランス
- Percent Balance Improvement：マッチング前後でバランスの各指標が何パーセント改善したか
- Sample Sizes：matched cohortに採用された症例数
バランスの指標として表示されているものは以下のとおり。標準化平均差（standardized mean difference, SMD）が一番よく使われる。経験的累積分布関数（empirical cumulative distribution function, eCDF）に関する指標については明確な基準がない模様。
Means Treated：治療群の平均Means Control：対照群の平均Std. Mean Diff.：治療群間のSMD。<0.1〜0.15ならOK。Var. Ratio：治療群間での分散比。連続変数に対してのみ計算される。上手くバランスがとれていれば、分散比は1に近くなるはず（0.5-2の間なら許容）。eCDF Mean：群間でeCDFが平均してどれくらい乖離しているかeCDF Max：群間でeCDFが最大でどれくらい乖離しているかStd. Pair Dist.：各ペア内の距離
マッチングの様子を図示する
治療群・対照群のそれぞれで、マッチングサンプルに抽出された人、されなかった人の傾向スコアを図示する。
plot(matched, type="jitter")
引数typeを指定して他のグラフを描ける。例えば、type="density"とすれば、各共変量の密度分布を比較する。
plot(matched_1, type="density")
QQプロットやeCDFは連続変数に対してのみ作図されるので、
plot(matched_1, type="qq", which.xs=c("height","ejecfrac","ves1proc"))
plot(matched_1, type="ecdf", which.xs=c("height","ejecfrac","ves1proc"))
他にもcobaltパッケージのlove.plot( )を使ってSMDをプロットすることができる。
library(cobalt) love.plot(matched_1, threshold = 0.15, abs = TRUE, stars="raw")
治療効果を推定する
作成されたマッチング後データは、match.data( )で抽出する。
m.out_1 <- match.data(matched_1)
中身を見てみると、
> head(m.out_1) lifepres cardbill abcix stent height female diabetic acutemi ejecfrac ves1proc death ps_logit ps distance weights subclass 367 0.0 14990 1 1 157 1 1 0 45 2 1 1.215812 0.7713257 0.7713257 1 1 369 0.0 96741 1 1 171 0 0 0 29 1 1 1.242675 0.7760293 0.7760293 1 2 370 11.6 7616 1 1 168 1 0 0 51 2 0 1.364860 0.7965484 0.7965484 1 3 371 11.6 7664 1 0 173 0 0 1 55 1 0 1.453961 0.8106073 0.8106073 1 4 374 11.6 8321 1 0 180 0 0 1 60 1 0 1.272454 0.7811625 0.7811625 1 5 375 11.6 8389 1 0 170 1 0 1 55 1 0 1.141000 0.7578631 0.7578631 1 6
distance, weights, subclassという変数が追加されている。 subclassはペアを示すID。
このデータを使ってglm( )を行うと、次のようになる（今回のマッチングではweightsは全て1なので、指定しなくても結果は同じ）。
fit <- glm(death ~ abcix, data=m.out_1, family = binomial(link = "logit"), weights = weights)
> summary(fit)$coef Estimate Std. Error z value Pr(>|z|) (Intercept) -2.9373967 0.2649533 -11.086471 1.459321e-28 abcix -0.7900164 0.4652834 -1.697925 8.952185e-02
マッチペアを考慮してロバスト分散推定を行うことが勧められているようなので、limtestパッケージおよびsandwichパッケージを使う。
llibrary(lmtest) library(sandwich) coeftest(fit, vcov. = vcovCL, cluster = ~subclass)
Estimate Std. Error z value Pr(>|z|) (Intercept) -2.93740 0.26540 -11.0679 < 2e-16 *** abcix -0.79002 0.45153 -1.7497 0.08018 .
ロバスト分散推定を使わない場合と少しだけ推定値・標準誤差が異なっている。
他にはbootstrappingを用いる方法があるが、本記事ではここまで。
おわりに
- matchit( )には色々なPS推定方法が指定できるので、最終的にマッチングを使わない場合でもmatchit( )でPS推定してしまうのもアリかもしれない。
参考資料
MatchItパッケージの解説：Matching Methods cran.r-project.org
マハラノビス距離が気になる人用 ja.wikipedia.org
*1:Kereiakes DJ, Obenchain RL, Barber BL, Smith A, McDonald M, Broderick TM, Runyon JP, Shimshak TM, Schneider JF, Hattemer CR, Roth EM, Whang DD, Cocks D, Abbottsmith CW. Abciximab provides cost-effective survival advantage in high-volume interventional practice. Am Heart J. 2000 Oct;140(4):603-10.
*2:Peter C. Austin. A comparison of 12 algorithms for matching on the propensity score. Statist. Med. 2014, 33 1057–1069.
