前回、前編・後編と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 1486table関数で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では珍しいと思います。
お楽しみに！
なお理論編でも紹介しましたが、傾向スコアについて詳しく深いところまで理解したい、という方はこのあたりが参考になるかと思います。
本ブログを読んだ後に読むと、入りやすいかもしれません。
調査観察データの統計科学―因果推論・選択バイアス・データ融合 (シリーズ確率と情報の科学)
- 作者: 星野崇宏
- 出版社/メーカー: 岩波書店
- 発売日: 2009/07/29
- メディア: 単行本
- 購入: 29人 クリック: 285回
- この商品を含むブログ (25件) を見る
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2016/06/10
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (4件) を見る
