1
/
5

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

はじめに

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

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

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

▼前回の記事はこちら

傾向スコアの各手法概要|APOLLO
はじめに 傾向スコアとは主に観察研究でよく使われる統計解析の手法です。観察研究では、共変量によるバイアスを小さくするランダム化のような操作を行うことができません。そこで、バイアスを小さくするために使われる手法が傾向スコアです。 今回は傾向スコアについて紹介し、さらに、その中で代表的な手法についてまとめたいと思います。 また、第二弾では、これらの手法を実装し、比較を行う予定です。 傾向スコア 交絡因子$${X_i}$$から予測された処置群$${D=1}$$に割り当てられる確率を傾向スコアと呼び、 $$ P(
https://note.com/apollo132/n/n0b6afb57e542

データセット

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

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

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

データ前処理

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

傾向スコアでセレクションバイアスを補正する - Qiita
はじめに  前回、『Uplift-Modelingで介入効果を最適化する』では、RCTによって得られたデータからある介入の効果を平均的に求めるのではなくて、ある特徴量を持った集団にのみ介入するという戦略を試しました。その結果、「...
https://qiita.com/usaito/items/09daccdd91bc98c21dff
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を利用して、傾向スコアの各手法を実装し、結果の比較を行いました。

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

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

アポロ株式会社's job postings
1 Likes
1 Likes

Weekly ranking

Show other rankings
Invitation from アポロ株式会社
If this story triggered your interest, have a chat with the team?