Np-Urのデータ分析教室

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

トップ > Pythonで実践編 > Pythonで傾向スコア(Propensity score)マッチングとIPWを実装してみた

Pythonで傾向スコア(Propensity score)マッチングとIPWを実装してみた

Pythonで実践編 傾向スコア

さて、今回は傾向スコアマッチングのPythonによる実践編です。
傾向スコアって何？という方は、まずはこちらの記事を参考にしてみてください。
www.randpy.tokyo
www.randpy.tokyo

今回の趣旨としては、Pythonでの実装という部分に重きを置いていますので、手法の細かいは説明は致しません。
というのも、傾向スコア関連の分析でPythonを使ったものが中々見当たらず、パッケージもほぼ皆無なので（誰か知っていたら教えてくださいー）、まずはPythonでもできるんだぞ！というところをお見せしたく、本記事では実装がメインとなっております。

傾向スコアマッチング、IPWなどの手法についてはRによる実践編で触れております。また、分析に使用するデータ、流れも基本的にはRによる実践編に習い今回進めていきますので、以下記事に一度目を通してから本記事に入っていただくといいかと思います。
www.randpy.tokyo

データセット

Rによる実践編でも使われたデータセットを使用します。
http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/

データの中身はこんな感じ。

#coding:utf-8
import numpy as np
import pandas as pd
import statsmodels.api as sm

#データのインポート
df = pd.read_csv("http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/rhc.csv")
df.head()

Unnamed: 0	cat1	cat2	ca	sadmdte	dschdte	dthdte	lstctdte	death	cardiohx	...	meta	hema	seps	trauma	ortho	adld3p	urin1	race	income	ptid
0	1	COPD	NaN	Yes	11142	11151.0	NaN	11382	No	0	...	No	No	No	No	No	0.0	NaN	white	Under $11k	5
1	2	MOSF w/Sepsis	NaN	No	11799	11844.0	11844.0	11844	Yes	1	...	No	No	Yes	No	No	NaN	1437.0	white	Under $11k	7
2	3	MOSF w/Malignancy	MOSF w/Sepsis	Yes	12083	12143.0	NaN	12400	No	0	...	No	No	No	No	No	NaN	599.0	white	$25-$50k	9
3	4	ARF	NaN	No	11146	11183.0	11183.0	11182	Yes	0	...	No	No	No	No	No	NaN	NaN	white	$11-$25k	10
4	5	MOSF w/Sepsis	NaN	No	12035	12037.0	12037.0	12036	Yes	0	...	No	No	No	No	No	NaN	64.0	white	Under $11k	11
5 rows × 63 columns

このデータは、RHCの施術を受けた患者、受けなかった患者のデータセットになります。
データセットの詳細は以下URLを参考にしてみてください。
Right Heart Catheterization Dataset

今回の目的は、RHCの施術が死亡率に与える影響について分析することですが、単純に施術を受けた患者とそうでない患者の死亡率を比較しても意味がありません。
例えば、RHCを受ける患者というのは、

  • そもそも重症患者で回復する見込みが低い人
  • 施術を受けて回復する見込みが高い人

の２通り考えられ、施術を受ける人というのは何らかのバイアスを受けている可能性が高いです。
なので、前回学んだ傾向スコアを使って、施術を受ける確率が等しい人同士で死亡率を比較していきます。

傾向スコアで比較する前に、単純に施術した人としなかった人の死亡率の差を確認してみます。
クロス集計を求めるにはpandasのcrosstabを使うと便利です。

pd.crosstab(df.death,df.swang1)

swang1	No RHC	RHC
death		
No	1315	698
Yes	2236	1486

1486/(698+1486) - 2236/(1315+2236) = 0.0507となり、施術した人の方が5％ほど死亡率が高いようです。

傾向スコアの算出

前回の都議選の分析でも使用したstatsmodelsライブラリを使用して、ロジスティック回帰を行います。

# データの整形
df = df.iloc[:,1:]
df1 = df[['cat1','sex','death','age','swang1','meanbp1','aps1']]
df1.head()

#傾向スコアを求めるのに使う説明変数
X = df1.drop(['death','swang1'],1)
dummy_X = pd.get_dummies(X[['cat1','sex']],drop_first=True)
X2 = pd.merge(X, dummy_X, left_index=True, right_index=True)
X = X2.drop(['cat1','sex'],1)

#施術したかどうかダミー変数化
ps_target = df1['swang1']
dummy_ps = pd.get_dummies(ps_target)
ps_target = dummy_ps['RHC']

#傾向スコアを求める
glm = sm.Logit(ps_target, X)
result = glm.fit()
ps = pd.Series(result.predict(X))

0    0.176820
1    0.494097
2    0.474606
3    0.368799
4    0.602842
dtype: float64

目的変数は、RHCの有無(データのswang1にあたります）で、説明変数に使ったデータは、性別、年齢、症状、血圧、重症度評価です。
（本当はAICなどで、全ての変数を使ってモデルの評価を行った方がいいとは思いますが、全変数を使ったらなぜか推定できませんでしたので、今回はいくつかの変数をピックアップしております。。。）
カテゴリ変数はダミー変数化してあげて、いつも通りsm.Logitでモデルを定義し、fit()で推定してあげます。
そのあとのpredictメソッドを使うことで施術を受ける確率というのを求めることができます。めでたし、めでたし！

さて、propensity score matchingとIPWの実装に移りますが、今一度理論について確認したい方は前回のRによる実践編の記事を再度読んでいただければと思います。

Propensity score matching

propensity score matchingを実装していきます。
Pythonでは、ライブラリが見当たらなかったので、
http://tomoshige-n.hatenablog.com/
こちらの記事を参考にさせていただきました。（Pythonで書き直した感じです。。。）
施術を受けた人と受けなかった人の傾向スコアを0.05刻みでマッチングさせています。

#death ダミー変数化
y = df1['death']
dummy_y = pd.get_dummies(y)
y = dummy_y['Yes']

ps = pd.Series(result.predict(X))
z1 = ps_target

#Propensity score matching
table = pd.concat([ps,z1,y],axis=1)
table.columns = ['ps','treat','died']

interval = np.arange(0,1.05,0.05)
match_list = []
for i in range(0,len(interval)-1):
    temp0 = table[(table['treat']==0) & (interval[i] < table['ps']) & (table['ps'] < interval[i+1])]
    temp1 = table[(table['treat']==1) & (interval[i] < table['ps']) & (table['ps'] < interval[i+1])]
    if (len(temp0) > 0) & (len(temp1) > 0):
        match_list.append(temp1['died'].mean()-temp0['died'].mean())

np.mean(match_list)

こちら計算してみると0.0402になりました。RHCを受けた人と受けなかった群の死亡率の差は4.02%ほど。
単純に比較した時と比べると差は小さくなりましたが、死亡率はRHCを受けた人の方が高いという結果です。

IPW

重み付け推定法もやってみましょう。こちらの実装は統計的因果推論(2): 傾向スコア(Propensity Score)の初歩をRで実践してみるを参考にさせていただきました。

ps = pd.Series(result.predict(X))
z1 = ps_target
y = y
ipwe1 = sum((z1*y)/ps)/sum(z1/ps)
ipwe0 = sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps))
ipwe1 - ipwe0

ipwe1-ipwe0を計算してみると、0.044となり、マッチングの結果とだいたい同じような結果になりました。
傾向スコアを考慮すると、死亡率の差は小さくなりましたが、依然としてRHCを受けた人の方が死亡率は高くなりました。

ちなみに、Rによる実践編では、傾向スコアを使った方が死亡率の差は大きくなっていました。この辺りは傾向スコアをどんな変数を使って求めているのか、GLMの解法、マッチングの仕方が異なることが原因と思われます。

終わりに

今回は傾向スコアの手法をPythonでもやってみるということで何とか実装できてよかったと思います。
とは言っても、ほぼほぼ偉大な先人の方々のコードをPythonで書き直しただけですが、、笑
ちなみに、PythonでRを使うためのライブラリ(pyper)があるため、無理してPythonで全て書かなくてもいいのでは？と思ったりもしましたが、
そうなるとこのブログの趣旨が少し変わってきてしまうので、できる限りライブラリがない場合もPythonで頑張って書いていこうと思います。
(pyperも便利なので、今度紹介できればと思います！）

次回は機械学習分野の記事を書く予定ですので、お楽しみにー！！