Introduction à la vérification de l'efficacité Chapitre 3 écrit en Python

introduction

Introduction à la vérification de l'efficacité - Raisonnement causal pour une comparaison correcte / Bases de l'économie quantitative Reproduire le code source en Python Faire.

J'ai déjà un Exemple d'implémentation de grand ancêtre, mais je le laisserai comme mémo pour mon étude.

Cet article couvre le chapitre 3. Le code est également posté sur github. De plus, les noms de variables et le contenu du traitement sont essentiellement implémentés dans le livre.

Retour logistique

Vous pouvez implémenter une régression logistique avec des modèles scicit-learn ou stats.

scikit-learn

scikit-learn doit être factice.

Retour logistique de sklearn

from sklearn.linear_model import LogisticRegression

X = pd.get_dummies(biased_data[['recency', 'history', 'channel']])  #Faire du canal une variable factice
treatment = biased_data['treatment'] 

model = LogisticRegression().fit(X, treatment)

statsmodels

Utilisez les classes glm ou logit.

statsmodels n'a pas besoin d'être factice, mais il semble que certaines valeurs de catégorie ne soient pas utilisées dans le modèle. (Channel = Multichannel n'est pas utilisé dans l'exemple ci-dessous. La cause n'est pas claire.) Si vous êtes inquiet, vous pouvez en faire un mannequin et il correspondra au résultat de sklearn.

Analyse de régression des modèles de statistiques

from statsmodels.formula.api import glm, logit

model = glm('treatment ~ history + recency + channel', data=biased_data, family=sm.families.Binomial()).fit()
# model = logit('treatment ~ history + recency + channel', data=biased_data).fit()  #Même résultat que ci-dessus

Estimation à l'aide du score de propension

«L'appariement des scores de propension» et «IPW», qui sont mentionnés comme méthodes d'estimation utilisant des scores de propension, peuvent être mis en œuvre sur DoWhy.

Avant d'effectuer chaque analyse, définissez une instance pour l'inférence causale commune. À ce moment-là, notez que la variable d'intervention doit être booléenne.

Préparation

from dowhy import CausalModel

biased_data = biased_data.astype({"treatment":'bool'}, copy=False)  #traitement booléen

model=CausalModel(
    data = biased_data,
    treatment='treatment',
    outcome='spend',
    common_causes=['recency', 'history', 'channel']
)

Appariement du score de propension

Correspondance la plus proche

L'appariement le plus proche est une méthode d'appariement 1: 1 entre les échantillons avec intervention et les échantillons sans intervention, qui sont expliqués dans le livre, avec des scores de propension similaires. Dans ce livre, seul ATT est estimé, mais DoWhy peut également estimer ATE. Notez que la valeur par défaut est ATE. Puisque les deux sont calculés par l'implémentation interne, ce serait bien si vous pouviez obtenir les deux.

De plus, cet ATE semble être calculé à partir de l'ATT et de l'ATC (valeur attendue de l'effet d'intervention sur l'échantillon de non-intervention).

Correspondance la plus proche

nearest_ate = model.estimate_effect(
    identified_estimand, 
    method_name="backdoor.propensity_score_matching",
    target_units='ate',
)
nearest_att = model.estimate_effect(
    identified_estimand, 
    method_name="backdoor.propensity_score_matching",
    target_units='att', 
)

print(f'ATE: {nearest_ate.value}')
print(f'ATT: {nearest_att.value}')

Appariement stratifié

Bien que cela ne soit pas mentionné dans ce livre, l'appariement stratifié est également mis en œuvre.

En regardant le code source, il semble que les données soient divisées dans la couche num_strata par le même nombre et estimées. clipping_threshold compare le nombre de cas avec et sans intervention dans chaque couche, et exclut les couches inférieures à cette valeur. (Probablement pour exclure la couche avec extrêmement peu de données d'intervention)

Appariement stratifié

stratification_ate = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.propensity_score_stratification",
    target_units='ate',
    method_params={'num_strata':50, 'clipping_threshold':5},
)
stratification_att = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.propensity_score_stratification",
    target_units='att',
    method_params={'num_strata':50, 'clipping_threshold':5},
)

print(f'ATE: {stratification_ate.value}')
print(f'ATT: {stratification_att.value}')

IPW

IPW

ipw_estimate = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.propensity_score_weighting",
    target_units='ate',
)

print(f'ATE: {ipw_estimate.value}')

(Bonus) Analyse de régression

L'analyse de régression effectuée au chapitre 2 est également possible.

python

estimate = model.estimate_effect(
    identified_estimand,
    method_name="backdoor.linear_regression", #test_significance=True
)

print('Causal Estimate is ' + str(estimate.value))

Visualisation de la différence moyenne standardisée

La visualisation des différences moyennes standardisées pour confirmer l'équilibre des covariables n'est pas prise en charge. Par conséquent, c'est gênant, mais il est calculé par vous-même.

Appariement du score de propension

Il semble que vous ne puissiez pas vous référer à la correspondance de l'appariement par score de propension, alors mettez en œuvre l'appariement vous-même. L'implémentation est basée sur Around here.

Je ne sais pas comment trouver la «distance» et la calculer correctement, alors ne l'utilisez pas comme référence.

Différence moyenne standard d'appariement du score de propension

from sklearn.neighbors import NearestNeighbors

#Calcul de l'ASAM
def calculate_asam(data, treatment_column, columns):
  treated_index = data[treatment_column]
  data = pd.get_dummies(data[columns])

  asam = data.apply(lambda c: abs(c[treated_index].mean() - c[~treated_index].mean()) / c.std())
  asam['distance'] = np.sqrt(np.sum(asam**2))  #Je ne comprends pas la définition
  return asam

#Appariement du score de propension
def get_matching_data(data, treatment_column, propensity_score):
  data['ps'] = propensity_score
  treated = data[data[treatment_column] == 1]
  control = data[data[treatment_column] == 0]

  #Faire correspondre les échantillons avec l'intervention
  control_neighbors = NearestNeighbors(n_neighbors=1, algorithm='ball_tree').fit(control['ps'].values.reshape(-1, 1))
  distances, indices = control_neighbors.kneighbors(treated['ps'].values.reshape(-1, 1))
  #Vous pouvez définir le seuil de correspondance avec les distances(Pas compatible)

  matching_data = pd.concat([treated, control.iloc[indices.flatten()]])
  return matching_data

matching_data = get_matching_data(biased_data[['treatment', 'recency', 'channel', 'history']], 'treatment', nearest_ate.propensity_scores)
unadjusted_asam = calculate_asam(biased_data, 'treatment', ['recency', 'history', 'channel'])
matching_adjusted_asam = calculate_asam(matching_data, 'treatment', ['recency', 'history', 'channel'])

IPW

Avec DescrStatsW, vous pouvez facilement calculer des statistiques pondérées. Cependant, la valeur non pondérée est-elle légèrement différente du résultat de numpy? Soyez donc prudent lorsque vous le manipulez.

Différence moyenne normalisée IPW

from statsmodels.stats.weightstats import DescrStatsW

def calculate_asam_weighted(data, treatment_column, columns, propensity_score):
  data = pd.get_dummies(data[[treatment_column] + columns])
  data['ps'] = propensity_score
  data['weight'] = data[[treatment_column, 'ps']].apply(lambda x: 1 / x['ps'] if x[treatment_column] else 1 / (1 - x['ps']), axis=1)

  asam_dict = dict()
  for column_name, column_value in data.drop([treatment_column, 'ps', 'weight'], axis=1).iteritems():
    treated_stats = DescrStatsW(column_value[data[treatment_column]], weights=data[data[treatment_column]]['weight'])
    control_stats = DescrStatsW(column_value[~data[treatment_column]], weights=data[~data[treatment_column]]['weight'])
    asam_dict[column_name] = abs(treated_stats.mean - control_stats.mean) / treated_stats.std

  asam = pd.Series(asam_dict)
  asam['distance'] = np.sqrt(np.sum(asam**2))  #Je ne comprends pas la définition
  return asam

ipw_adjusted_asam = calculate_asam_weighted(biased_data, 'treatment', ['recency', 'history', 'channel'], ipw_estimate.propensity_scores)

Visualisation

Visualisation

%matplotlib inline
import matplotlib.pyplot as plt

plt.scatter(unadjusted_asam, unadjusted_asam.index, label='unadjusted')
plt.scatter(matching_adjusted_asam, matching_adjusted_asam.index, label='adjusted(matching)')
plt.scatter(ipw_adjusted_asam, ipw_adjusted_asam.index, label='adjusted(ipw)')
plt.vlines(0.1, 0, len(unadjusted_asam)-1, linestyles='dashed')
plt.legend()
plt.show()

En quelque sorte, l'IPW est meilleur, mais probablement parce que l'appariement par score de propension permet la duplication des échantillons correspondants dans l'appariement du plus proche voisin et est biaisé en faveur de certains échantillons.

Relation

• J'ai écrit "Introduction à la vérification des effets" en Python
• DoWhy | Making causal inference easy — DoWhy | Making Causal Inference Easy documentation