Cet article est basé sur "Introduction à la vérification de l'efficacité - Raisonnement causal pour une comparaison correcte / Bases de l'économie quantitative" Ceci est un résumé de ce que j'ai appris et ressenti en lisant «Chapitre 3 Analyse à l'aide des scores de propension».

Basiques

Qu'est-ce qu'un score de propension?

** Observé ** Probabilité que la covariable $ X $ soit attribuée au groupe d'intervention dans des conditions données $ e (X) = P (Z = 1 | X) $ [^ px]
Toutes les covariables ne sont pas observées

[^ px]: Dans ce livre, le score de propension s'écrit $ P (X) $, mais pour plus de commodité, il est défini comme $ e (X) $.

Analyse régressive par rapport au score de propension [^ ra-vs-ps]

[^ ra-vs-ps]: Une petite comparaison étrange, mais à partir de ce livre, on peut la lire comme ceci

Hypothèses d'analyse de régression:
- CIA(Conditional Independence Assumption): \\{Y^{(1)}, Y^{(0)}\\}\perp Z|X
Linéarité des résultats et des covariables (omis dans ce livre?) --Hypothèse du score de propension --Version CIA: $ \ {Y ^ {(1)}, Y ^ {(0)} \} \ perp Z | e (X) $
Je pense que c'est la même chose d'expliquer l'allocation avec la covariable $ X $ et le score de propension $ e (X) $ pour équilibrer la covariable [^ balancing-score]
** SUTVA ** (Hypothèse de valeur de traitement unitaire stable) [^ sutva-krsk-phs] [^ sutva-pira-nino] [^ sutva-in-ra](Omis dans ce livre?)

CIA c.s. Affectation fortement négligeable

--CIA est utilisé dans ce livre, mais à proprement parler, ** affectation fortement négligeable ** est correcte [^ rosenbaum_1983]

Une allocation fortement négligeable ajoute les deux hypothèses suivantes à la CIA:
Les variables $ X $ ne sont
- 0 «En bref, je pense que c'est ** assez important ** parce que c'est NG s'il n'y a aucune possibilité qu'un échantillon soit attribué à la fois au groupe d'intervention et au groupe de non-intervention. --Par exemple, si tous les hommes sont affectés au groupe d'intervention, $ P (Z = 1 | hommes) = 1 $, donc NG

Estimation du score de propension

La régression logistique est souvent utilisée pour estimer les scores de propension, mais l'apprentissage automatique est également OK --Le modèle comprend non seulement des covariables, mais également d'autres variables qui affectent les résultats [^ tsugawa-ps2]
L'étalonnage des probabilités est requis lors de l'utilisation de l'apprentissage automatique [^ rmizutaa]
Forêt aléatoire et SVM? --D'un autre côté, si la fonction de perte est logloss, l'étalonnage n'est pas nécessaire car il converge vers la probabilité [^ étalonnage-inutile]

Estimation de l'effet de l'intervention à l'aide du score de propension

ATT v.s. ATE

ATT(Average Treatment effect on Treated):
- E[Y^{(1)}-Y^{(0)}|Z=1]
Valeur attendue de l'effet de l'intervention dans l'échantillon ayant reçu l'intervention --Dans ce livre, il est calculé en faisant correspondre
ATE(Average Treatment Effect):
- E[Y^{(1)}-Y^{(0)}]
Valeur attendue de l'effet d'intervention pour l'ensemble de l'échantillon --Dans ce livre, il est calculé par IPW

Correspondance v.s. IPW [^ matching-vs-ipw]

[^ matching-vs-ipw]: Il existe d'autres méthodes pour estimer l'effet de l'intervention en utilisant des scores de propension autres que l'appariement et l'IPW, mais ce livre traite de ces deux.

-** Correspondant à **: --Une paire d'échantillons avec des valeurs de score de propension similaires est prélevée dans le groupe d'intervention et le groupe sans intervention, et une mesure manquante (contrefaçon) est complétée par l'autre.

IPW(Inverse Probability Weighting):
Peser l'échantillon avec l'inverse du score de propension à se rapprocher de la distribution de la population --IPW ne considère pas le support commun (décrit plus loin) et ne doit pas être utilisé? [^ tsugawa-ps2] ――Je pense que vous pouvez l'utiliser si vous considérez un support commun

ATT en correspondance

Appariement d'échantillons avec des scores de propension similaires à ceux des échantillons du groupe d'intervention du groupe sans intervention, et en laissant la moyenne des différences de paires être ATT

ATE en correspondance

De même pour les échantillons du groupe sans intervention, soit ATE la moyenne de la différence entre les paires

ATT in IPW

--Calculé par la formule suivante

\frac{\sum_{i=1}^{N}Z_{i}Y_{i}}{\sum_{i=1}^{N}Z_{i}}
-\frac{\sum_{i=1}^{N}\frac{(1-Z_{i})\hat{e}(X_{i})}{1-\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{(1-Z_{i})\hat{e}(X_{i})}{1-\hat{e}(X_{i})}}

ATE in IPW

--Calculé par la formule suivante

\frac{\sum_{i=1}^{N}\frac{Z_{i}}{\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{Z_{i}}{\hat{e}(X_{i})}}
-\frac{\sum_{i=1}^{N}\frac{1-Z_{i}}{1-\hat{e}(X_{i})}Y_{i}}{\sum_{i=1}^{N}\frac{1-Z_{i}}{1-\hat{e}(X_{i})}}

À propos du support commun

Le soutien commun est la fourchette dans laquelle la distribution des scores de propension du groupe d'intervention et du groupe de non-intervention se chevauchent. --Si vous utilisez un échantillon qui n'est pas pris en charge en commun dans le calcul de l'IPW, le nombre inverse (poids) du score de propension devient trop grand pour bien estimer [^ tsugawa-ps1] -Est-ce que "hors support commun" et "poids trop grand" équivalent?
Est-il courant d'exclure des échantillons qui ne sont pas dans le support commun? --Ce livre présente la propriété partition trim mig / clipping ―― N'est-il pas possible de créer un nouveau biais en excluant le support non commun? ――Il est donc important de concevoir le plan expérimental de manière à pouvoir observer correctement deux groupes comparables afin de ne pas avoir à l'exclure ** ――Le biais généré à ce sujet doit être ajusté par le score de tendance, etc. ―― Puisqu'il n'y a pas d'échantillon qui puisse compléter le contre-fait en dehors du soutien commun, ne satisfait-il pas en premier lieu «l'allocation fortement négligeable»?
- P(Z=1|X)=0OuP(Z=1|X)=1Pas parce qu'il utilise la régression logistique

Évaluation du score de propension

—— Il est important que le score de propension équilibre la distribution des covariables entre les groupes d'intervention et de non-intervention -Autrement dit, faire des hypothèses\\{Y^{(1)}, Y^{(0)}\\}\perp Z|e(X)De\\{Y^{(1)}, Y^{(0)}\\}\perp Z|XJe pense que cela signifie revenir à

Dans ce livre, la différence moyenne standardisée (ASAM) est utilisée pour l'évaluation. --OK si ASAM est de 0,1 ou moins (pourquoi?) ――Si l'ancienne AUC est de 0,8 ou plus, la théorie OK est-elle absurde? [^ ps-auc] --C'est OK si AUC est égal ou supérieur à 0,8 [^ iwanami-ds3]
Qu'est-ce que cela signifie que l'ASC est petite (proche de 0,5):
L'allocation ne dépend pas des covariables ($ Z \ perp X $) → Pas besoin d'ajuster avec le score de propension → Est-ce que cela n'a pas de sens si l'AUC n'est pas grande dans une certaine mesure?
Que signifie avoir une grande AUC (proche de 1): ――Le soutien commun est faible → L'effet de l'intervention ne peut être estimé → N'est-il pas bon que l'AUC soit grande?

Édition pratique

MineThatData E-Mail Analytics And Data Mining Challenge dataset^[1]

Sera ajouté à une date ultérieure

LaLonde(1986)

Reportez-vous à ce manuel pour l'explication des données
Voir ici pour l'explication de l'article.

Lecture des données

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

df_cps1 = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls.dta')
df_cps3 = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls3.dta')
df_nsw = pd.read_stata('https://users.nber.org/~rdehejia/data/nsw_dw.dta')

Création de l'ensemble de données

--NSW: ensemble de données RCT --CPS1 + NSW: ensemble de données avec groupe de non-intervention NSW remplacé par CPS1 --CPS3 + NSW: ensemble de données avec groupe de non-intervention NSW remplacé par CPS3

treat_label = {0: 'untreat', 1: 'treat'}
X_columns = ['age', 'education', 'black', 'hispanic', 'married', 'nodegree', 're74', 're75']
y_column = 're78'
z_column = 'treat'

CPS1+NSW

df_cps1_nsw = pd.concat([df_nsw[df_nsw[z_column] == 1], df_cps1], ignore_index=True)

CPS3+NSW

df_cps3_nsw = pd.concat([df_nsw[df_nsw[z_column] == 1], df_cps3], ignore_index=True)

Confirmer la répartition du groupe d'intervention et du groupe de non-intervention

NSW

Les données sont bien équilibrées car elles ne changent que RCT.

df_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_nsw[c].groupby(df_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS1+NSW

――Ce sont des données assez déséquilibrées

df_cps1_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS1+NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

Pas aussi bon que NSW, mais plus équilibré que CPS1 + NSW

df_cps3_nsw[z_column].map(treat_label).value_counts().plot.bar(title='# of treatment in CPS3+NSW')

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

analyse de régression

import statsmodels.api as sm

NSW

results = sm.OLS(df_nsw[y_column], sm.add_constant(df_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

CPS1+NSW

results = sm.OLS(df_cps1_nsw[y_column], sm.add_constant(df_cps1_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

CPS3+NSW

results = sm.OLS(df_cps3_nsw[y_column], sm.add_constant(df_cps3_nsw[[z_column] + X_columns])).fit()
results.summary().tables[1]

résultat

base de données	Effet d'intervention[$]
NSW	1676.3
CPS1+NSW	699.1
CPS3+NSW	1548.2

Début du jeu

L'effet d'intervention de NSW est d'environ 1600 $ --Défiez-vous à quel point vous pouvez vous rapprocher de 1600 $ en utilisant le score de propension d'un état biaisé

Estimation du score de propension

D'où viennent «je (re74 ^ 2)» et «je (re75 ^ 2)» dans ce livre?

CPS1+NSW

--Estimation du score de propension par régression logistique

ps_model = sm.Logit(df_cps1_nsw[z_column], sm.add_constant(df_cps1_nsw[X_columns])).fit()
ps_model.summary().tables[1]

df_cps1_nsw['ps'] = 1 / (1 + np.exp(-ps_model.fittedvalues))

Vérifier la distribution des scores de propension

(df_cps1_nsw.groupby(df_cps1_nsw[z_column].map(treat_label))['ps']
 .plot.hist(density=True, alpha=0.5, title=f'propensity score densities', legend=True))

--Calculer AUC

from sklearn.metrics import roc_auc_score
roc_auc_score(df_cps1_nsw[z_column], df_cps1_nsw['ps'])

0.9708918648513447

CPS3+NSW

--Estimation du score de propension par régression logistique

ps_model = sm.Logit(df_cps3_nsw[z_column], sm.add_constant(df_cps3_nsw[X_columns])).fit()
ps_model.summary().tables[1]

df_cps3_nsw['ps'] = 1 / (1 + np.exp(-ps_model.fittedvalues))

Vérifier la distribution des scores de propension

(df_cps3_nsw.groupby(df_cps3_nsw[z_column].map(treat_label))['ps']
 .plot.hist(density=True, alpha=0.5, title=f'propensity score densities', legend=True))

--Calculer AUC

from sklearn.metrics import roc_auc_score
roc_auc_score(df_cps3_nsw[z_column], df_cps3_nsw['ps'])

0.8742266742266742

Estimation ATT

CPS1+NSW

correspondant à

--Définissez une fonction (générateur récursif) qui crée une paire

Ici a été utilisé comme référence.

from sklearn.neighbors import NearestNeighbors
def pairs_generator(s0: pd.Series, s1: pd.Series, threshold: np.float) -> pd.DataFrame:
    """
Un générateur récursif qui crée des paires à l'aide de la méthode K-voisinage.
    
    Parameters
    ----------
    s0, s1 : pd.Series
    threshold : np.float
Le seuil de la distance pour faire une paire. Si le point voisin le plus proche dépasse également le seuil, il ne sera pas apparié.
    
    Returns
    -------
    pairs : pd.DataFrame
Index d'une paire composée de l'échantillon s0 et de l'échantillon s1 le plus proche.
La colonne l est l'indice de s0 et la colonne r est l'indice de s1.
    
    Notes
    -----
Le même échantillon n'est jamais utilisé pour plusieurs paires.
    
    Examples
    --------
    s0 = df[df[z_column] == 0]['ps']
    s1 = df[df[z_column] == 1]['ps']
    threshold = df['ps'].std() * 0.1
    pairs = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
    """
    neigh_dist, neigh_ind = NearestNeighbors(1).fit(s0.values.reshape(-1, 1)).kneighbors(s1.values.reshape(-1, 1))
    pairs = pd.DataFrame(
        {'d': neigh_dist[:, 0], 'l': s0.iloc[neigh_ind[:, 0]].index},
        index=s1.index,
    ).query(f'd < {threshold}').groupby('l')['d'].idxmin().rename('r').reset_index()
    if len(pairs) > 0:
        yield pairs
        yield from pairs_generator(s0.drop(pairs['l']), s1.drop(pairs['r']), threshold)

--Faire une paire

s1 = df_cps1_nsw[df_cps1_nsw[z_column] == 1]['ps']
s0 = df_cps1_nsw[df_cps1_nsw[z_column] == 0]['ps']
threshold = df_cps1_nsw['ps'].std() * 0.1
pairs1_cps1_nsw = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
pairs0_cps1_nsw = pd.concat(pairs_generator(s0, s1, threshold), ignore_index=True) #Utilisé dans le calcul ATE

--Confirmer une paire de groupes d'intervention --Orange est un échantillon jeté ―― À l'origine, il est souhaitable de faire une paire afin que l'intersection comme indiqué sur la figure ne se produise pas (je pense)

ax = df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs on treated')
df_cps1_nsw.loc[~df_cps1_nsw.index.isin(pairs1_cps1_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps1_nsw['ps'].iloc[pairs1_cps1_nsw['l']].values, df_cps1_nsw['ps'].iloc[pairs1_cps1_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

--Calculer ATT

L'écart type est assez grand

(pairs1_cps1_nsw['l'].map(df_cps1_nsw[y_column]) - pairs1_cps1_nsw['r'].map(df_cps1_nsw[y_column])).agg(['mean', 'std'])

mean    1929.083008
std     9618.346680
dtype: float64

IPW

--Calculer le poids

df_cps1_nsw['w'] = df_cps1_nsw[z_column] / df_cps1_nsw['ps'] + (1 - df_cps1_nsw[z_column]) / (1 - df_cps1_nsw['ps'])

--Estimation ATT avec régression linéaire pondérée

att_model = sm.WLS(df_cps1_nsw[y_column], df_cps1_nsw[[z_column]].assign(untreat=1-df_cps1_nsw[z_column]), weights=df_cps1_nsw['w'] * df_cps1_nsw['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1180.4078220960955

IPW (prise en charge commune)

Extraire des échantillons dans la plage où les distributions des scores de propension se chevauchent

df_cps1_nsw_cs = df_cps1_nsw[df_cps1_nsw['ps'].between(*df_cps1_nsw.groupby(z_column)['ps'].agg(['min', 'max']).agg({'min': 'max', 'max': 'min'}))].copy()

--Estimation ATT avec régression linéaire pondérée

att_model = sm.WLS(df_cps1_nsw_cs[y_column], df_cps1_nsw_cs[[z_column]].assign(untreat=1-df_cps1_nsw_cs[z_column]), weights=df_cps1_nsw_cs['w'] * df_cps1_nsw_cs['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1243.8640545001808

--Confirmer la différence moyenne standardisée entre les deux groupes de covariables ――Après ajustement, il se trouve généralement sous la ligne pointillée rouge.

Est-ce que cela correspond au graphique de ce livre? (Surtout noir avant réglage)

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps1_nsw.groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps1_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps1_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps1_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps1_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

--Confirmer la distribution des covariables dans l'appariement --La moyenne ne dit pas la vérité, alors regardez la distribution ――L'équilibrage est fait correctement par rapport à la distribution d'origine

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw.loc[pairs1_cps1_nsw.unstack()][c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW ――L'équilibrage est fait correctement par rapport à la distribution d'origine

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw.sample(1000000, replace=True, weights=df_cps1_nsw['w'] * df_cps1_nsw['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW (support commun) ――L'équilibrage est fait correctement par rapport à la distribution d'origine

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw_cs.sample(1000000, replace=True, weights=df_cps1_nsw_cs['w'] * df_cps1_nsw_cs['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

correspondant à

--Faire une paire

s1 = df_cps3_nsw[df_cps3_nsw[z_column] == 1]['ps']
s0 = df_cps3_nsw[df_cps3_nsw[z_column] == 0]['ps']
threshold = df_cps3_nsw['ps'].std() * 0.1
pairs1_cps3_nsw = pd.concat(pairs_generator(s1, s0, threshold), ignore_index=True)
pairs0_cps3_nsw = pd.concat(pairs_generator(s0, s1, threshold), ignore_index=True) #Utilisé dans le calcul ATE

--Confirmer une paire de groupes d'intervention

ax = df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs on treated')
df_cps3_nsw.loc[~df_cps3_nsw.index.isin(pairs1_cps3_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps3_nsw['ps'].iloc[pairs1_cps3_nsw['l']].values, df_cps3_nsw['ps'].iloc[pairs1_cps3_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

--Calculer ATT

(pairs1_cps3_nsw['l'].map(df_cps3_nsw[y_column]) - pairs1_cps3_nsw['r'].map(df_cps3_nsw[y_column])).agg(['mean', 'std'])

mean    1463.115845
std     8779.598633
dtype: float64

IPW

--Calculer le poids

df_cps3_nsw['w'] = df_cps3_nsw[z_column] / df_cps3_nsw['ps'] + (1 - df_cps3_nsw[z_column]) / (1 - df_cps3_nsw['ps'])

--Estimation ATT avec régression linéaire pondérée

att_model = sm.WLS(df_cps3_nsw[y_column], df_cps3_nsw[[z_column]].assign(untreat=1-df_cps3_nsw[z_column]), weights=df_cps3_nsw['w'] * df_cps3_nsw['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

1214.0711733143507

IPW (prise en charge commune)

Extraire des échantillons dans la plage où les distributions des scores de propension se chevauchent

df_cps3_nsw_cs = df_cps3_nsw[df_cps3_nsw['ps'].between(*df_cps3_nsw.groupby(z_column)['ps'].agg(['min', 'max']).agg({'min': 'max', 'max': 'min'}))].copy()

--Estimation ATT avec régression linéaire pondérée

att_model = sm.WLS(df_cps3_nsw_cs[y_column], df_cps3_nsw_cs[[z_column]].assign(untreat=1-df_cps3_nsw_cs[z_column]), weights=df_cps3_nsw_cs['w'] * df_cps3_nsw_cs['ps']).fit()
att_model.summary().tables[1]

att_model.params['treat'] - att_model.params['untreat']

988.8657151185471

--Confirmer la différence moyenne standardisée entre les deux groupes de covariables

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps3_nsw.groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps3_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps3_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps3_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w'] * x['ps']))) / df_cps3_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

--Confirmer la distribution des covariables dans l'appariement

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw.loc[pairs1_cps3_nsw.unstack()][c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw.sample(1000000, replace=True, weights=df_cps3_nsw['w'] * df_cps3_nsw['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW (support commun)

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw_cs.sample(1000000, replace=True, weights=df_cps3_nsw_cs['w'] * df_cps3_nsw_cs['ps'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

résultat

――Les covariables sont équilibrées, mais la valeur estimée de l'effet d'intervention semble être assez différente de 1600 $.

base de données	correspondant à	IPW	IPW(Support commun)
CPS1+NSW	1929.1	1180.4	1243.9
CPS3+NSW	1463.1	1214.1	988.9

Estimation ATE

CPS1+NSW

correspondant à

--Vérifiez la paire

pairs_cps1_nsw = pd.concat([pairs1_cps1_nsw, pairs0_cps1_nsw.rename({'l': 'r', 'r': 'l'}, axis=1)], ignore_index=True)

ax = df_cps1_nsw.loc[pairs_cps1_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs')
df_cps1_nsw.loc[~df_cps1_nsw.index.isin(pairs_cps1_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps1_nsw['ps'].iloc[pairs_cps1_nsw['l']].values, df_cps1_nsw['ps'].iloc[pairs_cps1_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

--Calculer ATE

(pairs_cps1_nsw['l'].map(df_cps1_nsw[y_column]) - pairs_cps1_nsw['r'].map(df_cps1_nsw[y_column])).agg(['mean', 'std'])

mean    1836.516968
std     9634.636719
dtype: float64

IPW

--Estimation ATE avec régression linéaire pondérée

ate_model = sm.WLS(df_cps1_nsw[y_column], df_cps1_nsw[[z_column]].assign(untreat=1-df_cps1_nsw[z_column]), weights=df_cps1_nsw['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

-6456.300804768925

IPW (prise en charge commune)

--Estimation ATE avec régression linéaire pondérée

ate_model = sm.WLS(df_cps1_nsw_cs[y_column], df_cps1_nsw_cs[[z_column]].assign(untreat=1-df_cps1_nsw_cs[z_column]), weights=df_cps1_nsw_cs['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

545.3149727568589

--Confirmer la différence moyenne standardisée entre les deux groupes de covariables

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps1_nsw.groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps1_nsw.loc[pairs_cps1_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps1_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps1_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps1_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps1_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps1_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

--Confirmer la distribution des covariables dans l'appariement

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps1_nsw.loc[pairs_cps1_nsw.unstack()][c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw.sample(1000000, replace=True, weights=df_cps1_nsw['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW (support commun)

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps1_nsw_cs.sample(1000000, replace=True, weights=df_cps1_nsw_cs['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps1_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

CPS3+NSW

correspondant à

--Vérifiez la paire

pairs_cps3_nsw = pd.concat([pairs1_cps3_nsw, pairs0_cps3_nsw.rename({'l': 'r', 'r': 'l'}, axis=1)], ignore_index=True)

ax = df_cps3_nsw.loc[pairs_cps3_nsw.unstack()].plot.scatter(x='ps', y=z_column, label='matched', alpha=0.3, figsize=(16, 4), title='pairs')
df_cps3_nsw.loc[~df_cps3_nsw.index.isin(pairs_cps3_nsw.unstack())].plot.scatter(x='ps', y=z_column, label='unmatched', alpha=0.3, c='tab:orange', ax=ax)
for x in zip(df_cps3_nsw['ps'].iloc[pairs_cps3_nsw['l']].values, df_cps3_nsw['ps'].iloc[pairs_cps3_nsw['r']].values):
    plt.plot(x, [1, 0], c='tab:blue', alpha=0.3)

--Calculer ATE

(pairs_cps3_nsw['l'].map(df_cps3_nsw[y_column]) - pairs_cps3_nsw['r'].map(df_cps3_nsw[y_column])).agg(['mean', 'std'])

mean    1486.608276
std     8629.639648
dtype: float64

IPW

--Estimation ATE avec régression linéaire pondérée

ate_model = sm.WLS(df_cps3_nsw[y_column], df_cps3_nsw[[z_column]].assign(untreat=1-df_cps3_nsw[z_column]), weights=df_cps3_nsw['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

224.6762634665365

IPW (prise en charge commune)

--Estimation ATE avec régression linéaire pondérée

ate_model = sm.WLS(df_cps3_nsw_cs[y_column], df_cps3_nsw_cs[[z_column]].assign(untreat=1-df_cps3_nsw_cs[z_column]), weights=df_cps3_nsw_cs['w']).fit()
ate_model.summary().tables[1]

ate_model.params['treat'] - ate_model.params['untreat']

801.0736196669177

--Confirmer la différence moyenne standardisée entre les deux groupes de covariables

pd.DataFrame([{
    'unadjusted': np.subtract(*df_cps3_nsw.groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(matching)': np.subtract(*df_cps3_nsw.loc[pairs_cps3_nsw.unstack()].groupby(z_column)[c].mean()) / df_cps3_nsw[c].std(),
    'adjusted(weighting)': np.subtract(*df_cps3_nsw.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps3_nsw[c].std(),
    'adjusted(weighting in common support)': np.subtract(*df_cps3_nsw_cs.groupby(z_column).apply(lambda x: np.average(x[c], weights=x['w']))) / df_cps3_nsw_cs[c].std(),
} for c in X_columns], index=X_columns).abs().plot(marker='o', alpha=0.5, grid=True, ylim=(0, 100), logy=True, title='covariates balance')
plt.axhline(y=0.1, linestyle='--', c='r')

--Confirmer la distribution des covariables dans l'appariement

_, ax = plt.subplots(3, 3, figsize=(16, 12))
for i, c in enumerate(X_columns + [y_column]):
    (df_cps3_nsw.loc[pairs_cps3_nsw.unstack()][c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw.sample(1000000, replace=True, weights=df_cps3_nsw['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

--Confirmer la distribution des covariables dans IPW (support commun)

_, ax = plt.subplots(3, 3, figsize=(16, 12))
df = df_cps3_nsw_cs.sample(1000000, replace=True, weights=df_cps3_nsw_cs['w'], random_state=2020)
for i, c in enumerate(X_columns + [y_column]):
    (df[c].groupby(df_cps3_nsw_cs[z_column].map(treat_label))
     .plot.hist(density=True, alpha=0.5, title=f'{c} densities', legend=True, ax=ax[i // 3, i % 3]))

résultat

base de données	correspondant à	IPW	IPW(Support commun)
CPS1+NSW	1836.5	-6456.3	545.3
CPS3+NSW	1486.6	224.7	801.0

Impressions

Je ne peux pas du tout me débarrasser du "biais" ...

https://blog.minethatdata.com/2008/03/minethatdata-e-mail-analytics-and-data.html ↩︎

"Introduction à la vérification des effets Chapitre 3 Analyse utilisant le score de propension" + α est essayé en Python

Basiques

Qu'est-ce qu'un score de propension?

Analyse régressive par rapport au score de propension [^ ra-vs-ps]

CIA c.s. Affectation fortement négligeable

Estimation du score de propension

Estimation de l'effet de l'intervention à l'aide du score de propension

Correspondance v.s. IPW [^ matching-vs-ipw]

ATT en correspondance

ATE en correspondance

À propos du support commun

Évaluation du score de propension

Édition pratique

Lecture des données

Création de l'ensemble de données

Confirmer la répartition du groupe d'intervention et du groupe de non-intervention

analyse de régression

résultat

Début du jeu

Estimation du score de propension

Estimation ATT

correspondant à

IPW (prise en charge commune)

correspondant à

IPW (prise en charge commune)

résultat

Estimation ATE

correspondant à

IPW (prise en charge commune)

correspondant à

IPW (prise en charge commune)

résultat

Impressions