Mögliche Ergebnisse (potenzielle Ergebnisse) Hinweis auf kausale Inferenz in Python Teil 1
Dies ist eine japanische Übersetzung von diesem Notizbuch.
from __future__ import division
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
sns.set_palette("colorblind")
%matplotlib inline
import datagenerators as dg
import warnings
warnings.filterwarnings('ignore')
In diesem Beitrag ein Framework [potenzielles Ergebnis]( potenzielles Ergebnis]), das versucht, mithilfe des hervorragenden Pakets CausalInference auf kausale Zusammenhänge für Situationen zu schließen, in denen nur Beobachtungsdaten vorliegen. Dies ist eine Übersicht über die Verwendung von https://en.wikipedia.org/wiki/Rubin_causal_model. Der Autor hat eine Reihe von Blog-Posts über seine Funktionen geschrieben.
Die meisten Datensätze, die heruntergeladen werden können, sind statisch. In diesem Beitrag werde ich meine Funktion zum Generieren der Daten verwenden. Dies hat zwei Vorteile. Die Fähigkeit, Datensätze mit bestimmten Eigenschaften zu generieren und zu generieren, und die Fähigkeit, direkt in das Datengenerierungssystem einzugreifen, um zu überprüfen, ob die Inferenz korrekt ist. Alle diese Datengeneratoren generieren i.i.d. Samples aus mehreren Distributionen und geben die Ergebnisse als Pandas-Datenrahmen zurück. Die Funktionen, die diese Datensätze generieren, finden Sie in der angehängten Datei datagenerators.py auf github hier.
Schauen wir uns zunächst ein Motivationsbeispiel an.
Einführung
Eines Tages bemerkte ein Teamleiter, dass einige Mitglieder des Teams coole Hüte trugen, was die Produktivität dieser Teammitglieder tendenziell verringerte. Basierend auf den Daten basiert der Teamleiter darauf, ob Teammitglieder coole Hüte tragen ($ X = 1 $ für coole Hüte, $ X = 0 $ für nicht coole Hüte) und ob sie produktiv sind ($). Aufnahme starten (Y = 1 $ ist produktiv, $ Y = 0 $ ist unproduktiv).
Nach einer Woche Beobachtung erhielt das Team den folgenden Datensatz:
observed_data_0 = dg.generate_dataset_0()
observed_data_0.head()
| x | y | |
|---|---|---|
| 0 | 0 | 0 |
| 1 | 1 | 1 |
| 2 | 0 | 0 |
| 3 | 0 | 1 |
| 4 | 0 | 0 |
Die erste Frage, die der Teamleiter stellt, lautet: "Ist eine Person, die einen coolen Hut trägt, produktiver als jemand, der dies nicht tut?" Dies bedeutet die Schätzung der folgenden Mengen:
Und Sie können es direkt aus den Daten abschätzen:
def estimate_uplift(ds):
"""
Estiamte the difference in means between two groups.
Parameters
----------
ds: pandas.DataFrame
a dataframe of samples.
Returns
-------
estimated_uplift: dict[Str: float] containing two items:
"estimated_effect" - the difference in mean values of $y$ for treated and untreated samples.
"standard_error" - 90% confidence intervals arround "estimated_effect"
"""
base = ds[ds.x == 0]
variant = ds[ds.x == 1]
delta = variant.y.mean() - base.y.mean()
delta_err = 1.96 * np.sqrt(
variant.y.var() / variant.shape[0] +
base.y.var() / base.shape[0])
return {"estimated_effect": delta, "standard_error": delta_err}
estimate_uplift(observed_data_0)
{'estimated_effect': -0.15738429037833684,
'standard_error': 0.08639322378194732}
Leute mit coolen Hüten scheinen weniger produktiv zu sein.
Sie können auch einen statistischen Test durchführen, um dies sicherzustellen.
from scipy.stats import chi2_contingency
contingency_table = (
observed_data_0
.assign(placeholder=1)
.pivot_table(index="x", columns="y", values="placeholder", aggfunc="sum")
.values
)
_, p, _, _ = chi2_contingency(contingency_table, lambda_="log-likelihood")
# p-value
p
0.0005626147113456373
Dies ist ein kleiner p-Wert. Statistiker werden stolz sein
Sie können diese Informationen verwenden, um zu erklären, was Sie über die Wahrscheinlichkeit denken, dass jemand einen coolen Hut trägt. Solange wir glauben, dass es "aus derselben Verteilung abgeleitet" ist wie frühere Beobachtungen, wird erwartet, dass dieselbe Korrelation besteht.
Das Problem ist, wenn Teamleiter versuchen, diese Informationen als Debatte darüber zu verwenden, ob Menschen gezwungen werden sollen, coole Hüte zu tragen. Wenn ein Teamleiter dies tut, kann er das System, das wir abtasten, radikal ändern, ändern oder umgekehrt von zuvor beobachteten Korrelationen.
Der sauberste Weg, um die Auswirkungen einiger Änderungen in Ihrem System tatsächlich zu messen, besteht darin, einen randomisierten Vergleichstest durchzuführen (https://en.wikipedia.org/wiki/Randomized_controlled_trial). Stellen Sie insbesondere fest, wer einen coolen Hut bekommt und wer nicht, und sehen Sie den Unterschied im Wert $ y $, den Sie erhalten. Dadurch werden die Auswirkungen von Interlaced-Variablen (https://en.wikipedia.org/wiki/Confounding) beseitigt, die sich möglicherweise auf die Metriken auswirken, die Sie interessieren.
Nachdem wir den Datensatz aus einem bekannten Prozess (in diesem Fall der von mir erstellten Funktion) generiert haben, können wir direkt eingreifen, um die Wirksamkeit des A / B-Tests zu messen.
def run_ab_test(datagenerator, n_samples=10000, filter_=None):
"""
Generates n_samples from datagenerator with the value of X randomized
so that 50% of the samples recieve treatment X=1 and 50% receive X=0,
and feeds the results into `estimate_uplift` to get an unbiased
estimate of the average treatment effect.
Returns
-------
effect: dict
"""
n_samples_a = int(n_samples / 2)
n_samples_b = n_samples - n_samples_a
set_X = np.concatenate([np.ones(n_samples_a), np.zeros(n_samples_b)]).astype(np.int64)
ds = datagenerator(n_samples=n_samples, set_X=set_X)
if filter_ != None:
ds = ds[filter_(ds)].copy()
return estimate_uplift(ds)
run_ab_test(dg.generate_dataset_0)
{'estimated_effect': 0.18679999999999997,
'standard_error': 0.019256613018207393}
Plötzlich scheint sich die Richtung der Wirkung des Tragens eines kühlen Hutes umgekehrt zu haben.
Was ist los?
Hinweis: Im obigen Beispiel und in allen nachfolgenden Beispielen lautet das Beispiel iid und SUTVA. Es wird davon ausgegangen, dass Sie / wiki / Rubin_causal_model # Stable_unit_treatment_value_assumption_% 28SUTVA% 29 folgen. Grundsätzlich bedeutet dies, dass eine Person, die einen wirklich coolen Hut wählt oder gezwungen ist, diesen zu tragen, die Entscheidungen oder Auswirkungen einer anderen Person, die wirklich einen coolen Hut trägt, nicht beeinflusst. Ich werde. Strukturell haben alle von mir verwendeten synthetischen Datengeneratoren diese Eigenschaft. In Wirklichkeit muss man davon ausgehen, dass es wahr ist.
Definition des Kausalzusammenhangs
Das vorige Beispiel zeigt die alte statistische Proklamation:
** Korrelation bedeutet nicht Kausalzusammenhang **
"Kausale Beziehung" ist ein vages und philosophisches Wort. Im aktuellen Kontext bedeutet dies: "Wie wirkt sich die Änderung von $ X $ auf $ Y $ aus?"
Um genau zu sein, sind $ X $ und $ Y $ Zufallsvariablen, und der "Effekt", den Sie wissen möchten, besteht darin, $ X $ einen bestimmten Wert zuzuweisen. Dann, wie sich die Verteilung von $ Y $ ändert. Dieser Vorgang, bei dem eine Variable gezwungen wird, einen bestimmten Wert zu haben, wird als "Intervention" bezeichnet.
Wenn wir im vorherigen Beispiel nicht in das System eingreifen würden, würden wir eine Beobachtungsverteilung von $ Y $ erhalten, vorausgesetzt, wir haben $ X $ beobachtet.
Wir greifen ein, wenn wir die Leute dazu zwingen, cool auszusehen. Und die Verteilung von $ Y $ ist durch die Interventionsverteilung gegeben.
Im Allgemeinen sind die beiden nicht gleich.
Die Frage, die diese Notizen zu beantworten versuchen, ist, wie wir auf die Interventionsverteilung schließen können, wenn nur Beobachtungsdaten verfügbar sind. Dies ist eine nützliche Frage, da es viele unrealistische, machbare oder unethische Situationen gibt, in denen es unpraktisch, machbar oder unethisch ist, einen A / B-Test durchzuführen, um die Wirksamkeit einer Intervention direkt zu messen. Auch in dieser Situation möchte ich etwas über die Wirkung der Intervention sagen können.
Mögliches Ergebnis
Eine Möglichkeit, dieses Problem anzugehen, besteht darin, zwei neue Zufallsvariablen in das System einzuführen: $ Y_ {0} $ und $ Y_ {1} $, potenzielle Ergebnisse (http: // www. Bekannt als stat.unipg.it/stanghellini/rubinjasa2005.pdf). Wir gehen davon aus, dass diese Variablen existieren und wie jede andere Zufallsvariable behandelt werden können.
$ Y $ ist wie folgt definiert:
Y = Y_{1} whenX=1 Y = Y_{0} whenX=0
Dies ist eine zugrunde liegende Verteilung mit fehlenden Werten (https://en.wikipedia.org/wiki/Missing_data), aus der hervorgeht, wie sich die Verteilung ändert, wenn in das Problem eingegriffen wird. Verschiebt sich von zu ungefähr zu den mit iid gezeichneten Daten. Unter bestimmten Annahmen darüber, warum Werte fehlen, gibt es gut entwickelte Theorien darüber, wie fehlende Werte geschätzt werden können.
Ziel
Oft ist uns die vollständige Interventionsverteilung $ P (Y | \ hbox {do} (X)) $ egal, und eine Schätzung der Differenz zwischen den Durchschnittswerten zwischen den beiden Gruppen ist ausreichend. Dies ist eine Zahl, die als Average Treatment Effect (ATE) bezeichnet wird (https://en.wikipedia.org/wiki/Average_treatment_effect).
Wenn Sie A / B-Tests durchführen und die Durchschnittswerte für jede Gruppe vergleichen, hängt dies direkt mit den Zahlen zusammen, die wir messen.
Das Schätzen dieses Wertes aus der beobachteten Verteilung ergibt Folgendes:
Dies entspricht im Allgemeinen nicht der tatsächlichen ATE. weil
Die zwei verwandten Größen
ATT = E[Y_{1} - Y_{0}|X=1] , the "Average Treatment effect of the Treated"ATC = E[Y_{1} - Y_{0}|X=0] , the "Average Treatment effect of the Control"
Deshalb.
Eine Möglichkeit zur Interpretation von ATC besteht darin, den Effekt der Verarbeitung nur unbehandelter Proben zu messen und umgekehrt für ATT. Abhängig vom Anwendungsfall können sie ein natürlicheres Maß dafür sein, was Sie interessiert. Sie können sie alle mithilfe der folgenden Techniken schätzen:
Vermutungen aufstellen
Randomisieren Sie beim A / B-Test die Zuweisung von $ X $ nach dem Zufallsprinzip. Auf diese Weise können Sie auswählen, ob die Variablen $ Y_ {1} $ oder $ Y_ {0} $ angezeigt werden sollen. Dies macht das Ergebnis unabhängig vom Wert von $ X $.
Schreiben Sie dies wie folgt:
Das heißt, die Verteilung von $ X, Y_ {0}, Y_ {1} $ wird wie folgt berücksichtigt:
Wenn diese Unabhängigkeit garantiert ist, dann:
Wenn Sie die ATE anhand von Beobachtungsdaten schätzen möchten, müssen Sie andere Informationen zu der Stichprobe verwenden. Insbesondere muss davon ausgegangen werden, dass genügend zusätzliche Informationen vorhanden sind, um die Behandlungsoptionen für jedes Subjekt vollständig zu erläutern.
Wenn man diese zusätzlichen Informationen als Zufallsvariable $ Z $ bezeichnet, kann diese Annahme wie folgt geschrieben werden:
Oder
Dies bedeutet, dass die beobachtete Behandlung $ X $, die die Probe erhält, vollständig durch $ Z $ erklärt wird. Dies wird manchmal als "Annahme ignorieren" bezeichnet (https://en.wikipedia.org/wiki/Ignorability).
In diesem Beispiel ist das Beispiel für die Motivation für einen coolen Hut, dass es andere Faktoren gibt (nennen wir es "Fähigkeiten"), die sowohl die menschliche Produktivität als auch die Frage beeinflussen, ob sie einen coolen Hut tragen oder nicht. Meint Im obigen Beispiel ist ein Fachmann produktiver und trägt weniger wahrscheinlich einen coolen Hut. Diese Tatsachen könnten zusammen erklären, warum sich die Auswirkungen des Kühlhutes beim Ausführen des A / B-Tests umzukehren schienen.
Wenn Sie die Daten danach teilen, ob die Person kompetent ist oder nicht, können Sie feststellen, dass in jeder Untergruppe ein positiver Zusammenhang zwischen dem Tragen eines coolen Hutes und der Produktivität besteht.
observed_data_0_with_confounders = dg.generate_dataset_0(show_z=True)
print(estimate_uplift(observed_data_0_with_confounders.loc[lambda df: df.z == 0]))
print(estimate_uplift(observed_data_0_with_confounders.loc[lambda df: df.z == 1]))
{'estimated_effect': 0.1647058823529412, 'standard_error': 0.15924048578984673}
{'estimated_effect': 0.14041297935103236, 'standard_error': 0.17328031714939449}
Leider können wir $ Y_ {0} $ und $ Y_ {1} $ in derselben Stichprobe nicht beobachten, sodass wir die folgenden Annahmen nicht testen können.
Es muss anhand des Wissens des Systems bewertet werden, das wir untersuchen.
Die Qualität der Vorhersagen, die Sie treffen, hängt davon ab, wie gut diese Annahme gehalten wird. Das Simpson-Paradoxon (http://www.degeneratestate.org/posts/2017/Oct/22/generating-examples-of-simpsons-paradox/) besagt, dass $ Z $ nicht alle verschränkten Variablen enthält. Ein extremes Beispiel für die Tatsache, dass die Schlussfolgerungen, die wir ziehen, falsch sein können. Facebook veröffentlicht ein gutes Papier, in dem verschiedene kausale Inferenzansätze mit direkten A / B-Tests verglichen werden. Es zeigt, wie die Wirksamkeit überschätzt wird, wenn die bedingte Unabhängigkeit nicht gewahrt bleibt. Ich bin.
Sobald Sie diese Annahme getroffen haben, gibt es verschiedene Techniken, um sich ihr zu nähern. Der Rest dieses Artikels wird einige der einfacheren Ansätze skizzieren, aber denken Sie daran, dass dies ein fortlaufender Forschungsbereich ist.
Modellierung von Gegenfakten
Aus dem Obigen ist klar, dass ATE geschätzt werden kann, wenn $ Y_ {0} $ und $ Y_ {1} $ bekannt sind. Warum also nicht versuchen, diese direkt zu modellieren?
Insbesondere können Sie die folgenden Schätzungen vornehmen.
\hat{Y}_{0}(Z) = E[Y|Z, X=0] \hat{Y}_{1}(Z) = E[Y|Z, X=1]
Wenn wir diese beiden Größen modellieren können, können wir die ATE wie folgt schätzen.
Der Erfolg dieses Ansatzes hängt davon ab, wie gut die potenziellen Ergebnisse modelliert werden können. Um es in Aktion zu sehen, verwenden wir den folgenden Datengenerierungsprozess.
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False);
Bevor wir uns mit der Modellierung der Gegenfakten befassen, werfen wir einen Blick auf die Daten. Wenn man sich ansieht, wie $ Y $ verteilt ist, sieht es so aus, als gäbe es einen kleinen Unterschied zwischen den beiden Gruppen.
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 0].y, label="untreated")
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 1].y, label="treated")
<matplotlib.axes._subplots.AxesSubplot at 0x1a1a6c3590>
Dies kann durch Betrachtung des durchschnittlichen Unterschieds zwischen den beiden Gruppen bestätigt werden.
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
Observed ATE: 0.113 (0.103)
Wenn Sie sich jedoch die Verteilung von $ Z $ ansehen, die die Kovarianz darstellt, können Sie feststellen, dass es Unterschiede zwischen den Gruppen gibt.
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 0].z, label="untreated")
sns.kdeplot(observed_data_1.loc[lambda df: df.x == 1].z, label="treated")
<matplotlib.axes._subplots.AxesSubplot at 0x1a1a7bced0>
Dies ist besorgniserregend, wenn Sie der Meinung sind, dass $ Z $ einen gewissen Einfluss auf $ Y $ hat. Wir müssen einen Weg finden, um zwischen den Auswirkungen von $ X $ auf $ Y $ und den Auswirkungen von $ Z $ auf $ Y $ zu unterscheiden.
Sie können die tatsächliche ATE mit einem A / B-Test überprüfen, der sie simuliert, und bestätigen, dass es sich um die Differenz zum beobachteten Wert handelt.
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Real ATE: -0.479 (0.026)
Aber was ist, wenn Sie diesen A / B-Test nicht ausführen können? Sie müssen Ihr System modellieren.
Der einfachste Modelltyp ist das lineare Modell. Nehmen Sie insbesondere Folgendes an:
Wenn dies korrekt ist, lassen Sie das Modell zu den Daten passen
Sie können die lineare Regression verwenden, um eine Schätzung der ATE zu erhalten.
Das "Causalinference" -Paket bietet hierfür eine einfache Schnittstelle.
from causalinference import CausalModel
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_via_ols(adj=1)
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.488 0.032 -15.122 0.000 -0.551 -0.425
"Kausalinferenz" gibt eine Schätzung von ATE und die statistischen Eigenschaften dieser Schätzung zurück. Das für die Schätzung angegebene Konfidenzintervall ist das Konfidenzintervall unter der Annahme, dass das Modell das Kontrafakt genau beschreibt, und das Konfidenzintervall, wie gut das Modell das Kontrafakt beschreibt. Es ist wichtig zu erkennen, dass es so etwas nicht gibt.
In diesem Fall ist es dem Paket gelungen, die richtige ATE zu identifizieren - was gut ist, aber der Datengenerierungsprozess ist speziell darauf ausgelegt, die Annahmen zu erfüllen. Schauen wir uns einige Fälle an, die fehlschlagen können.
Der erste Fall ist, wenn der Effekt nicht einfach additiv ist.
observed_data_2 = dg.generate_dataset_2()
observed_data_2.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_2)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_2)))
Observed ATE: 0.689 (0.104)
Real ATE: 0.563 (0.031)
cm = CausalModel(
Y=observed_data_2.y.values,
D=observed_data_2.x.values,
X=observed_data_2.z.values)
cm.est_via_ols(adj=1)
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 0.246 0.088 2.805 0.005 0.074 0.419
Dies kann normalerweise durch stärkere Schätzungen überwunden werden. Ein einfacher, nicht parametrischer Ansatz ist die Matching-Methode (https://en.wikipedia.org/wiki/Matching_%28statistics%29). Die Idee ist, für jede behandelte Probe ähnliche unbehandelte Proben zu finden und diese Werte direkt zu vergleichen. Was Sie unter "ähnlich" verstehen, hängt von Ihrem speziellen Anwendungsfall ab.
Das "Kausalferenz" -Paket implementiert das Matching für jede Einheit, indem die ähnlichste Einheit aus den anderen Behandlungsgruppen durch Ersetzen ausgewählt und die Differenz zwischen diesen beiden Einheiten zur Berechnung der ATE verwendet wird. .. Standardmäßig werden Übereinstimmungsentscheidungen unter Verwendung der Methode des nächsten Nachbarn im kovariaten Raum $ Z $ getroffen, und Abstände werden mit der inversen Varianz jeder Dimension gewichtet.
Sie haben die Möglichkeit, die Anzahl der zu vergleichenden Einheiten und die Gewichtung jeder Dimension in der Übereinstimmung zu ändern. Weitere Informationen finden Sie in der Dokumentation (http://laurence-wong.com/software/matching).
Die übereinstimmende Schätzung kann mit dem folgenden Code berechnet werden.
cm = CausalModel(
Y=observed_data_2.y.values,
D=observed_data_2.x.values,
X=observed_data_2.z.values)
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 0.480 0.148 3.239 0.001 0.190 0.770
ATC 0.973 0.142 6.854 0.000 0.695 1.251
ATT 0.036 0.211 0.169 0.866 -0.379 0.450
Das Konfidenzintervall der Schätzung enthält die wahre ATE.
Covariate Ungleichgewicht
Ein schwieriger zu behandelndes Problem besteht darin, dass die verwendeten Kovariaten unausgeglichen sind: Der Kovariatenraum hat einen Bereich, der nur verarbeitete oder unverarbeitete Proben enthält. Hier müssen wir den Effekt der Behandlung extrapolieren - dies hängt stark von dem hypothetischen Modell ab, das wir verwenden.
Das folgende Beispiel zeigt dies.
observed_data_3 = dg.generate_dataset_3()
observed_data_3.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_3)))
Observed ATE: 1.362 (0.080)
Real ATE: 2.449 (0.033)
# OLS estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_via_ols()
print(cm.estimates)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.969 0.055 35.707 0.000 1.861 2.078
ATC 1.971 0.056 35.117 0.000 1.861 2.081
ATT 1.968 0.073 26.923 0.000 1.825 2.112
# Matching estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.913 0.156 12.269 0.000 1.607 2.218
ATC 1.965 0.277 7.093 0.000 1.422 2.509
ATT 1.876 0.177 10.626 0.000 1.530 2.222
Die OLS-Schätzung erfasst nicht den tatsächlichen Effekt, und die Übereinstimmungsschätzung verbessert sich ein wenig, aber die Daten enthalten nicht genügend Informationen, um vollständig auf einen nicht überlappenden Bereich zu extrapolieren.
Dieses Beispiel mag inkonsistent erscheinen, aber sobald Sie sich mit höherdimensionalen Kovariaten befassen, tritt das Problem häufiger auf.
causalinference bietet ein praktisches Tool zum schnellen Auswerten doppelter Variablen mithilfe der Eigenschaft 'summary_stats`.
print(cm.summary_stats)
Summary Statistics
Controls (N_c=206) Treated (N_t=294)
Variable Mean S.d. Mean S.d. Raw-diff
--------------------------------------------------------------------------------
Y -0.151 0.438 1.210 0.467 1.362
Controls (N_c=206) Treated (N_t=294)
Variable Mean S.d. Mean S.d. Nor-diff
--------------------------------------------------------------------------------
X0 0.267 0.185 0.681 0.205 2.116
Hier ist die normalisierte Differenz definiert als:
Dies ist kein strenger statistischer Test, liefert jedoch einige Indikatoren dafür, wie stark sich die einzelnen Kovariaten überschneiden. Werte größer als 1 deuten darauf hin, dass nicht viel dupliziert wird.
Neigungsbewertung
Propensity Scores sind Schätzungen darüber, wie wahrscheinlich es ist, dass ein Subjekt bei einer Kovariate behandelt wird.
Sie können dies nach Belieben schätzen, aber sobald Sie es geschätzt haben, können Sie einige Dinge tun.
Gewichtung durch umgekehrte Tendenzbewertung
Das Problem bei der Messung der kausalen Inferenz besteht darin, dass ich die Größe $ E [Y_ {i}] $ wissen möchte, aber bedenke, dass ich nur auf eine Stichprobe von $ E [Y_ {i} | X = i] $ zugreifen kann. Bitte gib mir.
Die Wahrscheinlichkeit eines möglichen Ergebnisses kann wie folgt erweitert werden:
Das ist wahr
Schlägt vor, dass geschätzt werden kann.
Daher können Sie potenzielle Ergebnisse wiederherstellen, indem Sie jeden Punkt mit einem umgekehrten Neigungswert gewichten. Das Ergebnis ist die Reverse Trend Score Weighting-Methode (https://en.wikipedia.org/wiki/Inverse_probability_weighting).
Mal sehen, was mit einem der vorherigen Datensätze zu tun ist.
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Observed ATE: 0.147 (0.102)
Real ATE: -0.490 (0.025)
Sie können die Methoden "est_propensity_s" und "est_propensity" im "CausalInference" -Paket verwenden, um Trends mithilfe der logistischen Regression für Kovariaten zu schätzen.
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
propensity = cm.propensity["fitted"]
df = observed_data_1
df["ips"] = np.where(
df.x == 1,
1 / propensity,
1 / (1 - propensity))
df["ipsw"] = df.y * df.ips
ipse = (
df[df.x == 1]["ipsw"].sum()
- df[df.x == 0]["ipsw"].sum()
) / df.shape[0]
ipse
-0.4613297031604915
Mit dem in diesem Beispiel verwendeten Datengenerator können Sie die Beziehung mit einer einfachen logistischen Regression gut beschreiben. In dem in diesem Beispiel verwendeten Datengenerator kann die Beziehung durch eine einfache logistische Regression gut beschrieben werden. Wenn wir jedoch die logistische Regressionsfunktion von sklean verwenden (die Standardisierung wird standardmäßig verwendet). Wenn Sie versuchen, den Neigungswert zu schätzen, erhalten Sie die falsche Antwort.
from sklearn.linear_model import LogisticRegression
lg = LogisticRegression()
X = df.z.values.reshape(-1,1)
y = df.x.values
lg.fit(X,y)
propensity = lg.predict_proba(X)[:,1]
df["ips"] = np.where(
df.x == 1,
1 / propensity,
1 / (1 - propensity))
df["ipsw"] = df.y * df.ips
ipse = (
df[df.x == 1]["ipsw"].sum()
- df[df.x == 0]["ipsw"].sum()
) / df.shape[0]
ipse
-0.3291860760196773
Besser als unsere naiven Schätzungen, aber nicht korrekt.
Robuster gewichteter Schätzer
Auf diese Weise können Sie eine gewichtete Schätzung des inversen Neigungswerts mit einer linearen Schätzung des Effekts kombinieren, um zu versuchen, die Fehler in beiden Modellen zu reduzieren. Dies erfolgt durch Vorab-Durchführung einer gewichteten linearen Regression der Daten, wobei jeder Punkt mit einem Reverse-Propensity-Score gewichtet wird. Das Ergebnis ist doppelt robuster gewichteter Schätzer.
Die Idee ist, dass es eine Verzerrung gibt, bei der Proben in den Beobachtungsdaten behandelt werden, so dass die Proben, die behandelt wurden, aber weniger wahrscheinlich behandelt werden, wichtiger sind und mehr Gewichte haben. Sollte gegeben sein.
Es gibt die folgenden Methoden, um dies anzuwenden.
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Observed ATE: 0.145 (0.106)
Real ATE: -0.494 (0.026)
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
cm.est_via_weighting()
print(cm.estimates)
Treatment Effect Estimates: Weighting
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.521 0.034 -15.344 0.000 -0.588 -0.455
cm.estimates
{'weighting': {'ate': -0.5214252844257842, 'ate_se': 0.03398288940994425}}
Unbestimmtheit und Neigung Punktzahl
Im vorherigen Abschnitt haben wir angenommen, dass das Ergebnis und die Behandlung unabhängig sind, da die Kovariaten angegeben wurden.
Sie können auch eine etwas stärkere Annahme treffen. Das heißt, das Ergebnis ist unabhängig von der Behandlung und hängt von der Wahrscheinlichkeit der Neigung ab.
Mit dieser Annahme können Sie die Dimensionen der Verschränkungsvariablen reduzieren. Dies ermöglicht es, einige Techniken auszuführen, die in höherdimensionalen Einstellungen möglicherweise nicht funktionieren.
Trimmen
Wir haben zuvor bestätigt, dass kovariate Ungleichgewichte Probleme verursachen können. Eine einfache Lösung besteht darin, eine Kontrafaktvorhersage nur in Bereichen mit guter Überlappung durchzuführen oder Punkte zu "trimmen", an denen keine gute Überlappung vorliegt. Für hochdimensionale Daten kann es schwierig sein, "gute Überlappung" zu definieren. Die Verwendung von Neigungsbewertungen zur Definition von Überlappung ist eine Möglichkeit, dies zu lösen.
Schauen wir uns einen Datensatz mit weniger Überlappung an.
observed_data_3 = dg.generate_dataset_3()
observed_data_3.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_3)))
Observed ATE: 1.361 (0.080)
Real ATE: 2.437 (0.033)
CausalInference bietet eine Möglichkeit, Daten basierend auf Neigungswerten zu trimmen.
# OLS estimator
cm = CausalModel(
Y=observed_data_3.y.values,
D=observed_data_3.x.values,
X=observed_data_3.z.values)
cm.est_propensity_s()
cm.trim_s()
cm.est_via_matching()
print(cm.estimates)
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 2.068 0.127 16.350 0.000 1.820 2.316
ATC 2.249 0.229 9.843 0.000 1.802 2.697
ATT 1.920 0.102 18.739 0.000 1.719 2.121
Sie können den Rest der Daten wie folgt sehen:
# mask out data ignored by the
propensity = cm.propensity["fitted"]
cutoff = cm.cutoff
mask = (propensity > cutoff) & (propensity < 1 - cutoff)
# plot the data
observed_data_3[mask].plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False)
# actual response curves
z = np.linspace(0,1,100)
y0 = np.where(z >= 0.4, -4*(z - 0.4), 0)
y1 = np.where(z < 0.6, -4*(z - 0.6), 0) + 1
plt.plot(z,y0, "b")
plt.plot(z,y1, "r")
# filter out data in regions we have trimmed when we calculate the true uplift
filter_ = lambda df: (df.z > 0.2) & (df.z < 0.7)
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_3[mask])))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(
**run_ab_test(dg.generate_dataset_3, filter_= filter_)))
Observed ATE: 1.708 (0.125)
Real ATE: 1.984 (0.031)
Dies funktioniert in solchen Fällen ziemlich gut.
Bei der Anwendung des Trimmens wird ausdrücklich angegeben, dass eine kausale Folgerung nur für einige Stichproben des kovariaten Raums möglich ist. Über ATE für Proben außerhalb dieser Regionen kann nichts gesagt werden.
Schichtung
Eine andere Verwendung für Neigungsbewertungen sind geschichtete oder blockierte Schätzungen. Es besteht aus der Gruppierung von Datenpunkten in Gruppen mit ähnlichen Trends und der Schätzung der ATE innerhalb dieser Gruppen. Auch hier bietet CausalInference eine großartige Schnittstelle, um dies zu tun.
Verwenden Sie die Methode "stratify" (für benutzerdefinierte Statusgrenzen) oder die Methode "stratify_s" (Grenzen automatisch auswählen), um die Ebene zu bestimmen.
observed_data_1 = dg.generate_dataset_1()
cm = CausalModel(
Y=observed_data_1.y.values,
D=observed_data_1.x.values,
X=observed_data_1.z.values)
cm.est_propensity_s()
cm.stratify_s()
print(cm.strata)
Stratification Summary
Propensity Score Sample Size Ave. Propensity Outcome
Stratum Min. Max. Controls Treated Controls Treated Raw-diff
--------------------------------------------------------------------------------
1 0.093 0.260 108 18 0.163 0.182 -0.544
2 0.266 0.327 23 9 0.296 0.292 -0.639
3 0.327 0.411 14 17 0.368 0.376 -0.441
4 0.412 0.555 31 31 0.473 0.490 -0.428
5 0.557 0.769 45 80 0.661 0.667 -0.499
6 0.770 0.910 18 106 0.838 0.859 -0.536
Verwenden Sie dann die Methode "est_via_blocking", um die Schätzungen für diese Schichten zu einer Gesamt-ATE zu kombinieren.
cm.est_via_blocking()
print(cm.estimates)
Treatment Effect Estimates: Blocking
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.559 0.033 -16.722 0.000 -0.625 -0.494
ATC -0.571 0.039 -14.797 0.000 -0.646 -0.495
ATT -0.548 0.038 -14.304 0.000 -0.624 -0.473
Das funktioniert auch gut.
Das Aufteilen von Daten in Gruppen nach Neigungsbewertung ist nützlich, wenn Sie keine Vorkenntnisse darüber haben, was eine "ähnliche" Einheit ausmacht, aber dies ist nicht der einzige Weg. Wenn Sie im Voraus wissen, dass verschiedene Gruppen von Proben wahrscheinlich auf ähnliche Weise von Eingriffen betroffen sind, teilen Sie die Stichprobe vor der Schätzung der ATE in diese Gruppen ein und bündeln Sie die Ergebnisse global. Es ist sinnvoll, eine ATE zu bekommen.
Welche Methode soll verwendet werden?
Dies umfasst die meisten gängigen Methoden zur kausalen Folgerung aus beobachteten Daten. Die verbleibende Frage lautet: "Wie und wie soll entschieden werden, welche Methode verwendet werden soll?" Dies ist keine einfache Frage. Es gibt auch automatisierte Techniken wie dieses Dokument, aber ich habe sie nicht ausprobiert.
Um Ihre Technik zu wählen, müssen Sie letztendlich einige Annahmen darüber treffen, wie Sie die Gegenfakten aufbauen. Matching ist ein guter Ansatz, wenn Sie darauf vertrauen, dass sich die Daten im kovariaten Bereich gut überlappen. Wenn dies nicht der Fall ist, müssen Sie ein zuverlässiges Modell verwenden, um erfolgreich auf den unerforschten Bereich zu extrapolieren, oder die Annahme treffen, dass so etwas wie ein Neigungswert genügend Informationen erfasst, um eine Vernachlässigbarkeit anzunehmen. Es gibt.
Um zu betonen, dass alle diese Methoden fehlschlagen können, hier ein weiteres Beispiel. Im Gegensatz zum vorherigen Beispiel gibt es eine oder mehrere Kovariaten. Wie alle vorherigen Datengeneratoren folgt auch dieser Datengenerator den folgenden Annahmen:
Lassen Sie uns blind die Methoden ausprobieren, die wir bisher besprochen haben, und sehen, was passiert.
data_gen = dg.generate_exercise_dataset_2
ds = data_gen()
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(ds)))
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(data_gen)))
zs = [c for c in ds.columns if c.startswith("z")]
cm = CausalModel(
Y=ds.y.values,
D=ds.x.values,
X=ds[zs].values)
cm.est_via_ols()
cm.est_via_matching()
cm.est_propensity_s()
cm.est_via_weighting()
cm.stratify_s()
cm.est_via_blocking()
print(cm.estimates)
Observed ATE: -0.423 (0.711)
Real ATE: 4.570 (0.184)
Treatment Effect Estimates: OLS
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.165 0.363 -0.455 0.649 -0.878 0.547
ATC 0.438 0.402 1.088 0.276 -0.350 1.226
ATT -0.468 0.384 -1.218 0.223 -1.220 0.285
Treatment Effect Estimates: Matching
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 1.241 0.581 2.137 0.033 0.103 2.379
ATC 1.314 0.695 1.890 0.059 -0.049 2.677
ATT 1.204 0.667 1.804 0.071 -0.104 2.512
Treatment Effect Estimates: Weighting
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE 3.190 0.478 6.667 0.000 2.252 4.128
Treatment Effect Estimates: Blocking
Est. S.e. z P>|z| [95% Conf. int.]
--------------------------------------------------------------------------------
ATE -0.003 0.378 -0.009 0.993 -0.745 0.738
ATC -0.003 0.378 -0.009 0.993 -0.745 0.738
ATT -0.003 0.378 -0.009 0.993 -0.745 0.738
y = []
yerr = []
x_label = []
for method, result in dict(cm.estimates).items():
y.append(result["ate"])
yerr.append(result["ate_se"])
x_label.append(method)
x = np.arange(len(y))
plt.errorbar(x=x, y=y, yerr=yerr, linestyle="none", capsize=5, marker="o")
plt.xticks(x, x_label)
plt.title("Estimated Effect Size", fontsize=18)
plt.hlines(4.6, -0.5, 3.5, linestyles="dashed")
plt.xlim(-0.5,3.5);
Die gestrichelte Linie daneben zeigt die wahre ATE dieses Datensatzes.
Die Ergebnisse jeder Methode sind nicht nur unterschiedlich, sondern alle Methoden verfehlen den wahren Wert.
Dies sollte eine Warnung vor den Einschränkungen dieser Art von Technik sein. Es kann für den Leser eine sehr interessante Übung sein, herauszufinden, welche Eigenschaften des Datensatzes dazu führen, dass diese Techniken von ihren wahren Werten abweichen.
Struktur des kausalen Denkens
Hoffentlich sind Sie sich zu diesem Zeitpunkt der Bedeutung der Vernachlässigungsannahme bewusst.
Worüber ich bisher noch nicht gesprochen habe, ist, wie man $ Z $ wählt, damit dies wahr ist. Letztendlich sollte es aus der Kenntnis des Bereichs des untersuchten Systems stammen. Es gibt eine Reihe leistungsstarker Tools namens Causal Relationship Model (https://en.wikipedia.org/wiki/Causal_graph). Sie können es verwenden, um Ihr Wissen über das System, das Sie studieren, als grafisches Modell Ihres Systems zu kodieren oder um auf die oben beschriebenen Annahmen der bedingten Unabhängigkeit zu schließen.
Der nächste Beitrag zum kausalen Denken wird diese diskutieren.
Eine weitere Frage, die in diesem Beitrag aufgeworfen wird, ist, ob der einzige Weg, kausale Schlussfolgerungen zu ziehen, darin besteht, die Verschränkungsvariablen anzupassen. Nicht so - Sie können es an anderer Stelle in [Posts later] verwenden (http://www.degeneratestate.org/posts/2018/Sep/03/causal-inference-with-python-part-3-frontdoor-adjustment/). Werfen wir einen Blick auf die Technik.
Code
Die Hinweise zu diesem Beitrag finden Sie unter github hier.
Recommended Posts