ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

傾向スコア解析(2):傾向スコアを使った逆確率重み付け法(IPW)で交絡を調整する [R]

傾向スコア(propensity score, PS)を使った交絡調整方法には、

  1. マッチング
  2. 層別化
  3. PSを共変量とした回帰モデル
  4. 逆確率重み付け法

がある。

ここでは「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を受ける傾向スコアを p1とすると、それぞれの患者に対する重みwは次の式で計算される。


w_i = \frac{Z_i}{p_{1i}} + \frac{1-Z_i}{1-p_{1i}}

Z=1のときは第2項が消えて、治療を受ける傾向スコアの逆数になる。 また、Z=0のときは第1項が消えて、治療を受けない確率( = 1-p_1)の逆数になる。

上記は正確に言えばaverage treatment effect (ATE)を推定するときの重みで、average treatment effect on the treated(ATT)を推定したいときは次の式で計算される重みを使う。


w_i = Z_i + (1-Z_i)\frac{p_{1i}}{1-p_{1i}}

ここではATE推定のための重みを計算して、データに加える。重みは上記の式で定義してもよいし、if_else( )を使ってもよい。

d <- d %>%
  mutate(wt = if_else(abcix==1,1/ps,1/(1-ps)))

重みが極端な値の場合の対処方法

傾向スコアが0や1に近い値の人が、本来受けると予想される治療と反対の群に属していると、重みがとてつもなく大きな値になってしまう。とてつもなく大きな重みが割り付けられた人に引っ張られて、結果がなんだかおかしなことになってしまうことがある。そのようなときは、次のような対処方法を試してみる。

  1. PSモデル再考
  2. PSを利用した他の解析に変更
  3. weight truncationを行う
  4. 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は、それぞれの群に割り付けられる確率(=治療群なら p_1, 対照群なら 1-p_1)の平均をかけたもの。確率の平均値(0〜1)をかけるので、元のweightよりも小さくなる。


w_i = c_1 \frac{Z_i}{p_{1i}} + c_0 \frac{1-Z_i}{1-p_{1i}}

ここで c_1は「共変量を考慮しない場合に治療を受ける確率(つまり p_1の平均値)」で、 c_0は「共変量を考慮しない場合に治療を受けない確率(つまり 1-p_1の平均値)(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]]

図1:逆確率をそのままweightにした場合

図2:99パーセンタイル以上をtruncationした場合

図3:stabilizationした場合

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)は


ESS = \frac{(\sum w_i)^2}{\sum w_i^2}

で計算されるもの。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パッケージの解説:

ngreifer.github.io

  • 矢内先生@高知工科大学の計量経済学応用の授業スライド。矢内先生の資料はRの環境準備から懇切丁寧に解説してくださっていてとても助かります。

yukiyanai.github.io

*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.