傾向スコアの各手法をPythonで実装

はじめに

傾向スコアとは主に観察研究でよく使われる統計解析の手法です。

観察研究では、共変量によるバイアスを小さくするランダム化のような操作を行うことができないため、バイアスを小さくするために傾向スコアが使われます。前回傾向スコアの各手法についてまとめた記事を書きました。

今回はpythonで実装を行い、結果の比較を行っていきます。

▼前回の記事はこちら

データセット

傾向スコアでよく用いられる、Right Heart Catheterization Datasetというデータセットを使用します。

このデータセットでは、心臓カテーテル（RHC）を受けた患者は処置群に、受けてない患者はコントロール群に割り当てられます。また、アウトカムは入院後180日以内に死亡したかどうかです。

今回データセットに含まれる5735人中、RHCを受けた患者は2184人、受けてない患者は3551人です。また、RHCを受けた患者の死亡率は68.0%、受けてない患者の死亡率は63.0%となっています。
ここから傾向スコアの各手法による補正を行い、どの程度結果に差が出るのかをみていきます。

データ前処理

下記の記事を参考にデータの前処理を行います。今回、傾向スコアをはじめ、以降出てくる予測モデルにはロジスティクス回帰を利用します。

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statistics import stdev

data_df = pd.read_csv("rhc.csv")
cols = ["cat1", "sex", "race", "edu", "income",
        "resp", "card", "neuro", "gastr", "renal", "meta", "hema", "seps", "trauma", "ortho",
        "das2d3pc", "dnr1", "ca", "surv2md1", "aps1", "scoma1", "wtkilo1", "temp1",
        "resp1", "hrt1", "pafi1", "paco21", "ph1", "wblc1", "hema1", "sod1", "pot1", "crea1",
        "bili1", "alb1", "cardiohx", "chfhx", "dementhx", "psychhx", "chrpulhx", "renalhx",
        "liverhx", "gibledhx", "immunhx", "transhx", "amihx",
        "age", "meanbp1"]
categorical_columns = ["cat1", "sex", "race", "edu", "income", "ca", "dnr1",
                       "resp", "card", "neuro", "gastr", "renal", "meta", "hema", "seps", "trauma", "ortho"]
# 交絡因子
X = data_df[cols]
dummy = pd.get_dummies(X[categorical_columns], drop_first=True)
X = pd.concat([X, dummy], axis=1)
X = X.drop(categorical_columns, axis=1)
X.loc[:, "Intercept"] = 1
# 処置群とコントロール群
D = pd.get_dummies(data_df["swang1"])["RHC"]
# アウトカム
Y = pd.get_dummies(data_df["death"])["Yes"]
# 傾向スコア
model = sm.Logit(D, X).fit(method="newton")
e = model.predict(X)

実装

マッチングでは、差が標準偏差の1/100以下の場合のみマッチングするように制限を与えています。また、層別化では、0.05区切りで20層に層別化を行いました。

# マッチング
target_list = e[e.index.isin(D[D==1].index)]
ate = 0
k = 0
for i in D[D==0].index:
    num = e[e.index.isin(D[D==0].index)][i]
    idx = np.abs(np.asarray(target_list) - num).argmin()
    if np.abs(target_list.iloc[idx] - num) > stdev(e)/100:
        continue
    ate += int(Y.loc[target_list.index[idx]]) - int(Y.loc[i])
    k += 1
ate = ate / k

# 層別化
interval = np.arange(0,1.05,0.05)
ate = 0
for i in range(0,len(interval)-1):
    temp0 = e[(e.index.isin(D[D==0].index)) & 
              (e>=interval[i]) &
              (e<interval[i+1])].index
    temp1 = e[(e.index.isin(D[D==1].index)) & 
              (e>=interval[i]) &
              (e<interval[i+1])].index
    if (len(temp0) > 0) & (len(temp1) > 0):
        ate += (Y[temp1].mean() - Y[temp0].mean())*(len(temp1) + len(temp0))
ate = ate / len(e)

# カーネルマッチング
f = sm.Logit(Y[Y.index.isin(D[D==0].index)], e[e.index.isin(D[D==0].index)]).fit(method="newton")
Y_hat = f.predict(e[e.index.isin(D[D==1].index)])
ate = (Y[Y.index.isin(D[D==1].index)] - Y_hat).sum() / len(Y_hat)

# DR
model0 = sm.Logit(Y[Y.index.isin(D[D==0].index)], X[X.index.isin(D[D==0].index)]).fit(method="ncg")
Y_hat0 = model0.predict(X)
model1 = sm.Logit(Y[Y.index.isin(D[D==1].index)], X[X.index.isin(D[D==1].index)]).fit(method="ncg")
Y_hat1 = model1.predict(X)
ate = ((D * Y / e - (1 - (D / e) * Y_hat1)).sum() - ((1-D) * Y / (1-e) - (1 - ((1-D) / (1-e)) * Y_hat0)).sum()) / len(Y)

# OW
ato = (D * (1-e) * Y).sum() / (D * (1-e)).sum() - ((1-D) * e * Y).sum() / ((1-D) * e).sum()

結果

マッチングのATEは+9.6%、層別化+5.1%、カーネルマッチング+2.4%、IPW+5.9%、DR+5.1%、OWのATOは+6.0%と、数値に差はあるものの、全ての手法で、処置群の方が死亡率が高くなるという結果が出ました。

また、カーネルマッチングにおいて傾向スコアからアウトカムを予測するモデルを作成した際のAUCは0.49、DRで利用した交絡因子からアウトカムを予測するモデルのAUCはそれぞれ0.65と0.67でした。
これらの手法は、より精度の高い予測モデルを作成するまでは、実際に利用することは難しそうです。

安定性

ここからは、ランダムにデータをサンプリングして結果がどのように変わるかを見ていきます。5735人の中からランダムに4000人抽出して算出した結果は以下の通りです。

今回利用したデータサイズでは、マッチング系の手法は結果があまり安定しないことが分かりました。
また、極端な傾向スコアを取るデータが少なかったこともあり、IPWとOWの安定性が高くなっています。

まとめ

今回は、Right Heart Catheterization Datasetを利用して、傾向スコアの各手法を実装し、結果の比較を行いました。

最後まで読んでいただき、ありがとうございます。

アポロならではの技術的課題に対する取り組みやプロダクト開発の試行錯誤で得た学びなどを定期的に発信していきます。
少しでも業界へ貢献できれば嬉しいです。今後ともよろしくお願いいたします。