Résultats potentiels (résultats potentiels) Note d'inférence causale dans Python Partie 1
Ceci est une traduction japonaise de ce cahier.
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')
Dans cet article, j'utilise l'excellent package CausalInference pour essayer de déduire des relations causales pour des situations où il n'y a que des données d'observation [ résultats potentiels]( Voici un aperçu de la façon d'utiliser https://en.wikipedia.org/wiki/Rubin_causal_model). L'auteur a écrit une série de articles de blog sur ses fonctionnalités.
La plupart des ensembles de données pouvant être téléchargés sont statiques, donc dans cet article, j'utiliserai ma fonction pour générer les données. Cela présente deux avantages. La capacité de générer et de générer des ensembles de données avec des propriétés spécifiques, et la capacité «d'intervenir» directement dans le système de génération de données pour vérifier si l'inférence est correcte. Tous ces générateurs de données génèrent des échantillons i.i.d. à partir de plusieurs distributions et renvoient les résultats sous forme de trames de données pandas. Les fonctions qui génèrent ces ensembles de données se trouvent dans le fichier joint datagenerators.py sur github ici.
Tout d'abord, regardons un exemple de motivation.
introduction
Un jour, un chef d'équipe a remarqué que certains membres de l'équipe portaient des casquettes cool, ce qui avait tendance à réduire la productivité de ces membres de l'équipe. Sur la base des données, le chef d'équipe est basé sur le fait que les membres de l'équipe portent des chapeaux cool ($ X = 1 $ pour les chapeaux cool, $ X = 0 $ pour les chapeaux non cool), et s'ils sont productifs ($). Démarrer l'enregistrement (Y = 1 $ est productif, $ Y = 0 $ est improductif).
Après une semaine d'observation, l'équipe a obtenu l'ensemble de données suivant:
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 |
La première question posée par le chef d'équipe est: "Une personne portant un chapeau cool est-elle plus susceptible d'être plus productive que quelqu'un qui ne le fait pas?" Cela signifie estimer les quantités suivantes:
Et vous pouvez l'estimer directement à partir des données:
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}
Les gens avec des chapeaux sympas semblent moins productifs.
Vous pouvez également exécuter un test statistique pour vous en assurer.
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
C'est une petite valeur p. Les statisticiens seront fiers
Vous pouvez utiliser ces informations pour expliquer ce que vous pensez de la probabilité que quelqu'un porte un chapeau cool. Tant que nous pensons qu'il est «dérivé de la même distribution» que les observations précédentes, on s'attend à ce que la même corrélation existe.
Le problème est lorsque les chefs d'équipe essaient d'utiliser cette information pour débattre de l'opportunité de forcer les gens à porter des chapeaux sympas. Si un chef d'équipe fait cela, il peut changer radicalement le système que nous échantillonnons, changer ou vice versa des corrélations précédemment observées.
Le moyen le plus propre de mesurer réellement l'impact de certains changements dans votre système est d'exécuter un test comparatif aléatoire (https://en.wikipedia.org/wiki/Randomized_controlled_trial). Plus précisément, choisissez au hasard qui obtient un chapeau cool et qui ne le fait pas, et voyez la différence dans la valeur $ y $ que vous obtenez. Cela élimine les effets des variables entrelacées (https://en.wikipedia.org/wiki/Confounding) qui peuvent affecter les métriques qui vous intéressent.
Maintenant que nous avons généré le jeu de données à partir d'un processus connu (dans ce cas la fonction que j'ai créée), nous pouvons directement y intervenir pour mesurer l'efficacité du test A / B.
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}
Soudain, la direction de l'effet du port d'un chapeau cool semble s'être inversée.
Que se passe-t-il?
Remarque: Dans l'exemple ci-dessus et tous les exemples suivants, l'exemple est iid et SUTVA. On suppose que vous suivez / wiki / Rubin_causal_model # Stable_unit_treatment_value_assumption_% 28SUTVA% 29). Fondamentalement, cela signifie que lorsqu'une personne choisit ou est obligée de porter un chapeau vraiment cool, cela n'affecte pas les choix ou les effets d'une autre qui porte vraiment un chapeau cool. Je vais. Structurellement, tous les générateurs de données synthétiques que j'utilise ont cette propriété. En réalité, c'est une autre chose que vous devez supposer être vraie.
Définition de la relation causale
L'exemple précédent montre l'ancienne proclamation statisticienne:
** La corrélation ne signifie pas une relation causale **
"Relation causale" est un mot vague et philosophique. Dans le contexte actuel, nous l'utilisons pour signifier "quel est l'effet de changer $ X $ sur $ Y $?"
Pour être précis, $ X $ et $ Y $ sont des variables aléatoires, et l '«effet» que vous voulez savoir est d'attribuer une valeur spécifique à $ X $. Alors comment la distribution de $ Y $ change. Cet acte de forcer une variable à avoir une valeur spécifique est appelé «intervention».
Dans l'exemple précédent, si nous n'intervenions pas dans le système, nous obtiendrions une distribution d'observation de $ Y $, à condition d'observer $ X $.
Nous intervenons lorsque nous forçons les gens à avoir l'air cool. Et la distribution de $ Y $ est donnée par la distribution intervention.
En général, les deux ne sont pas les mêmes.
La question à laquelle ces notes essaient de répondre est de savoir comment pouvons-nous déduire la distribution des interventions lorsque seules des données d'observation sont disponibles. C'est une question utile car il existe de nombreuses situations irréalistes, réalisables ou contraires à l'éthique dans lesquelles il est impossible, faisable ou contraire à l'éthique d'effectuer un test A / B pour mesurer directement l'efficacité d'une intervention. Même dans cette situation, j'aimerais pouvoir dire quelque chose sur l'effet de l'intervention.
Résultat potentiel
Une façon d'aborder ce problème consiste à introduire deux nouvelles variables aléatoires dans le système: $ Y_ {0} $ et $ Y_ {1} $, Résultats potentiels (http: // www. Connu sous le nom de stat.unipg.it/stanghellini/rubinjasa2005.pdf). Nous supposons que ces variables existent et peuvent être traitées comme toute autre variable aléatoire.
$ Y $ est défini comme suit:
Y = Y_{1} whenX=1 Y = Y_{0} whenX=0
Il s'agit d'une distribution sous-jacente avec des valeurs manquantes (https://en.wikipedia.org/wiki/Missing_data), sur la façon dont la distribution change en cas d'intervention dans le problème. Décale de à environ les données dessinées avec iid. Selon certaines hypothèses sur les raisons pour lesquelles les valeurs sont manquantes, il existe des théories bien développées sur la façon d'estimer les valeurs manquantes.
Cible
Souvent, nous ne nous soucions pas de la distribution complète des interventions, $ P (Y | \ hbox {do} (X)) $, et une estimation de la différence entre les moyennes entre les deux groupes est suffisante. Il s'agit d'un nombre connu sous le nom d'effet de traitement moyen (ATE) (https://en.wikipedia.org/wiki/Average_treatment_effect).
Lorsque vous effectuez des tests A / B et comparez les moyennes de chaque groupe, cela est directement lié aux nombres que nous mesurons.
L'estimation de cette valeur à partir de la distribution observée donne les résultats suivants:
Cela ne correspond généralement pas au vrai ATE. parce que
Les deux grandeurs liées
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"
Voilà pourquoi.
Une façon d'interpréter l'ATC est de mesurer l'effet du traitement uniquement des échantillons non traités, et vice versa pour l'ATT. Selon le cas d'utilisation, ils peuvent être une mesure plus naturelle de ce qui vous tient à cœur. Vous pouvez tous les estimer à l'aide des techniques suivantes:
Faire des hypothèses
Lorsque vous effectuez des tests A / B, randomisez l'affectation $ X $. Cela a pour effet de vous permettre de choisir de révéler les variables $ Y_ {1} $ ou $ Y_ {0} $. Cela rend le résultat indépendant de la valeur de $ X $.
Écrivez ceci comme suit:
Autrement dit, la distribution de $ X, Y_ {0}, Y_ {1} $ est factorisée comme suit:
Si cette indépendance est garantie, alors:
Si vous souhaitez estimer l'ATE à l'aide de données d'observation, vous devez utiliser d'autres informations dont vous disposez sur l'échantillon. En particulier, il faut ** supposer ** qu'il existe suffisamment d'informations supplémentaires pour expliquer pleinement les choix de traitement pour chaque sujet.
En appelant cette information supplémentaire la variable aléatoire $ Z $, cette hypothèse peut s'écrire:
Ou
Cela signifie que le traitement observé $ X $ que reçoit le spécimen est entièrement expliqué par $ Z $. Ceci est parfois appelé "Ignorer l'hypothèse" (https://en.wikipedia.org/wiki/Ignorability).
Dans cet exemple, l'exemple de motivation pour un chapeau cool est qu'il y a d'autres facteurs (appelons cela des «compétences») qui affectent à la fois la productivité humaine et le fait qu'ils portent ou non un chapeau cool. Veux dire Dans l'exemple ci-dessus, un homme du métier est plus productif et moins susceptible de porter un chapeau cool. Ces faits peuvent expliquer ensemble pourquoi les effets de Cool Hat semblaient s'inverser lors de l'exécution des tests A / B.
La division des données selon que la personne est compétente ou non révèle une relation positive entre le port de chapeaux sympas et la productivité dans chaque sous-groupe.
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}
Malheureusement, nous ne pouvons pas observer $ Y_ {0} $ et $ Y_ {1} $ dans le même échantillon, nous ne pouvons donc pas tester les hypothèses suivantes.
Il doit être évalué en utilisant la connaissance du système que nous étudions.
La qualité des prédictions que vous faites dépend de la qualité de cette hypothèse. Le paradoxe de Simpson (http://www.degeneratestate.org/posts/2017/Oct/22/generating-examples-of-simpsons-paradox/) dit si $ Z $ ne contient pas toutes les variables intriquées. , Un exemple extrême du fait que les inférences que nous faisons peuvent être erronées. Facebook publie un bon article comparant différentes approches d'inférence causale utilisant des tests directs A / B, montrant comment l'efficacité est surestimée lorsque l'indépendance conditionnelle n'est pas maintenue. Je suis.
Une fois que vous faites cette hypothèse, il existe plusieurs techniques pour l'aborder. Le reste de cet article décrira certaines des approches les plus simples, mais gardez à l'esprit qu'il s'agit d'un domaine de recherche en cours.
Modélisation des contre-faits
D'après ce qui précède, il est clair que l'ATE peut être estimé si $ Y_ {0} $ et $ Y_ {1} $ sont connus. Alors pourquoi ne pas essayer de les modéliser directement?
Plus précisément, vous pouvez faire les estimations suivantes.
\hat{Y}_{0}(Z) = E[Y|Z, X=0] \hat{Y}_{1}(Z) = E[Y|Z, X=1]
Si nous pouvons modéliser ces deux grandeurs, nous pouvons estimer l'ATE comme suit.
Le succès de cette approche dépend de la manière dont les résultats potentiels peuvent être modélisés. Pour le voir en action, utilisons le processus de génération de données suivant.
observed_data_1 = dg.generate_dataset_1()
observed_data_1.plot.scatter(x="z", y="y", c="x", cmap="rainbow", colorbar=False);
Avant de modéliser les contre-faits, examinons les données. En regardant comment $ Y $ est distribué, il semble qu'il y ait une petite différence entre les deux groupes.
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>
Cela peut être confirmé en regardant la différence moyenne entre les deux groupes.
print("Observed ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**estimate_uplift(observed_data_1)))
Observed ATE: 0.113 (0.103)
Cependant, si vous regardez la distribution de $ Z $, qui est la covariance, vous pouvez voir qu'il existe des différences entre les groupes.
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>
Ceci est un problème si vous pensez que $ Z $ a un effet sur $ Y $. Nous avons besoin d'un moyen de faire la distinction entre l'impact de $ X $ sur $ Y $ et l'impact de $ Z $ sur $ Y $.
Vous pouvez vérifier l'ATE réel avec un test A / B qui le simule et confirmer qu'il s'agit de la différence par rapport à la valeur observée.
print("Real ATE: {estimated_effect:.3f} ({standard_error:.3f})".format(**run_ab_test(dg.generate_dataset_1)))
Real ATE: -0.479 (0.026)
Mais que faire si vous ne pouvez pas exécuter ce test A / B? Vous devez modéliser votre système.
Le type de modèle le plus simple est le modèle linéaire. Plus précisément, supposez ce qui suit:
Si cela est exact, laissez le modèle s'adapter aux données
Vous pouvez utiliser la régression linéaire pour obtenir une estimation de l'ATE.
Le paquet causalinference fournit une interface simple pour ce faire.
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
«causalinference» renvoie une estimation de l'ATE et les propriétés statistiques de cette estimation. L'intervalle de confiance indiqué pour l'estimation est l'intervalle de confiance en supposant que le modèle décrit avec précision le contrefact, et l'intervalle de confiance quant à la qualité de la description du contrefact par le modèle. Il est important de reconnaître qu'une telle chose n'existe pas.
Dans ce cas, le package a réussi à identifier le bon ATE - ce qui est bon, mais le processus de génération de données est spécialement conçu pour répondre aux hypothèses. Regardons quelques cas qui peuvent échouer.
Le premier cas est celui où l'effet n'est pas simplement additif.
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
Cela peut généralement être surmonté en utilisant des estimations plus solides. Une approche simple et non paramétrique est la méthode Matching (https://en.wikipedia.org/wiki/Matching_%28statistics%29). L'idée est de trouver des échantillons similaires non traités pour chaque échantillon traité et de comparer ces valeurs directement. Ce que vous entendez par «similaire» dépend de votre cas d'utilisation particulier.
Le package «causalinference» met en œuvre la mise en correspondance pour chaque unité en sélectionnant l'unité la plus similaire parmi les autres groupes de traitement par remplacement et en utilisant la différence entre ces deux unités pour calculer l'ATE. .. Par défaut, les choix d'appariement sont effectués à l'aide de la méthode du plus proche voisin dans l'espace des covariables $ Z $, et les distances sont pondérées par la variance inverse de chaque dimension.
Vous avez la possibilité de modifier le nombre d'unités à comparer et la pondération de chaque dimension dans la correspondance. Consultez la documentation (http://laurence-wong.com/software/matching) pour plus d'informations.
L'estimation correspondante peut être calculée avec le code suivant.
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
L'intervalle de confiance de l'estimation contient le véritable ATE.
Déséquilibre covariable
Un problème plus difficile à résoudre est lorsque les covariables utilisées sont déséquilibrées: l'espace des covariables a une zone qui ne contient que des échantillons traités ou non. Ici, nous devons extrapoler l'effet du traitement - cela dépend fortement du modèle hypothétique que nous utilisons.
L'exemple ci-dessous le démontre.
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
L'estimation MCO ne saisit pas l'effet réel, et l'estimation d'appariement s'améliore un peu, mais il n'y a pas assez d'informations dans les données pour extrapoler complètement à une zone sans chevauchement.
Cet exemple peut sembler incohérent, mais une fois que vous commencez à examiner les covariables de dimension supérieure, le problème devient plus courant.
causalinference fournit un outil pratique pour évaluer rapidement les variables en double en utilisant la propriété 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
Ici, la différence normalisée est définie comme:
Il ne s'agit pas d'un test statistique rigoureux, mais il fournit des indicateurs du chevauchement entre chaque covariable. Les valeurs supérieures à 1 suggèrent qu'il n'y a pas beaucoup de duplication.
Score de propension
Les scores de propension sont des estimations de la probabilité qu'un sujet soit traité compte tenu d'une covariable.
Vous pouvez estimer cela comme vous le souhaitez, mais une fois que vous l'avez estimé, vous pouvez faire plusieurs choses.
Pondération par score de tendance inverse
Le problème avec la mesure de l'inférence causale est que je veux connaître la quantité $ E [Y_ {i}] $, mais gardez à l'esprit que je ne peux accéder qu'à un échantillon de $ E [Y_ {i} | X = i] $. S'il te plait donne moi.
La probabilité de résultat potentiel peut être étendue comme suit:
C'est vrai
Suggestions qui peuvent être estimées.
Par conséquent, vous pouvez récupérer des résultats potentiels en pondérant chaque point avec un score de propension inverse. Le résultat est la méthode de pondération des scores de tendance inversée (https://en.wikipedia.org/wiki/Inverse_probability_weighting).
Voyons quoi faire avec l'un des ensembles de données précédents.
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)
Les méthodes du package CausalInference ʻest_propensity_s et ʻest_propensity peuvent être utilisées pour estimer les tendances en utilisant la régression logistique sur les covariables.
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
Le générateur de données utilisé dans cet exemple vous permet de bien décrire la relation avec une simple régression logistique. Dans le générateur de données utilisé dans cet exemple, la relation peut être bien décrite par une simple régression logistique, mais si nous utilisons la fonction de régression logistique de «sklean» (la normalisation est utilisée par défaut). Si vous essayez d'estimer le score de propension, vous obtiendrez la mauvaise réponse.
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
Mieux que nos estimations naïves, mais pas correctes.
Estimateur pondéré robuste
De cette façon, vous pouvez combiner une estimation pondérée du score de propension inverse avec une estimation linéaire de l'effet pour essayer de réduire les défauts dans l'un ou l'autre modèle. Cela se fait en pré-exécutant une régression linéaire pondérée sur les données, chaque point étant pondéré par un score de propension inverse. Le résultat est estimateur pondéré doublement robuste.
L'idée est qu'il y a un biais dans lequel les échantillons sont traités dans les données d'observation, de sorte que les échantillons qui ont été traités mais moins susceptibles d'être traités sont plus importants et ont plus de poids. Devrait être donné.
Il existe les méthodes suivantes pour appliquer cela.
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}}
Score de non-confusion et de propension
Dans la section précédente, nous avons supposé que le résultat et le traitement étaient indépendants parce que les covariables étaient données.
Vous pouvez également faire une hypothèse légèrement plus forte. Autrement dit, le résultat est indépendant du traitement et est conditionné par la probabilité de propension.
Avec cette hypothèse, vous pouvez réduire les dimensions de la variable d'intrication. Cela vous permet d'exécuter certaines techniques qui peuvent ne pas fonctionner dans des paramètres dimensionnels plus élevés.
garniture
Nous avons précédemment confirmé que les déséquilibres de covariables peuvent causer des problèmes. Une solution simple consiste à faire une prédiction contrefactuelle uniquement dans les zones avec un bon chevauchement, ou à "rogner" les points où il n'y a pas de bon chevauchement. Pour les données de grande dimension, il peut être difficile de définir un «bon chevauchement» - l'utilisation de scores de propension pour définir un chevauchement est une façon de résoudre ce problème.
Regardons un ensemble de données avec moins de chevauchement.
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 fournit un moyen de découper les données en fonction des scores de propension.
# 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
Vous pouvez voir le reste des données comme suit:
# 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)
Cela fonctionne plutôt bien dans de tels cas.
Lors de l'application de l'ajustement, il est explicitement indiqué que l'inférence causale ne peut être faite que pour certains échantillons de l'espace des covariables. Rien ne peut être dit sur l'ATE pour les spécimens hors de ces régions.
Stratification
Les estimations stratifiées ou bloquées constituent une autre utilisation des scores de propension. Il consiste à regrouper les points de données en groupes de tendances similaires et à estimer l'ATE au sein de ces groupes. Encore une fois, CausalInference fournit une excellente interface pour ce faire.
Utilisez la méthode stratify (pour les limites de stata définies par l'utilisateur) ou la méthode stratify_s (sélection automatique des limites) pour déterminer la couche.
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
Utilisez ensuite la méthode ʻest_via_blocking` pour combiner les estimations de ces couches en un ATE global.
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
Cela fonctionne également bien.
La stratification des données en groupes par score de propension est utile si vous n'avez aucune connaissance préalable de ce qui constitue une unité «similaire», mais ce n'est pas la seule façon. Si vous savez à l'avance que différents groupes d'échantillons sont susceptibles d'être affectés par les interventions de la même manière, divisez l'échantillon en ces groupes avant d'estimer l'ATE et regroupez les résultats au niveau mondial. Il est logique d'obtenir un ATE.
Quelle méthode utiliser?
Cela couvre la plupart des méthodes courantes d'inférence causale à partir des données observées. La question restante est "comment et comment décider de la méthode à utiliser?" Ce n'est pas une question simple. Il existe également des techniques automatisées comme cet article, mais je ne les ai pas essayées.
En fin de compte, afin de choisir votre technique, vous devez faire des hypothèses sur la façon dont vous construisez le contre-fait. L'appariement est une bonne approche si vous faites confiance aux données pour avoir un bon chevauchement dans l'espace des covariables. Si ce n'est pas le cas, vous devez utiliser un modèle fiable pour réussir à extrapoler à la région inexplorée, ou faire l'hypothèse que quelque chose comme un score de propension capture suffisamment d'informations pour présumer de la négligence. Il y a.
Pour souligner que toutes ces méthodes peuvent échouer, voici un autre exemple. Contrairement à l'exemple précédent, il existe une ou plusieurs covariables. Comme tous les générateurs de données précédents, ce générateur de données suit également les hypothèses suivantes par conception:
Essayons aveuglément les méthodes dont nous avons discuté jusqu'à présent et voyons ce qui se passe.
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);
La ligne pointillée à côté montre le véritable ATE de cet ensemble de données.
Non seulement les résultats de chaque méthode sont différents, mais toutes les méthodes manquent la vraie valeur.
Cela devrait être un avertissement sur les limites de ce type de technique. Il peut être très intéressant pour le lecteur d'essayer quelles propriétés de l'ensemble de données font que ces techniques s'écartent de leurs vraies valeurs.
Structure du raisonnement causal
J'espère qu'à ce stade, vous serez conscient de l'importance de l'hypothèse de négligence.
Ce dont je n'ai pas encore parlé, c'est de savoir comment choisir $ Z $ pour que cela soit vrai. En fin de compte, il doit provenir de la connaissance du domaine du système étudié. Il existe un ensemble d'outils puissants appelés le modèle de relation causale (https://en.wikipedia.org/wiki/Causal_graph). Il peut être utilisé pour coder les connaissances sur le système à l'étude en tant que modèle graphique du système, ou pour déduire les hypothèses d'indépendance conditionnelle comme décrit ci-dessus.
Le prochain article sur le raisonnement causal en discutera.
Une autre question soulevée par cet article est de savoir si la seule façon de faire des inférences causales est d'ajuster les variables d'intrication. Pas si - vous pouvez l'utiliser ailleurs dans Posts later Jetons un coup d'œil à la technique.
code
Les notes pour cet article peuvent être trouvées sur github ici.
Recommended Posts