傾向スコアのAUCの良さは推定量の性質の良さにつながるのか？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#%%
import numpy as np
import pandas as pd
import plotly.express as px
import statsmodels.api as sm 
from sklearn import metrics

#-------------
# DGP
#-------------
def dgp():
    N = 10000
    x1 = np.random.uniform(0, 1, N)
    x2 = np.random.uniform(0, 1, N)
    x3 = np.random.uniform(0, 1, N)

    logit_model = np.exp(-0.5 + 2 * x1 - 2 * x2 - 2 * x1*x2 + 1 * x3)
    assignment_prob = logit_model / (1 + logit_model)
    assigment = np.random.binomial(1, assignment_prob)

    potential_untreated = np.exp(2 * x1) + np.random.normal(0, 0.5, N)
    uplift = 1
    potential_treated = potential_untreated + uplift


    df = pd.DataFrame({
        'const': 1,
        'x1': x1,
        'x2': x2,
        'x3': x3,
        'x1x2': x1*x2,
        'prob': assignment_prob,
        'D': assigment,
        'uplift':uplift, 
        'y0': potential_untreated,
        'y1': potential_treated
    })
    df['observed'] = (1 - df['D']) * df['y0'] + (df['D']) * df['y1']

    return df

def logit_compute(df):
    y = df.loc[:, 'D'].values
    x_simple = df.loc[:, ['const', 'x1']]
    x_true = df.loc[:, ['const', 'x1', 'x2', 'x1x2', 'x3']]

    logit_reg_simple = sm.Logit(y, x_simple).fit()
    logit_reg_true = sm.Logit(y, x_true).fit()

    df['logit_cia'] = logit_reg_simple.predict(x_simple)
    df['logit_full'] = logit_reg_true.predict(x_true)
    return df

def get_auc(df, col):
    fpr, tpr = metrics.roc_curve(df['D'], df[col], pos_label=1)[:2]
    return metrics.auc(fpr, tpr)

def compute_normalized_diff(df, prob='logit_cia'):
    features = ['x1', 'x2', 'x3']
    snd = []

    wtavg = lambda x: np.average(x.loc[:,features], weights = x[prob].values, axis = 0)
    tmp_df = df.loc[:, features + ['D', prob]]
    weighted_mean = tmp_df.groupby(['D']).apply(wtavg)
    weighted_var = tmp_df.groupby(['D']).var()
    
    for i, j in enumerate(features):
        mean_dif = abs(weighted_mean[0][i] - weighted_mean[1][i])
        denominator = np.sqrt(np.sum(weighted_var[j]))
        smd = mean_dif/denominator
        snd.append(smd)
    
    return pd.DataFrame({'features': features, 'snd': snd})

df = dgp()
df = logit_compute(df)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#--------------------
# True ATE, ATT, ATU
#--------------------
def get_true_effect(df):
    stats_df = df.groupby('D')[['y0', 'y1']].mean().reset_index(drop=True)
    stats_df = pd.concat([stats_df, 
                          pd.DataFrame(df.loc[:, ['y0', 'y1']].mean(axis=0)).T], 
                          axis=0).set_index(pd.Index(['ATU', 'ATT', 'ATE']))
    stats_df['effect'] = stats_df['y1'] - stats_df['y0']
    return stats_df

def estimate_effect(df, col='logit_cia'):
    y0_treated_weight = ((1 - df['D']) * df[col]) / (1 - df[col])
    y1_untreated_weight = (df['D']) * (1 - df[col]) / df[col]
    y1_ate_weight = df['D'] / df[col]
    y0_ate_weight = (1 - df['D'])/(1 - df[col])

    y0_untreated, y1_treated = df.groupby('D')['observed'].mean().values
    y0_treated = np.sum(y0_treated_weight * df['observed']) / np.sum(y0_treated_weight)
    y1_untreated = np.sum(y1_untreated_weight * df['observed']) / np.sum(y1_untreated_weight)
    
    y1_ate = np.sum(y1_ate_weight * df['observed'])/ np.sum(y1_ate_weight)
    y0_ate = np.sum(y0_ate_weight * df['observed'])/ np.sum(y0_ate_weight)

    res_df = pd.DataFrame({'y0':[y0_untreated, y0_treated, y0_ate],
                           'y1':[y1_untreated, y1_treated, y1_ate]})
    res_df = res_df.set_index(pd.Index(['ATU', 'ATT', 'ATE']))
    res_df['effect'] = res_df['y1'] - res_df['y0']

    return res_df

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#-----------------
# simulation
#-----------------
np.random.seed(42)

res_true = []
res_estimates_cia = []
res_estimates_full = []
cia_auc = []
cia_full = []    

for i in range(1000):
    tmp = dgp()
    tmp = logit_compute(tmp)
    res_true.append(get_true_effect(tmp).loc['ATE', 'effect'])
    res_estimates_cia.append(estimate_effect(tmp).loc['ATE', 'effect'])
    res_estimates_full.append(estimate_effect(tmp, col='logit_full').loc['ATE', 'effect'])
    cia_auc.append(get_auc(tmp, 'logit_cia'))
    cia_full.append(get_auc(tmp, 'logit_full'))

compare_df = pd.DataFrame({'true':res_true, 
                           'cia_model_ATE':res_estimates_cia, 
                           'full_model_ATE':res_estimates_full,
                           'cia_model_AUC': cia_auc,
                           'full_model_AUC': cia_full})
compare_df.describe()