傾向スコア(propensity score, PS)を使った交絡調整方法には、
- マッチング
- 層別化
- PSを共変量とした回帰モデル
- 逆確率重み付け法
がある。
ここでは「4. 逆確率重み付け法(inverted probability weighting, IPW)」について説明するが、傾向スコアを推定するところまでは前回と同じ。
パッケージとデータの準備 〜 傾向スコアの推定
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
:病変血管数
ロジスティック回帰モデルを使って傾向スコアを推定し、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"))
各症例の重み(weight)を計算する
治療の割り付け変数をZ(0/1)、治療Z=1を受ける傾向スコアをとすると、それぞれの患者に対する重みwは次の式で計算される。
Z=1のときは第2項が消えて、治療を受ける傾向スコアの逆数になる。 また、Z=0のときは第1項が消えて、治療を受けない確率()の逆数になる。
上記は正確に言えばaverage treatment effect (ATE)を推定するときの重みで、average treatment effect on the treated(ATT)を推定したいときは次の式で計算される重みを使う。
ここではATE推定のための重みを計算して、データに加える。重みは上記の式で定義してもよいし、if_else( )を使ってもよい。
d <- d %>% mutate(wt = if_else(abcix==1,1/ps,1/(1-ps)))
重みが極端な値の場合の対処方法
傾向スコアが0や1に近い値の人が、本来受けると予想される治療と反対の群に属していると、重みがとてつもなく大きな値になってしまう。とてつもなく大きな重みが割り付けられた人に引っ張られて、結果がなんだかおかしなことになってしまうことがある。そのようなときは、次のような対処方法を試してみる。
- PSモデル再考
- PSを利用した他の解析に変更
- weight truncationを行う
- stabilized weightを使う
3つ目のtruncationは、重みの上限を設定してしまう方法。例えば99パーセンタイルよりも大きなweightは全て99パーセンタイル値と同じにしてしまう。
d <- d %>% mutate(wt_trun = if_else(wt>quantile(wt,prob=0.99),quantile(wt,prob=0.99),wt))
4つ目のstabilized weightは、それぞれの群に割り付けられる確率(=治療群なら, 対照群なら)の平均をかけたもの。確率の平均値(0〜1)をかけるので、元のweightよりも小さくなる。
ここでは「共変量を考慮しない場合に治療を受ける確率(つまりの平均値)」で、は「共変量を考慮しない場合に治療を受けない確率(つまりの平均値)(sampling weightがある場合はそれを加味した加算平均だが、ここでは割愛)。
d <- d %>% mutate(c1 = mean(ps), c0 = mean(1-ps), wt_stab = if_else(abcix==1,wt*c1,wt*c0))
比較のために図示してみる。
res <- d %>% pivot_longer(cols=c(wt,wt_trun,wt_stab), names_to = "wt_type", values_to = "weight") %>% group_by(wt_type) %>% do(plots = ggplot(data=., aes(x=weight)) + geom_density(aes(colour=factor(abcix))) + scale_x_continuous(limits = c(0,10)) + ggtitle(.$wt_type) ) res[[2]]
WeightItパッケージを使う方法
どんな計算をしているのか確認するためにここまで書いてきたが、WeightItパッケージを使えばall-in-oneでとても便利。
library(WeightIt) wt_out <- weightit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, estimand="ATE", method="ps", stabilize = FALSE)
estimand
:デフォルトは"ATE"。他にもATT (average treatment effect on the treated), ATC (average treatment effect on the control), ATO (average treatment effect in the overlap), ATM (average treatment effect in the matched sample), and ATOS (average treatment effect in the optimal subset)が指定できるみたい(名前だけ言われても困るが)。method
:デフォルトは"ps"で、一般化線形モデルを使って推定した傾向スコアを使って重み付けする。他にもオプションがあるが主に傾向スコアをどういったアルゴリズムで推定するか。stabilize
:前述のstabilized weightを使うかどうか。デフォルトはFALSEになっている。
計算された重みはweights
という名前で保存されているので、wt_out$weights
で取り出して利用可能。
計算されたweightの概要はsummary( )で確認できる。
> summary(wt_out) Summary of weights - Weight ranges: Min Max treated 1.0204 |-| 3.2033 control 1.3027 |--------------------------| 23.9979 - Units with 5 greatest weights by group: 34 80 93 72 22 treated 2.4968 2.5399 2.6684 2.673 3.2033 988 995 978 993 979 control 14.8467 15.1977 15.6449 18.0062 23.9979 - Weight statistics: Coef of Var MAD Entropy # Zeros treated 0.200 0.154 0.019 0 control 0.703 0.368 0.153 0 - Effective Sample Sizes: Control Treated Unweighted 298. 698. Weighted 199.68 671.09
共変量のバランスを確認する
cobaltパッケージを使う
WeightItパッケージのreferenceにもcobaltパッケージのbal.tab( )を使う方法が紹介されている。stats
で表示する統計量を指定できる。ここでは
"m"
:標準化平均差(standardized mean difference)を表示。そのままだと2値変数については単純に割合の差を返すだけなので、2値変数についても標準化したけでばbinary="std"
をつける。"v"
:分散比(variance ratio)を表示。連続変数に対してのみ計算される。
library(cobalt) bal.tab(wt_out, stats = c("m","v"), binary="std")
Call weightit(formula = abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, method = "ps", estimand = "ATE", stabilize = FALSE) Balance Measures Type Diff.Adj V.Ratio.Adj prop.score Distance -0.0192 0.9373 stent Binary 0.0062 . height Contin. -0.0120 0.8409 female Binary 0.0217 . diabetic Binary -0.0517 . acutemi Binary -0.0032 . ejecfrac Contin. -0.0005 0.9475 ves1proc Contin. -0.0736 0.9061 Effective sample sizes Control Treated Unadjusted 298. 698. Adjusted 199.68 671.09
Effective sample size(ESS)は
で計算されるもの。Weightに偏りがある場合、例えば1つだけ他の症例よりも1000倍も重み付けされる症例があると、100例あっても実際は約4例分にしかならない。全く偏りがない(=全て同じweight)ならば、実際のサンプルサイズとESSは一致する。
MatchItパッケージを使う
MatchItパッケージのmatchit( )でもmethod=NULL
としてマッチング方法を指定しないことでバランスをチェックできる。
サンプリング確率を指定するs.weights
に、推定した傾向スコアを指定してやればよい。
library(MatchIt) wt_tab <- matchit(abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data=d, method=NULL, s.weights = wt_out$weights)
> summary(wt_tab) Call: matchit(formula = abcix ~ stent + height + female + diabetic + acutemi + ejecfrac + ves1proc, data = d, method = NULL, s.weights = wt_out$weights) Summary of Balance for All Data: Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max distance 0.4981 0.4964 0.0865 0.9000 0.0187 0.0492 stent 0.6702 0.6673 0.0063 . 0.0029 0.0029 height 171.4044 171.5317 -0.0120 0.8409 0.0087 0.0376 female 0.3438 0.3334 0.0219 . 0.0104 0.0104 diabetic 0.2228 0.2447 -0.0526 . 0.0219 0.0219 acutemi 0.1433 0.1443 -0.0029 . 0.0010 0.0010 ejecfrac 50.9431 50.9481 -0.0005 0.9475 0.0133 0.0422 ves1proc 1.3840 1.4284 -0.0670 0.9061 0.0078 0.0280 Sample Sizes: Control Treated All (ESS) 199.68 671.09 All 298. 698. Matched (ESS) 199.68 671.09 Matched 298. 698. Unmatched 0. 0. Discarded 0. 0.
Std.Mean Diff.やVar.Ratioを見比べてみると、cobaltパッケージを使ったときと概ね同じような値になっていそう。
治療効果を推定する
Weightを組み込んだ一般化線形モデルは基本のglm( )ではなく、surveyパッケージのsvyglm( )を使う。svyシリーズは基本的にロバスト標準誤差推定法を使ってくれるみたい(sandwich推定法の一般形であるHorvitz-Thompson型推定法というものが使われているらしい)。
family=binomial
とすると「整数にならないものに二項分布使ってんじゃねえよ!」って怒られるので、quasi-binomialなるものを指定します。Weightがかかったアウトカムは整数になってないからダメなんだろうと勝手に理解。
library(survey) dsgn_wt <- svydesign(~1, weights = wt_out$weights, data = d) fit <- svyglm(death ~ abcix, design = dsgn_wt, family = quasibinomial(link="logit"))
> summary(fit) Call: svyglm(formula = death ~ abcix, design = dsgn_wt, family = binomial(link = "logit")) Survey design: svydesign(~1, weights = wt_out$weights, data = d) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -2.4250 0.3918 -6.190 8.8e-10 *** abcix -1.7471 0.5031 -3.472 0.000538 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for binomial family taken to be 1.001005) Number of Fisher Scoring iterations: 7
他にはbootstrappingを用いる方法があるが、本記事ではここまで。
おわりに
- マッチングに比べて何をやってるかが直感的に理解しにくい
- 結局のところ実装に関しては、weightをどうするか、ロバスト分散推定をどの関数を使ってやるか、が問題。
- 先住猫(サビ・メス・3歳)予防接種前にオヤツを控えてましたが、体重は増えてました
参考資料
- WeightItパッケージの解説:
- 矢内先生@高知工科大学の計量経済学応用の授業スライド。矢内先生の資料はRの環境準備から懇切丁寧に解説してくださっていてとても助かります。
*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.