Dieser Artikel basiert auf "Einführung in die Überprüfung der Wirksamkeit - kausales Denken für einen korrekten Vergleich / Grundlagen der quantitativen Ökonomie" Dies ist eine Zusammenfassung dessen, was ich durch das Lesen von "Kapitel 3 Analyse mit Propensity Scores" gelernt und gefühlt habe.

Grundlagen

Was ist ein Neigungswert?

** Beobachtet ** Wahrscheinlichkeit, dass die Kovariate $ X $ unter den gegebenen Bedingungen der Interventionsgruppe zugeordnet wird $ e (X) = P (Z = 1 | X) $ [^ px]
Nicht alle Kovariaten werden beobachtet

[^ px]: In diesem Buch wird der Neigungswert als $ P (X) $ geschrieben, der Einfachheit halber wird er jedoch als $ e (X) $ festgelegt.

Regressive Analyse gegen Propensity Score [^ ra-vs-ps]

[^ ra-vs-ps]: Ein etwas seltsamer Vergleich, aber aus diesem Buch kann man es so lesen

Annahmen zur Regressionsanalyse:
- CIA(Conditional Independence Assumption): \\{Y^{(1)}, Y^{(0)}\\}\perp Z|X
Ergebnis und kovariate Linearität (in diesem Buch weggelassen?)
Annahme der Neigungsbewertung
CIA-Version: $ \ {Y ^ {(1)}, Y ^ {(0)} \} \ perp Z | e (X) $
Ich denke, es ist dasselbe, die Zuordnung mit der Kovariate $ X $ und der Neigungsbewertung $ e (X) $ zu erklären, um die Kovariate auszugleichen [^ balancing-score]
** SUTVA ** (Annahme eines stabilen Einheitsbehandlungswerts) [^ sutva-krsk-phs] [^ sutva-pira-nino] [^ sutva-in-ra](in diesem Buch weggelassen?)

CIA gegen stark vernachlässigbare Zuordnung

--CIA wird in diesem Buch verwendet, aber genau genommen ist ** stark vernachlässigbare Zuordnung ** korrekt [^ rosenbaum_1983]

Eine stark vernachlässigbare Allokation fügt der CIA die folgenden zwei Annahmen hinzu: --Kovariaten $ X $ werden nur beobachtet
- 0 Kurz gesagt, ich denke, dass es ** ziemlich wichtig ** ist, weil es NG ist, wenn es keine Möglichkeit gibt, dass eine Stichprobe sowohl der Interventionsgruppe als auch der Nicht-Interventionsgruppe zugewiesen wird.
Wenn zum Beispiel alle Männer der Interventionsgruppe zugeordnet sind, ist $ P (Z = 1 | Männer) = 1 $, also NG

Schätzung der Neigungsbewertung

--Logistische Regression wird häufig verwendet, um Neigungswerte zu schätzen, aber maschinelles Lernen ist auch in Ordnung

Das Modell enthält nicht nur Kovariaten, sondern auch andere Variablen, die die Ergebnisse beeinflussen [^ tsugawa-ps2]
Wahrscheinlichkeitskalibrierung ist erforderlich, wenn maschinelles Lernen verwendet wird [^ rmizutaa]
Zufälliger Wald und SVM?
Wenn andererseits die Verlustfunktion logarithmisch ist, ist eine Kalibrierung nicht erforderlich, da sie zur Wahrscheinlichkeit konvergiert [^ Kalibrierung nicht erforderlich]

Schätzung des Interventionseffekts anhand des Neigungsscores

ATT v.s. ATE

ATT(Average Treatment effect on Treated):
- E[Y^{(1)}-Y^{(0)}|Z=1]
Erwarteter Wert des Interventionseffekts in der Probe, die die Intervention erhalten hat
In diesem Buch wird es durch Matching berechnet
ATE(Average Treatment Effect):
- E[Y^{(1)}-Y^{(0)}]
Erwarteter Wert des Interventionseffekts für die gesamte Probe
In diesem Buch wird es von IPW berechnet

Matching gegen IPW [^ Matching-vs-IPw]

[^ Matching-vs-IPW]: Es gibt andere Methoden zur Abschätzung des Interventionseffekts unter Verwendung anderer Neigungswerte als Matching und IPW, aber dieses Buch befasst sich mit diesen beiden.

** Matching **:
Nehmen Sie ein Paar Proben mit ähnlichen Neigungswerten aus der Interventionsgruppe und der Nicht-Interventionsgruppe und ergänzen Sie eine fehlende Maßnahme (Fälschung) durch die andere
IPW(Inverse Probability Weighting):
Gewichten Sie die Stichprobe mit der Umkehrung des Neigungswerts, um der Verteilung der Bevölkerung näher zu kommen --IPW berücksichtigt keine allgemeine Unterstützung (später beschrieben) und sollte nicht verwendet werden? [^ tsugawa-ps2] ――Ich denke, Sie können es verwenden, wenn Sie gemeinsame Unterstützung in Betracht ziehen

ATT im Matching

Abgleichen von Stichproben mit ähnlichen Neigungsbewertungswerten wie die Stichproben in der Interventionsgruppe aus der Nicht-Interventionsgruppe, wobei der Durchschnitt der Paardifferenzen ATT sein soll.

ATE im Matching

Ähnlich wie bei Stichproben in der Gruppe ohne Intervention wird der Durchschnitt der Differenz zwischen den Paaren als ATE verwendet.

ATT in IPW

Berechnet nach folgender Formel

\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

Berechnet nach folgender Formel

\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})}}

Über gemeinsame Unterstützung

Gemeinsame Unterstützung ist der Bereich, in dem sich die Verteilung der Neigungswerte der Interventionsgruppe und der Nicht-Interventionsgruppe überschneiden.
Wenn Sie eine Stichprobe verwenden, die bei der Berechnung des IPW nicht gemeinsam unterstützt wird, wird die inverse Zahl (Gewicht) des Neigungswerts zu groß, um eine gute Schätzung vorzunehmen [^ tsugawa-ps1] -Ist "aus gemeinsamer Unterstützung" und "Gewicht zu groß" gleichwertig?
Ist es üblich, Proben auszuschließen, die nicht allgemein unterstützt werden?
In diesem Buch wird das Trimmen / Abschneiden von Immobilienwerten vorgestellt ――Ist es nicht möglich, eine neue Tendenz zu schaffen, indem nicht gemeinsame Unterstützung ausgeschlossen wird? ――So ist es wichtig, den Versuchsplan so zu gestalten, dass Sie zwei vergleichbare Gruppen richtig beobachten können, damit Sie ihn nicht ausschließen müssen ** ――Die dabei erzeugte Verzerrung sollte durch die Tendenzbewertung usw. angepasst werden. ――Da es keine Stichprobe gibt, die die Gegentatsache außerhalb der gemeinsamen Unterstützung ergänzen kann, erfüllt sie nicht überhaupt die „stark vernachlässigbare Zuteilung“?
- P(Z=1|X)=0OderP(Z=1|X)=1Nicht, weil es logistische Regression verwendet

Bewertung der Neigungsbewertung

――Es ist wichtig, ob die Verteilung der Kovariaten zwischen der Interventionsgruppe und der Nicht-Interventionsgruppe zwischen den beiden Gruppen durch den Neigungswert ausgeglichen wird. -Machen Sie also Annahmen\\{Y^{(1)}, Y^{(0)}\\}\perp Z|e(X)Von\\{Y^{(1)}, Y^{(0)}\\}\perp Z|XIch denke, es bedeutet, zurückzukehren

In diesem Buch wird die standardisierte durchschnittliche Differenz (ASAM) zur Bewertung verwendet. --OK wenn ASAM 0.1 oder weniger ist (warum?) ――Wenn die alte AUC 0,8 oder mehr beträgt, ist die OK-Theorie Unsinn? [^ ps-auc]
Es ist in Ordnung, wenn die AUC 0,8 oder mehr beträgt [^ iwanami-ds3] --Was bedeutet es, dass die AUC klein ist (nahe 0,5):
Die Zuordnung hängt nicht von den Kovariaten ab ($ Z \ perp X $) → Keine Anpassung an den Neigungswert erforderlich → Ist es bedeutungslos, wenn die AUC zu einem gewissen Grad nicht groß ist? --Was bedeutet es, eine große AUC (nahe 1) zu haben: ――Die gemeinsame Unterstützung ist gering → Der Interventionseffekt kann nicht geschätzt werden → Ist es nicht gut, wenn die AUC groß ist?

Praktische Ausgabe

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

Wird zu einem späteren Zeitpunkt hinzugefügt

LaLonde(1986)

Zur Erläuterung der Daten wird auf dieses Handbuch verwiesen
Siehe hier für die Erklärung des Papiers.

Daten gelesen

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')

Datensatzerstellung

--NSW: RCT-Datensatz --CPS1 + NSW: Datensatz mit NSW-Nichtinterventionsgruppe durch CPS1 ersetzt --CPS3 + NSW: Datensatz mit NSW-Nichtinterventionsgruppe durch CPS3 ersetzt

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)

Bestätigen Sie die Verteilung der Interventionsgruppe und der Nicht-Interventionsgruppe

NSW

Die Daten sind ausgewogen, da sie nur RCT singen.

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

――Es handelt sich um ziemlich unausgeglichene Daten

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

Nicht so gut wie NSW, aber ausgewogener als 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]))

Regressionsanalyse

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]

Ergebnis

Datensatz	Interventionseffekt[$]
NSW	1676.3
CPS1+NSW	699.1
CPS3+NSW	1548.2

Spielbeginn

Der Interventionseffekt von NSW beträgt ca. 1600 US-Dollar
Fordern Sie heraus, wie nahe Sie mit dem Neigungswert aus einem voreingenommenen Zustand an 1600 US-Dollar kommen können

Schätzung der Neigungsbewertung

Woher kommen das Ich (re74 ^ 2) und das Ich (re75 ^ 2) in diesem Buch?

CPS1+NSW

Schätzen Sie den Neigungswert durch logistische Regression

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))

--Überprüfen Sie die Verteilung der Neigungswerte

(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))

AUC berechnen

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

0.9708918648513447

CPS3+NSW

Schätzen Sie den Neigungswert durch logistische Regression

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))

--Überprüfen Sie die Verteilung der Neigungswerte

(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))

AUC berechnen

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

0.8742266742266742

ATT-Schätzung

CPS1+NSW

passend

Definieren Sie eine Funktion (rekursiver Generator), die ein Paar erstellt
Hier wurde als Referenz verwendet.

from sklearn.neighbors import NearestNeighbors
def pairs_generator(s0: pd.Series, s1: pd.Series, threshold: np.float) -> pd.DataFrame:
    """
Ein rekursiver Generator, der Paare mit der K-Nachbarschaftsmethode erstellt.
    
    Parameters
    ----------
    s0, s1 : pd.Series
    threshold : np.float
Die Schwelle der Entfernung, um ein Paar zu bilden. Wenn der nächste Nachbarpunkt ebenfalls den Schwellenwert überschreitet, wird er nicht gepaart.
    
    Returns
    -------
    pairs : pd.DataFrame
Index eines Paares bestehend aus der s0-Probe und der nächsten s1-Probe.
Spalte l ist der Index von s0 und Spalte r ist der Index von s1.
    
    Notes
    -----
Dieselbe Stichprobe wird niemals für mehrere Paare verwendet.
    
    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)

Bilden Sie ein Paar

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) #Wird bei der ATE-Berechnung verwendet

Bestätigen Sie zwei Interventionsgruppen --Orange ist eine weggeworfene Probe ――Ursprünglich ist es wünschenswert, ein Paar zu bilden, damit der in der Abbildung gezeigte Schnittpunkt nicht auftritt (glaube ich).

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)

ATT berechnen
Die Standardabweichung ist ziemlich groß

(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

Gewicht berechnen

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

ATT mit gewichteter linearer Regression schätzen

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 (gemeinsame Unterstützung)

Extrahieren Sie Stichproben in dem Bereich, in dem sich die Verteilungen der Neigungsbewertungen überschneiden

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()

ATT mit gewichteter linearer Regression schätzen

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

Bestätigen Sie die standardisierte durchschnittliche Differenz zwischen den beiden Gruppen von Kovariaten ――Nach der Einstellung liegt sie im Allgemeinen unter der rot gepunkteten Linie.
Entspricht es der Grafik in diesem Buch? (Besonders schwarz vor der Einstellung)

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')

Bestätigen Sie die Verteilung der Kovariaten beim Matching
Der Durchschnitt sagt nicht die Wahrheit, schauen Sie sich also die Verteilung an ――Der Ausgleich erfolgt ordnungsgemäß im Vergleich zur ursprünglichen Verteilung

_, 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]))

Bestätigen Sie die Verteilung der Kovariaten in IPW ――Der Ausgleich erfolgt ordnungsgemäß im Vergleich zur ursprünglichen Verteilung

_, 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]))

--Bestätigen Sie die Verteilung der Kovariaten in IPW (gemeinsame Unterstützung) ――Der Ausgleich erfolgt ordnungsgemäß im Vergleich zur ursprünglichen Verteilung

_, 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

passend

Bilden Sie ein Paar

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) #Wird bei der ATE-Berechnung verwendet

Bestätigen Sie zwei Interventionsgruppen

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)

ATT berechnen

(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

Gewicht berechnen

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

ATT mit gewichteter linearer Regression schätzen

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 (gemeinsame Unterstützung)

Extrahieren Sie Stichproben in dem Bereich, in dem sich die Verteilungen der Neigungsbewertungen überschneiden

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()

ATT mit gewichteter linearer Regression schätzen

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

Bestätigen Sie die standardisierte durchschnittliche Differenz zwischen den beiden Gruppen von Kovariaten

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')

Bestätigen Sie die Verteilung der Kovariaten beim Matching

_, 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]))

Bestätigen Sie die Verteilung der Kovariaten in 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]))

--Bestätigen Sie die Verteilung der Kovariaten in IPW (gemeinsame Unterstützung)

_, 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]))

Ergebnis

――Die Kovariaten sind ausgeglichen, aber der geschätzte Wert des Interventionseffekts scheint sich deutlich von 1600 USD zu unterscheiden.

Datensatz	passend	IPW	IPW(Gemeinsame Unterstützung)
CPS1+NSW	1929.1	1180.4	1243.9
CPS3+NSW	1463.1	1214.1	988.9

ATE-Schätzung

CPS1+NSW

passend

--Überprüfen Sie das Paar

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)

ATE berechnen

(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

ATE mit gewichteter linearer Regression schätzen

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 (gemeinsame Unterstützung)

ATE mit gewichteter linearer Regression schätzen

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

Bestätigen Sie die standardisierte durchschnittliche Differenz zwischen den beiden Gruppen von Kovariaten

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')

Bestätigen Sie die Verteilung der Kovariaten beim Matching

_, 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]))

Bestätigen Sie die Verteilung der Kovariaten in 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]))

--Bestätigen Sie die Verteilung der Kovariaten in IPW (gemeinsame Unterstützung)

_, 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

passend

--Überprüfen Sie das Paar

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)

ATE berechnen

(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

ATE mit gewichteter linearer Regression schätzen

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 (gemeinsame Unterstützung)

ATE mit gewichteter linearer Regression schätzen

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

Bestätigen Sie die standardisierte durchschnittliche Differenz zwischen den beiden Gruppen von Kovariaten

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')

Bestätigen Sie die Verteilung der Kovariaten beim Matching

_, 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]))

Bestätigen Sie die Verteilung der Kovariaten in 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]))

--Bestätigen Sie die Verteilung der Kovariaten in IPW (gemeinsame Unterstützung)

_, 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]))

Ergebnis

Datensatz	passend	IPW	IPW(Gemeinsame Unterstützung)
CPS1+NSW	1836.5	-6456.3	545.3
CPS3+NSW	1486.6	224.7	801.0

Impressionen

Ich kann die "Voreingenommenheit" überhaupt nicht loswerden ...

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

"Einführung in die Effektüberprüfung Kapitel 3 Analyse mit dem Neigungswert" + α wird in Python versucht

Grundlagen

Was ist ein Neigungswert?

Regressive Analyse gegen Propensity Score [^ ra-vs-ps]

CIA gegen stark vernachlässigbare Zuordnung

Schätzung der Neigungsbewertung

Schätzung des Interventionseffekts anhand des Neigungsscores

Matching gegen IPW [^ Matching-vs-IPw]

ATT im Matching

ATE im Matching

Über gemeinsame Unterstützung

Bewertung der Neigungsbewertung

Praktische Ausgabe

Daten gelesen

Datensatzerstellung

Bestätigen Sie die Verteilung der Interventionsgruppe und der Nicht-Interventionsgruppe

Regressionsanalyse

Ergebnis

Spielbeginn

Schätzung der Neigungsbewertung

ATT-Schätzung

passend

IPW (gemeinsame Unterstützung)

passend

IPW (gemeinsame Unterstützung)

Ergebnis

ATE-Schätzung

passend

IPW (gemeinsame Unterstützung)

passend

IPW (gemeinsame Unterstützung)

Ergebnis

Impressionen