これで無理なら諦めて!世界一やさしいデータ分析教室

オーブンソースデータなどWeb上から入手できるデータを用いて、RとPython両方使って分析した結果を書いていきます

傾向スコア(Propensity score)をRで実践 マッチングとIPWの結果を比較

前回、前編・後編と2回に傾向スコアの考え方について学びました。
www.randpy.tokyo
www.randpy.tokyo

今回は傾向スコアを使って実際に分析をしていきます。
これまでの実践編記事では、主にスクレイピングを使ってデータを集めていました。

今回は趣向を変えて、公開されているオープンソースデータを使って傾向スコアの実践をしていきます。

なお、こちらの記事のPython編も公開しました。
Pythonの方が興味あるぜ!って方はこちらご覧ください。
www.randpy.tokyo

扱うデータの説明

Rにはもともと傾向スコアを計算する用のlalondeというデータセットがあるのですが、それを使った記事が多く見られたので別のデータを使ってみることにしました。

http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets
こちらのサイトに分析に使えそうな多くのデータが掲載されています。

この中から傾向スコアの論文を読んでいるとしょっちゅう見かける、「Right heart catheterization dataset」という心臓カテーテルに関するデータを発見したのでそちらを使います。
(記事執筆後、ググってみたらこちらのデータ使った傾向スコア解説記事も普通にありました…無念……)

なお、いつも通りコードは私のGithub上にあげています。
github.com

まずはデータを読み込みましょう。

dataSet = read.csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")

ここでは直接URLから取得していますが、該当データをクリック&ダウンロードして、それを作業ディレクトリに保存して

dataSet = read.csv("rhc.csv")

としてももちろん問題ないです。

簡単にどんなデータか確認しておきましょう。

> head(dataSet)
  X              cat1          cat2  ca sadmdte dschdte dthdte lstctdte death cardiohx chfhx dementhx psychhx chrpulhx renalhx liverhx gibledhx malighx immunhx transhx amihx      age    sex       edu  surv2md1
1 1              COPD          <NA> Yes   11142   11151     NA    11382    No        0     0        0       0        1       0       0        0       1       0       0     0 70.25098   Male 12.000000 0.6409912
2 2     MOSF w/Sepsis          <NA>  No   11799   11844  11844    11844   Yes        1     1        0       0        0       0       0        0       0       1       1     0 78.17896 Female 12.000000 0.7549996
3 3 MOSF w/Malignancy MOSF w/Sepsis Yes   12083   12143     NA    12400    No        0     0        0       0        0       0       0        0       1       1       0     0 46.09198 Female 14.069916 0.3169999
4 4               ARF          <NA>  No   11146   11183  11183    11182   Yes        0     0        0       0        0       0       0        0       0       1       0     0 75.33197 Female  9.000000 0.4409790
5 5     MOSF w/Sepsis          <NA>  No   12035   12037  12037    12036   Yes        0     0        0       0        0       0       0        0       0       0       0     0 67.90997   Male  9.945259 0.4369998
6 6              COPD          <NA>  No   12389   12396     NA    12590    No        0     1        0       0        1       0       0        0       0       0       0     0 86.07794 Female  8.000000 0.6650000
  das2d3pc t3d30 dth30 aps1 scoma1 meanbp1       wblc1 hrt1 resp1    temp1    pafi1     alb1    hema1     bili1     crea1 sod1     pot1 paco21      ph1 swang1  wtkilo1 dnr1           ninsclas resp card neuro
1 23.50000    30    No   46      0      41 22.09765620  124    10 38.69531  68.0000 3.500000 58.00000 1.0097656 1.1999512  145 4.000000     40 7.359375 No RHC 64.69995   No           Medicare  Yes  Yes    No
2 14.75195    30    No   50      0      63 28.89843750  137    38 38.89844 218.3125 2.599609 32.50000 0.6999512 0.5999756  137 3.299805     34 7.329102    RHC 45.69998   No Private & Medicare   No   No    No
3 18.13672    30    No   82      0      57  0.04999542  130    40 36.39844 275.5000 3.500000 21.09766 1.0097656 2.5996094  146 2.899902     16 7.359375    RHC  0.00000   No            Private   No  Yes    No
4 22.92969    30    No   48      0      55 23.29687500   58    26 35.79688 156.6562 3.500000 26.29688 0.3999634 1.6999512  117 5.799805     30 7.459961 No RHC 54.59998   No Private & Medicare  Yes   No    No
5 21.05078     2   Yes   72     41      65 29.69921880  125    27 34.79688 478.0000 3.500000 24.00000 1.0097656 3.5996094  126 5.799805     17 7.229492    RHC 78.39996  Yes           Medicare   No  Yes    No
6 17.50000    30    No   38      0     115 18.00000000  134    36 39.19531 184.1875 3.099609 30.50000 1.0097656 1.3999023  138 5.399414     68 7.299805 No RHC 54.89999   No           Medicare  Yes   No    No
  gastr renal meta hema seps trauma ortho adld3p urin1  race     income ptid
1    No    No   No   No   No     No    No      0    NA white Under $11k    5
2    No    No   No   No  Yes     No    No     NA  1437 white Under $11k    7
3    No    No   No   No   No     No    No     NA   599 white   $25-$50k    9
4    No    No   No   No   No     No    No     NA    NA white   $11-$25k   10
5    No    No   No   No   No     No    No     NA    64 white Under $11k   11
6    No    No   No   No   No     No    No      0   242 white Under $11k   12

たくさん変数がありますね。ここでまず重要なのが「swang1」という変数です。これは、右心カテーテル(rhc : right heart catheterization)の施術を受けたかどうかを表しています。
そして「death」が亡くなってしまったかどうかを表しています。

あとは年齢だったり性別だったり、教育年数だったり、ガンなどその他の病気にかかっているか、といったデータが含まれています。
その他詳しい変数についての説明は変数説明をご覧ください。かなりざっくりと書いてあり、あとたまに変数名が異なることがあり混乱するかもしれませんが、まあ分かると思います。

データセットを眺めていると、施術が行われたか(swang1)と死亡率(death)を元に分析ができそうですね。
通常は施術が行われた方が死亡率は改善しそうですが、そもそも施術が行われるということは、回復する見込みがあるということかもしれません。もしくは回復見込みがどうしても無いから、最終手段として手術を行っているのかもしれません。

その場合、施術が行われた人は回復しやすい人、施術が行われなかった人は回復しづらい人(もしくはその逆)である可能性があるということになります。
そのような属性の偏りがあるまま、死亡率の比較をしても意味がありません。

そこで傾向スコアの出番!!というわけですね。
今回は、年齢や性別、他の病気にかかっているか、という変数をもとに施術が行われる確率(= 傾向スコア)を求めていきます。

傾向スコアとロジスティック回帰で求める

その前に、傾向スコアを導入する前の、施術と死亡率の関係を確認しておきましょう。

> table(dataSet[,c("swang1", "death")])
      death
swang1   No  Yes
     0 1315 2236
     1  698 1486

table関数で2変数の関係をまとめたものを取得しました。施術が行われた(swang1 == 1)ときと行われなかった(swang == 0)ときで死亡率を比較すると、

> 1486/(698+1486) - 2236/(1315+2236)
[1] 0.05072115

施術がおこなれたほうが5%ほど死亡率が高いという結果になりました。

この結果が傾向スコアを使って調整すると、どうなるのか?
一緒に分析を進めていきましょう。

まず、ロジスティック回帰を使って施術が行われそうな確率を求めます。

propensityScoreForAllVal <- glm(swang1 ~ age + sex + race + edu + income + ninsclas + cat1 + das2d3pc + dnr1 + ca + surv2md1 + aps1 + scoma1 + wtkilo1 + temp1 + meanbp1 + resp1 + hrt1 + pafi1 + paco21 + ph1 + wblc1 + hema1 + sod1 + pot1 + crea1 + bili1 + alb1 + resp + card + neuro + gastr + renal + meta + hema + seps + trauma + ortho + cardiohx + chfhx + dementhx + psychhx + chrpulhx + renalhx + liverhx + gibledhx + malighx + immunhx + transhx + amihx,
                            family  = binomial(link = "logit"),
                            data    = dataSet)
step(propensityScoreForAllVal)

とりあえず全変数を使っています。
その後、step関数を使って
www.randpy.tokyo
こちらの際に軽く解説した、AIC基準に基づいた最適な変数選択をしてみます。

そうして選ばれた変数を使ってもう一度ロジスティック回帰。

> propensityScoreModel = glm(swang1 ~ age + edu + ninsclas + cat1 + dnr1 + ca + surv2md1 + aps1 + scoma1 + wtkilo1 + meanbp1 + resp1 + hrt1 + pafi1 + paco21 + ph1 + hema1 + sod1 + pot1 + crea1 + alb1 + resp + card + neuro + gastr + renal + hema + seps + trauma + dementhx + psychhx + renalhx + transhx, 
                       family = binomial(link = "logit"), 
                       data = dataSet)

結果は以下のようになると思います。

> summary(propensityScoreModel)

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                 15.9413914  2.9914119   5.329 9.87e-08 ***
age                         -0.0040670  0.0027626  -1.472 0.140986    
edu                          0.0243429  0.0104295   2.334 0.019593 *  
ninsclasMedicare             0.2283698  0.1338679   1.706 0.088020 .  
ninsclasMedicare & Medicaid  0.4042857  0.1681289   2.405 0.016189 *  
ninsclasNo insurance         0.5227892  0.1657751   3.154 0.001613 ** 
ninsclasPrivate              0.4078453  0.1181147   3.453 0.000554 ***
ninsclasPrivate & Medicare   0.3393003  0.1367745   2.481 0.013111 *  
cat1CHF                      0.9208619  0.1378854   6.678 2.41e-11 ***
cat1Cirrhosis               -1.2241020  0.1987718  -6.158 7.35e-10 ***
cat1Colon Cancer            -0.1060948  1.1196155  -0.095 0.924505    
cat1Coma                    -0.5853801  0.1803182  -3.246 0.001169 ** 
cat1COPD                    -0.4264449  0.1744703  -2.444 0.014516 *  
cat1Lung Cancer             -0.3822752  0.5410519  -0.707 0.479852    
cat1MOSF w/Malignancy       -0.0134676  0.1584761  -0.085 0.932276    
cat1MOSF w/Sepsis            0.5286990  0.0901422   5.865 4.49e-09 ***
dnr1Yes                     -0.6694645  0.1155871  -5.792 6.96e-09 ***
caNo                         0.8059294  0.1722874   4.678 2.90e-06 ***
caYes                        0.2585397  0.1593823   1.622 0.104774    
surv2md1                    -1.9072657  0.3495673  -5.456 4.87e-08 ***
aps1                         0.0052484  0.0028034   1.872 0.061187 .  
scoma1                      -0.0043847  0.0015843  -2.768 0.005648 ** 
wtkilo1                      0.0062375  0.0011927   5.230 1.70e-07 ***
meanbp1                     -0.0061301  0.0010032  -6.111 9.93e-10 ***
resp1                       -0.0215001  0.0025674  -8.374  < 2e-16 ***
hrt1                         0.0051457  0.0008781   5.860 4.62e-09 ***
pafi1                       -0.0046667  0.0003400 -13.725  < 2e-16 ***
paco21                      -0.0257867  0.0036083  -7.147 8.90e-13 ***
ph1                         -1.6180934  0.3762769  -4.300 1.71e-05 ***
hema1                       -0.0112413  0.0046567  -2.414 0.015778 *  
sod1                        -0.0111386  0.0042674  -2.610 0.009050 ** 
pot1                        -0.1653928  0.0337861  -4.895 9.82e-07 ***
crea1                        0.0556529  0.0216191   2.574 0.010046 *  
alb1                        -0.0907288  0.0471720  -1.923 0.054435 .  
respYes                     -0.2838936  0.0818782  -3.467 0.000526 ***
cardYes                      0.5747452  0.0841161   6.833 8.33e-12 ***
neuroYes                    -0.5049139  0.1335887  -3.780 0.000157 ***
gastrYes                     0.3494236  0.1009571   3.461 0.000538 ***
renalYes                     0.3016878  0.1483822   2.033 0.042034 *  
hemaYes                     -0.5279391  0.1460169  -3.616 0.000300 ***
sepsYes                      0.2752581  0.0913221   3.014 0.002577 ** 
traumaYes                    1.2484435  0.3341051   3.737 0.000186 ***
dementhx                    -0.3978487  0.1197277  -3.323 0.000891 ***
psychhx                     -0.3927508  0.1370921  -2.865 0.004172 ** 
renalhx                     -0.3401781  0.1793038  -1.897 0.057799 .  
transhx                      0.4777575  0.0984829   4.851 1.23e-06 ***

傾向スコアを表示したい場合は、

propensityScores = propensityScoreModel$fitted.values

という形で、fitted.valuesで見ることができます。

傾向スコアマッチングを実践

傾向スコアマッチングについて、傾向スコア解説後編で簡単に説明しましたが、実践するにあたってもう少し追記しておきます。

傾向スコアマッチングのアルゴリズムとは、施術を受けた群からランダムに1つの被験者を取得し、その傾向スコアと最も近い値を持つ被験者を施術を受けなかった群から探してペアにするという考えです。どんどん作って、ペア同士で値(ここでは死亡したかどうか)の差を計算するという操作を全ペアで行い、その差の期待値を最終的に2群の差とするというものです。

傾向スコアマッチングの中には、ペアを選択する際に、1度別のペアで選択されたデータは外すという考え(非復元抽出)と、1度別のペアで選択されたデータも候補とするという考え(復元抽出)という2つの手法があります。またさらに、最も近い傾向スコアを持つものを選ぶといっても、離れすぎているものはNGにしましょう、ということで例えば傾向スコアが0.2離れたらペアを作らない、といった制約をかけることもあります。他にも1:1でペアを作るのではなく、1:多とすることもあります(これについてはあまりよく分かっていません…)。

今回は、非復元抽出、1:1、0.1より離れた傾向スコアはペアとしない、という条件でマッチングを行います。

傾向スコアマッチングの結果とコード

#install.packages("Matching")
library(Matching) 
propensityScoreMatching0.1 = Match(Y = as.integer(dataSet$death)-1 , Tr = (dataSet$swang1==1), X = propensityScores, M=1,caliper = 0.1, ties=FALSE, replace = FALSE)

マッチング自体はMachingというライブラリを使って計算を行います。
Yには傾向スコア調整後に比較したい値、Trには何をもって施術された群とするか、Xには各被験者の傾向スコアが入ります。

結果を確認すると、

> summary(propensityScoreMatching0.1)

Estimate...  0.059053 
SE.........  0.016992 
T-stat.....  3.4753 
p.val......  0.0005102 

Original number of observations..............  5735 
Original number of treated obs...............  2184 
Matched number of observations...............  1541 
Matched number of observations  (unweighted).  1541 

Caliper (SDs)........................................   0.1 
Number of obs dropped by 'exact' or 'caliper'  643 

となりました。
施術が行われた群と行われなかった群の死亡率の差を見ると、5.9%。

単純に比較した値よりも死亡率の差が大きくでました。
このことから、傾向スコア調整後も、やはりRHCという施術には効果がないのでは?むしろ悪化させるのでは?ということが言えそうです。

IPWを実践

次にIPW(傾向スコアによる重み付け推定法)もやってみます。

傾向スコアによる重み付け推定法では、傾向スコアの逆数によって重みづけをしてあげることで、各群の期待値を直接求めることができるという便利な手法です。
(傾向スコアマッチングでは、2群の差の期待値は求まるが各群の期待値を求めることはできない。)

IPWの結果とコード

Y  <- as.integer(dataSet$death)-1
z1 <- dataSet$swang1
ipwe1 <- sum((z1*Y/propensityScores)/sum(z1/propensityScores))
ipwe0 <- sum(((1-z1)*Y)/(1-propensityScores))/sum((1-z1)/(1-propensityScores))

こんな感じです。とてもシンプルにかけました。(統計的因果推論(2): 傾向スコア(Propensity Score)の初歩をRで実践してみるを参考にさせていただきました。)

> ipwe1
[1] 0.682175
> ipwe0
[1] 0.6249477
> ipwe1 - ipwe0
[1] 0.05722736

施術が行われた群と行われなかった群の死亡率の差を見ると、5.7%。

こちらも単純に比較した値よりも死亡率の差が大きくでました。

2つの傾向スコア調整法から考えると、やはりRHCって効果が薄い?(むしろ死亡率を高めている....)という考えさせられる結果が出てきました。
(もちろん、今回使ったデータセットのバイアスや分析手法の課題などの様々な問題はあり、断定して結論づけられるわけではないので注意してくださいませ。)

まとめ

傾向スコアを使って調整してあげることで、ただ比較するよりも厳密に「右心カテーテル法」の効果を確認してみました。
本来傾向スコアを使う際は、傾向スコアによる調整で本当にランダムなデータと言えるかどうかプロットなどを通して確認する必要がありますが、今回は省略しました。

元論文を読みたいところですが、お金がかかるので諸事情でちょっと断念してしまいましたが、元論文ではどんな結果になっているのでしょうか?気になります…。

なお、今回はオープンソースのデータを使って分析を行いました。
実は当初はスクレイピングでデータを集めようかと思っていましたが、傾向スコアに向くようなデータが無さそう…ということで断念しました。
機会があれば、スクレイピングしたデータを使って実践してみます。

もし、「このサイトのデータ使えるのでは?」的な意見があれば、ツイッターや本ブログのコメントにてもらえると大変嬉しいです。
さて、次回はPython編です!

Rで傾向スコアを実践している Webの記事はちょくちょくありますが、Pythonでは珍しいと思います。
お楽しみに!

なお理論編でも紹介しましたが、傾向スコアについて詳しく深いところまで理解したい、という方はこのあたりが参考になるかと思います。
本ブログを読んだ後に読むと、入りやすいかもしれません。

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)

岩波データサイエンス Vol.3

岩波データサイエンス Vol.3