倾向性评分匹配（Propensity Score Matchi

作者: 数科每日 | 来源:发表于2022-01-28 12:10 被阅读0次

本文参考自 Propensity Score Matching，Zolzaya Luvsandorj

Propensity Score Matching 是观察性因果研究中最常用的一种 Matching Sample 的手段，本文通过一个例子来介绍这种方法。

术语

在开始之前，我们先介绍三个术语

结果变量 Outcome variable: 代表实验结果的变量。
干预变量 Treatment variable: 我们希望通过控制，研究的变量, 是实验要研究的 “因”。
混淆变量 Confounding variable: 同时影响 Outcome variable 和Treatment variable 的变量，是实验的干扰因素。

image.png

观察数据与随机分组

为了研究因果关系，我们一般都会通过随机实验来控制干预变量，随机实验中的随机分组可以排除一切混淆变量，从而保证实验的有效性。但是有些时候，我们无法进行实验，可供研究的只有已经生成的数据。由于没有经过随机分组，没有排除混淆变量，这时 Treatment variable 是不可用的。

随机分组的数据

观察性数据，分组不随机

1 倾向性评分匹配 Propensity score matching

倾向得分匹配是一种非实验性的因果推理技术。它试图在混淆变量上平衡干预组合对照组，使它们具有可比性，以便我们可以使用观察数据得出干预变量的因果效用的结论，它一般分为5个步骤

收集数据
计算 Propensity score
匹配样本
评估匹配结果
得出因果结论

1.1 收集数据

在进行 Propensity score matching，最重要的一步就是收集数据。在收集数据的时候，一定要把所有可能的Confounding variable 都囊括进来。如果有重要的Confounding variable 没有被纳入到数据集中，那么最后的匹配就很可能是无效的。

本例采用泰坦尼克号数据, 为了简洁，我们只使用少量的变量作为混淆变量。

混淆变量

这里，我们试图调查：购买三等舱，是否会影响到最后的生存率。而我们认为年纪，性别可能会同时影响到是否购买三等舱（年纪，性别也会影响经济水平）和是否获救（老人，妇女，儿童被优待）。

import numpy as np
import pandas as pd
pd.options.display.float_format = "{:.2f}".format
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='darkgrid', context='talk')
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.metrics import roc_auc_score, f1_score
from causalinference import CausalModel
df = sns.load_dataset('titanic')
df['is_pclass3'] = df['pclass']==3
df['is_female'] = df['sex']=='female'
df = df.filter(['survived', 'is_pclass3', 'is_female', 'age'])\
       .dropna().reset_index(drop=True)
df

image.png

TREATMENT = 'is_pclass3'
OUTCOME = 'survived'
df.groupby(TREATMENT)[OUTCOME].describe()

首先，我们先直接看一下三等舱与生存率的关系:

image.png

三等舱生存率是 24%，非三等舱生存率为 57%。

接下来，我们看一下，在不同的混淆变量下，干预组（三等舱）与对照组（其他舱）获救率的分布。

C_COLOUR = 'grey'
T_COLOUR = 'green'
C_LABEL = 'Control'
T_LABEL = 'Treatment'
sns.kdeplot(data=df[~df[TREATMENT]], x='age', shade=True, 
            color=C_COLOUR, label=C_LABEL)
sns.kdeplot(data=df[df[TREATMENT]], x='age', shade=True, 
            color=T_COLOUR, label=T_LABEL)
plt.legend();

可以看到，买三等舱票的人中，年轻人居多。

image.png

接下来，我们看看性别：

F_COLOUR = 'magenta'
M_COLOUR = 'blue'
F_LABEL = 'Female'
M_LABEL = 'Male'
gender = 100 * pd.crosstab(df[TREATMENT].replace({True: T_LABEL, 
                                                  False: C_LABEL}), 
                           df['is_female'].replace({True: 'Female',
                                                    False: 'Male'}), 
                           normalize='index')
gender['All'] = 100
plt.figure(figsize=(5, 4))
sns.barplot(data=gender, x=gender.index.astype(str),  y="All", 
            color=M_COLOUR, label=M_LABEL)
sns.barplot(data=gender, x=gender.index.astype(str),  y='Female', 
            color=F_COLOUR, label=F_LABEL)
plt.legend(loc='center', bbox_to_anchor=(1.3, 0.8))
plt.xlabel('')
plt.ylabel('Percentage');

image.png

很明显，两组中性别组成成分不一致。

从上面的分析可以看出，干预组与控制组的年龄，性别组成成分差异比较大，因此我们不能不考虑这些因素，直接得出三等舱影响生存率的结论。

1.2 计算 Propensity Scores

Propensity Scores 其实就是用混淆变量(age, is_female)来回归是否进入干预组（is_pclass3），因为 “是否进入干预组” 是一个二分类变量（0或者1表示），我们可以用 Logistic Regression。

# Build a descriptive model
t = df[TREATMENT]
X = pd.get_dummies(df.drop(columns=[OUTCOME, TREATMENT]))
pipe = Pipeline([
    ('scaler', StandardScaler()),
    ('logistic_classifier', LogisticRegression())
])
pipe.fit(X, t)
# Predict
threshold = 0.5
df['proba'] = pipe.predict_proba(X)[:,1]
df['logit'] = df['proba'].apply(lambda p: np.log(p/(1-p)))
df['pred'] = np.where(df['proba']>=threshold, 1, 0)
df.head()

image.png

proba 就是 Propensity Scores，他代表该条记录被分配到干预组中的概率。

接着我们看一下得到的模型的性能：

print(f"Accuracy: {np.mean(df[TREATMENT]==df['pred']):.4f},\
 ROC AUC: {roc_auc_score(df[TREATMENT], df['proba']):.4f},\
 F1-score: {f1_score(df[TREATMENT], df['pred']):.4f}")
# Visualise confusion matrix
pd.crosstab(df[TREATMENT], df['pred']).rename(columns={0: False, 
                                                       1:True})

模型性能

我们看一下 prob 和对应的 logit 在 control 和 treatment 组上的分布：

fig, ax = plt.subplots(1,2, figsize=(10,4))
# Visualise propensity
sns.kdeplot(data=df[~df[TREATMENT]], x='proba', shade=True, 
            color=C_COLOUR, label=C_LABEL, ax=ax[0])
sns.kdeplot(data=df[df[TREATMENT]], x='proba', shade=True, 
            color=T_COLOUR, label=T_LABEL, ax=ax[0])
ax[0].set_title('Propensity')
ax[0].legend(loc='center', bbox_to_anchor=(1.1, -0.3))
# Visualise logit propensity
sns.kdeplot(data=df[~df[TREATMENT]], x='logit', shade=True, 
            color=C_COLOUR, label=C_LABEL, ax=ax[1])
sns.kdeplot(data=df[df[TREATMENT]], x='logit', shade=True, 
            color=T_COLOUR, label=T_LABEL, ax=ax[1])
ax[1].set_title('Logit Propensity')
ax[1].set_ylabel("");

image.png

我们可以看到，干预组和控制组有重叠部分，这对接下来的匹配来说，是个好消息。

1.3 匹配

我们依据 Propensity Score ，把干预组和对照组中的得分最接近样本进行一对一匹配，注意，在匹配中会重复使用样本，因此有些样本会在结果中出现多次。

# Sort by 'logit' so it's quicker to find match
df.sort_values('logit', inplace=True)
n = len(df)-1
for i, (ind, row) in enumerate(df.iterrows()): 
    # Match the most similar untreated record to each treated record
    if row[TREATMENT]:
        # Find the closest untreated match among records sorted 
        # higher. 'equal_or_above would' be more accurate but 
        # used 'above' for brevity        
        if i<n:
            above = df.iloc[i:]
            control_above = above[~above[TREATMENT]]
            match_above = control_above.iloc[0]
            distance_above = match_above['logit'] - row['logit']
            df.loc[ind, 'match'] = match_above.name
            df.loc[ind, 'distance'] = distance_above
        
        # Find the closest untreated match among records sorted 
        # lower. 'equal_or_below' would be more accurate but 
        # used 'below' for brevity  
        if i>0:
            below = df.iloc[:i-1]
            control_below = below[~below[TREATMENT]]
            match_below = control_below.iloc[-1]
            distance_below = match_below['logit'] - row['logit']
            if i==n:
                df.loc[ind, 'match'] = match_below.name
                df.loc[ind, 'distance'] = distance_below
            
            # Only overwrite if match_below is closer than match_above
            elif distance_below<distance_above:
                df.loc[ind, 'match'] = match_below.name
                df.loc[ind, 'distance'] = distance_below
df[df[TREATMENT]]

image.png

我们重新创建一个 data frame

indices = df[df['match'].notna()].index.\
          append(pd.Index(df.loc[df['match'].notna(), 'match']))
matched_df = df.loc[indices].reset_index(drop=True)
matched_df

image.png

1.4 评估匹配结果

COLUMNS = ['age', 'is_female', OUTCOME]
matches = pd.merge(df.loc[df[TREATMENT], COLUMNS+['match']], 
                   df[COLUMNS], left_on='match', 
                   right_index=True, 
                   how='left', suffixes=('_t', '_c'))
matches

image.png

for var in ['logit', 'age']:
    print(f"{var} | Before matching")
    display(df.groupby(TREATMENT)[var].describe())
    print(f"{var} | After matching")
    display(matched_df.groupby(TREATMENT)[var].describe())

image.png

可以看到，匹配后的干预组和对照组非常接近了。我们在看一下在年龄上的分布

for var in ['logit', 'age']:
    fig, ax = plt.subplots(1,2,figsize=(10,4))
    # Visualise original distribution
    sns.kdeplot(data=df[~df[TREATMENT]], x=var, shade=True, 
                color=C_COLOUR, label=C_LABEL, ax=ax[0])
    sns.kdeplot(data=df[df[TREATMENT]], x=var, shade=True, 
                color=T_COLOUR, label=T_LABEL, ax=ax[0])
    ax[0].set_title('Before matching')
    
    # Visualise new distribution
    sns.kdeplot(data=matched_df[~matched_df[TREATMENT]], x=var, 
                shade=True, color=C_COLOUR, label=C_LABEL, ax=ax[1])
    sns.kdeplot(data=matched_df[matched_df[TREATMENT]], x=var, 
                shade=True, color=T_COLOUR, label=T_LABEL, ax=ax[1])
    ax[1].set_title('After matching')
    ax[1].set_ylabel("")
    plt.tight_layout()
ax[0].legend(loc='center', bbox_to_anchor=(1.1, -0.3));

image.png

可以看大， matching 以后的数据集平衡性好了很多。我们在看一下性别：

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
# Visualise original distribution
sns.barplot(data=gender, x=gender.index.astype(str), y="All", 
            color=M_COLOUR, label=M_LABEL, ax=ax[0])
sns.barplot(data=gender, x=gender.index.astype(str), y='Female', 
            color=F_COLOUR, label=F_LABEL, ax=ax[0])
ax[0].legend(loc='center', bbox_to_anchor=(1.1, -0.3))
ax[0].set_xlabel('')
ax[0].set_ylabel('Percentage')
ax[0].set_title('Before matching')
# Visualise new distribution
gender_after = 100 * pd.crosstab(
    matched_df[TREATMENT].replace({True: T_LABEL, False: C_LABEL}), 
    matched_df['is_female'].replace({True: 'Female', False: 'Male'}), 
    normalize='index'
)
gender_after['All'] = 100
sns.barplot(data=gender_after, x=gender_after.index.astype(str), 
            y="All", color=M_COLOUR, label=M_LABEL, ax=ax[1])
sns.barplot(data=gender_after, x=gender_after.index.astype(str), 
            y='Female', color=F_COLOUR, label=F_LABEL, ax=ax[1])
ax[1].set_xlabel('')
ax[1].set_title('After matching')
ax[1].set_ylabel('');

image.png

1.5 在匹配后的数据上进行评估

summary = matched_df.groupby(TREATMENT)[OUTCOME]\
                    .aggregate(['mean', 'std', 'count'])
summary

image.png

接着，我们评估一下 ATT 和 ATE。关于这些指标的意义，请参考[因果推断常用评估指标]。(https://www.jianshu.com/p/5ebb3b2a8f55)

y = df[OUTCOME].values
t = df[TREATMENT].values
X = df[['is_female', 'age']]
X = pd.DataFrame(StandardScaler().fit_transform(X), 
                 columns=X.columns).values
model = CausalModel(y, t, X)
model.est_via_matching()
print(model.estimates)

image.png

通过ATT，我们可以得出三等舱确实会影响生存率的结论（降低）。虽然和我们一开始得出的结论相同，但是由于排除了混淆因素，现在的结论更加有效。

倾向性评分匹配（Propensity Score Matchi

术语

观察数据与随机分组

1 倾向性评分匹配 Propensity score matching

1.1 收集数据

1.2 计算 Propensity Scores

1.3 匹配

1.4 评估匹配结果

1.5 在匹配后的数据上进行评估

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

倾向性评分匹配 （Propensity Score Matchi

术语

观察数据与随机分组

1 倾向性评分匹配 Propensity score matching

1.1 收集数据

1.2 计算 Propensity Scores

1.3 匹配

1.4 评估匹配结果

1.5 在匹配后的数据上进行评估

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

倾向性评分匹配（Propensity Score Matchi