J'ai écrit "Introduction à la vérification des effets" en Python

TL;DR

• Le code source R du livre "Introduction à la vérification de l'efficacité / Raisonnement causal pour une comparaison correcte / Bases de l'économie quantitative" a été (presque) reproduit en Python.
• https://github.com/nekoumei/cibook-python
• Dans cet article, R est principalement une bibliothèque! Je vais expliquer la partie qui ne l'est pas en Python

Présentation des livres

https://www.amazon.co.jp/dp/B0834JN23Y

Le tableau ci-dessus est répertorié sur Amazon, donc je pense qu'il est rapide de le voir. .. C'est un très bon livre. Comment se débarrasser des préjugés pour prendre des décisions précises? Différentes méthodes d'inférence causale (score de propension / DiD / RDD, etc.) sont introduites avec une implémentation par code source R. Tout au long, comment pouvons-nous utiliser le raisonnement causal pour vérifier l'efficacité de problèmes réels? Il a été rédigé du point de vue de, et j'ai trouvé que c'était très pratique.

Écrit en Python

https://github.com/nekoumei/cibook-python J'ai créé un Jupyter Notebook qui correspond au code source public R d'origine (https://github.com/ghmagazine/cibook). De plus, lors de l'implémentation de ce Python, la bibliothèque de visualisation graphique utilise plotly.express, sauf pour certaines parties. Les graphiques tracés ne peuvent pas être affichés dans le rendu Github Notebook, donc si vous souhaitez vérifier en ligne, veuillez vérifier les pages Github répertoriées dans le README. (Affichage du html sous les images) En ce qui concerne plotly.express, l'article suivant sera très utile en japonais. Résumé de la méthode de dessin de base de Plotly Express, le standard de facto de la bibliothèque de dessin Python à l'ère Reiwa

Points clés lors de l'écriture en Python

analyse de régression

J'utilise les modèles OLS et WLS de statistiques.

(Extrait de ch2_regression.ipynb)

##Régression multiple sur données biaisées
y = biased_df.spend
#Dans R lm, les variables catégorielles sont automatiquement converties en variables factices, donc reproduisez-les.
X = pd.get_dummies(biased_df[['treatment', 'recency', 'channel', 'history']], columns=['channel'], drop_first=True)
X = sm.add_constant(X)
results = sm.OLS(y, X).fit()
nonrct_mreg_coef = results.summary().tables[1]
nonrct_mreg_coef

Ici, il existe deux différences par rapport à la fonction R lm. (1) Il semble que lm convertit automatiquement les variables catégorielles en variables fictives. Puisque sm.OLS ne peut pas le faire automatiquement, pd.get_dummies () est utilisé. (2) Dans lm, le terme de biais est ajouté automatiquement. Ceci est également ajouté manuellement. (Référence: Statistiques: essayez une analyse de régression multiple avec Python et R)

Lire l'ensemble de données au format RData

Certains des ensembles de données utilisés dans l'implémentation sont publiés sous forme de packages R tels que experimentdatar, ou sont publiés au format RData. il y a. ~~ Il n'y a aucun moyen de les lire en Python uniquement. ~~ (S'il vous plaît laissez-moi savoir si vous en avez un) Dans la section des commentaires, upura-san m'a appris! https://qiita.com/nekoumei/items/648726e89d05cba6f432#comment-0ea9751e3f01b27b0adb Le fichier .rda est un package appelé rdata qui peut être lu sans R. (Extrait de ch2_voucher.ipynb)

parsed = rdata.parser.parse_file('../data/vouchers.rda')
converted = rdata.conversion.convert(parsed)
vouchers = converted['vouchers']

Lors de la lecture de données en utilisant R via rpy2 et de la conversion en pandas DataFrame, c'est comme suit. (Extrait de ch2_voucher.ipynb)

from rpy2.robjects import r, pandas2ri
from rpy2.robjects.packages import importr
pandas2ri.activate()
experimentdatar = importr('experimentdatar')
vouchers = r['vouchers']

Comme décrit dans ch2_voucher.ipynb, le package R est installé à l'avance dans l'environnement interactif R. De plus, l'ensemble de données utilisé par ch3_lalonde.ipynb est un fichier .dta. Cela peut être lu par pandas read_stata (). (Extrait de ch3_lalonde.ipynb)

cps1_data = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls.dta')
cps3_data = pd.read_stata('https://users.nber.org/~rdehejia/data/cps_controls3.dta')
nswdw_data = pd.read_stata('https://users.nber.org/~rdehejia/data/nsw_dw.dta')

Appariement du score de propension (appariement le plus proche)

Il existe plusieurs méthodes d'appariement pour l'appariement des scores de propension, mais le livre semble utiliser l'appariement le plus proche de MatchIt. Il ne semblait pas y avoir de bonne bibliothèque pour Python, donc je l'implémente honnêtement. (Extrait de ch3_pscore.ipynb)

def get_matched_dfs_using_propensity_score(X, y, random_state=0):
    #Calculer le score de propension
    ps_model = LogisticRegression(solver='lbfgs', random_state=random_state).fit(X, y)
    ps_score = ps_model.predict_proba(X)[:, 1]
    all_df = pd.DataFrame({'treatment': y, 'ps_score': ps_score})
    treatments = all_df.treatment.unique()
    if len(treatments) != 2:
        print('Seuls 2 groupes peuvent être mis en correspondance. 2 groupes doivent être[0, 1]Veuillez exprimer avec.')
        raise ValueError
    # treatment ==1 au groupe1, treatment ==Soit 0 group2. Puisque le groupe2 qui correspond au groupe1 est extrait, il doit s'agir d'une estimation de ATT.
    group1_df = all_df[all_df.treatment==1].copy()
    group1_indices = group1_df.index
    group1_df = group1_df.reset_index(drop=True)
    group2_df = all_df[all_df.treatment==0].copy()
    group2_indices = group2_df.index
    group2_df = group2_df.reset_index(drop=True)

    #Écart type du score global de propension* 0.2 est le seuil
    threshold = all_df.ps_score.std() * 0.2

    matched_group1_dfs = []
    matched_group2_dfs = []
    _group1_df = group1_df.copy()
    _group2_df = group2_df.copy()

    while True:
        #Trouver et faire correspondre un voisin le plus proche dans les voisins les plus proches
        neigh = NearestNeighbors(n_neighbors=1)
        neigh.fit(_group1_df.ps_score.values.reshape(-1, 1))
        distances, indices = neigh.kneighbors(_group2_df.ps_score.values.reshape(-1, 1))
        #Supprimer les doublons
        distance_df = pd.DataFrame({'distance': distances.reshape(-1), 'indices': indices.reshape(-1)})
        distance_df.index = _group2_df.index
        distance_df = distance_df.drop_duplicates(subset='indices')
        #Supprimer les enregistrements qui dépassent le seuil
        distance_df = distance_df[distance_df.distance < threshold]
        if len(distance_df) == 0:
            break
        #Extraire et supprimer les enregistrements correspondants
        group1_matched_indices = _group1_df.iloc[distance_df['indices']].index.tolist()
        group2_matched_indices = distance_df.index
        matched_group1_dfs.append(_group1_df.loc[group1_matched_indices])
        matched_group2_dfs.append(_group2_df.loc[group2_matched_indices])
        _group1_df = _group1_df.drop(group1_matched_indices)
        _group2_df = _group2_df.drop(group2_matched_indices)

    #Renvoie un enregistrement correspondant
    group1_df.index = group1_indices
    group2_df.index = group2_indices
    matched_df = pd.concat([
        group1_df.iloc[pd.concat(matched_group1_dfs).index],
        group2_df.iloc[pd.concat(matched_group2_dfs).index]
    ]).sort_index()
    matched_indices = matched_df.index

    return X.loc[matched_indices], y.loc[matched_indices]

L'appariement de 1 point dans le groupe de traitement avec 1 point dans le groupe témoin avec le score de propension le plus proche est répété. A ce moment, un seuil est établi de sorte que seules les paires dont la distance est inférieure à 0,2 fois std soient extraites. Pour plus de détails à ce sujet, voir Concept de score de propension et sa pratique. Je compare avec le résultat correspondant par MatchIt dans ch3_pscore.ipynb, mais il n'y a pas de correspondance exacte. La conclusion est la même et l'équilibre des covariables est bon (voir ci-dessous), il est donc généralement bon. Encore une fois, il n'y a pas de bibliothèque pratique comme le cobalt love.plot () de R, donc je la visualise moi-même. (Depuis images / ch3_plot1.html)

Estimation pondérée par probabilité inverse (IPW)

L'avantage d'IPW est qu'il est simple à mettre en œuvre. Ceci est également comparé à WeightIt, mais cela semble généralement correct. (Extrait de ch3_pscore.ipynb)

def get_ipw(X, y, random_state=0):
    #Calculer le score de propension
    ps_model = LogisticRegression(solver='lbfgs', random_state=random_state).fit(X, y)
    ps_score = ps_model.predict_proba(X)[:, 1]
    all_df = pd.DataFrame({'treatment': y, 'ps_score': ps_score})
    treatments = all_df.treatment.unique()
    if len(treatments) != 2:
        print('Seuls 2 groupes peuvent être mis en correspondance. 2 groupes doivent être[0, 1]Veuillez exprimer avec.')
        raise ValueError
    # treatment ==1 au groupe1, treatment ==Soit 0 group2.
    group1_df = all_df[all_df.treatment==1].copy()
    group2_df = all_df[all_df.treatment==0].copy()
    group1_df['weight'] = 1 / group1_df.ps_score
    group2_df['weight'] = 1 / (1 - group2_df.ps_score)
    weights = pd.concat([group1_df, group2_df]).sort_index()['weight'].values
    return weights

CausalImpact Comme décrit dans ch4_did.ipynb, les deux bibliothèques suivantes sont utilisées à des fins de comparaison.

• dafiti/causalinpact: https://github.com/dafiti/causalimpact
• tcassou/causal_impact: https://github.com/tcassou/causal_impact

Comme le résultat est beau, j'utilise généralement l'impact causal de Dafiti, mais c'est l'impact causal de tcassou que le résultat de l'estimation était proche du livre (impact causal original de R). Les deux semblent être estimés par le modèle d'espace d'état des modèles de statistiques, mais je n'ai pas compris la différence de mise en œuvre. S'il vous plaît, faites-moi savoir. .. Dans les deux cas, l'erreur d'estimation est extrêmement faible par rapport à l'implémentation R. Ce n'est peut-être pas si bon qu'il existe de nombreuses covariables.

dafiti / diagramme d'impact causal

tracé de tcassou / causal_impact

Conception discontinue de régression (RDD)

Avec R, vous pouvez l'exécuter rapidement avec rddtools, mais avec Python ce n'est pas le cas. Tout d'abord, vérifiez l'équation de régression utilisée dans rdd_reg_lm de rddtools. (Référence: https://cran.r-project.org/web/packages/rddtools/rddtools.pdf P23)

Y = α + τD + β_1(X-c)+ β_2D(X-c) + ε

En passant, je pensais que RDD créerait des modèles de régression à gauche et à droite de la valeur de coupure, et prendrait la différence entre les valeurs estimées à la valeur de coupure, mais une régression. Il s'exprime par une expression. Est-ce la même chose en termes de sens? Où D est une variable binaire avec une valeur de 1 si X est supérieur ou égal à la valeur seuil et 0 si elle est antérieure. Le coef de D est la valeur que vous voulez vérifier comme quantité d'effet. De plus, c est la valeur limite. Sur la base de ce qui précède, nous allons le mettre en œuvre. Je me réfère également au code source de rddtools pour la mise en œuvre. (https://github.com/MatthieuStigler/RDDtools/blob/master/RDDtools/R/model.matrix.RDD.R) (Extrait de ch5_rdd.ipynb)

class RDDRegression:
#Paquet R rddtools rdd_reg_Reproduire lm
#Référence: https://cran.r-project.org/web/packages/rddtools/rddtools.pdf P23
    def __init__(self, cut_point, degree=4):
        self.cut_point = cut_point
        self.degree = degree
        
    def _preprocess(self, X):
        X = X - threshold_value
        X_poly = PolynomialFeatures(degree=self.degree, include_bias=False).fit_transform(X)
        D_df = X.applymap(lambda x: 1 if x >= 0 else 0)
        X = pd.DataFrame(X_poly, columns=[f'X^{i+1}' for i in range(X_poly.shape[1])])
        X['D'] = D_df
        for i in range(X_poly.shape[1]):
            X[f'D_X^{i+1}'] = X_poly[:, i] * X['D']
        return X
    
    def fit(self, X, y):
        X = X.copy()
        X = self._preprocess(X)
        self.X = X
        self.y = y
        X = sm.add_constant(X)
        self.model = sm.OLS(y, X)
        self.results = self.model.fit()
        coef = self.results.summary().tables[1]
        self.coef = pd.read_html(coef.as_html(), header=0, index_col=0)[0]
        
    def predict(self, X):
        X = self._preprocess(X)
        X = sm.add_constant(X)
        return self.model.predict(self.results.params, X)

Les caractéristiques polynomiales de scicit-learn sont utilisées comme prétraitement pour effectuer une régression polynomiale. (De images / ch5_plot2_3.html) J'ai pu bien l'estimer. Comme vous pouvez le voir en regardant le Notebook, l'effet estimé est presque le même que celui du livre. C'était bien.

nonparametric RDD (RDestimate) Je ne pouvais pas ... C'est presque la partie de "(presque) reproduit en Python" au début. RDestimate utilise la méthode d'Immunos et Kalyanaraman (2012) Optimal Bandwidth Choice for the Regression Discontinuity Estimator pour sélectionner la bande passante optimale, mais je ne savais pas trop comment estimer cette bande passante. Dans ch5_rdd.ipynb, j'ai essayé d'implémenter MSE pour trouver la bande passante optimale lors de la modification approximative de la bande passante, mais cela ne semble pas très bien fonctionner. (De ch5_plot4.html) Hmmm, n'est-ce pas? Il semble qu'il est impossible de prédire en dehors de la bande passante si elle est estimée en réduisant uniquement la bande passante, mais comment le faites-vous réellement? Ou peut-être ne suis-je intéressé que par l'estimation du quartier limite pour ne pas trop m'inquiéter.

en conclusion

Si vous trouvez quelque chose qui ne va pas dans votre compréhension ou votre mise en œuvre, ou si quelque chose ne va pas, veuillez nous en informer. Twitter