Geschrieben "Einführung in die Effektüberprüfung" in Python

TL;DR

• Der R-Quellcode des Buches "Einführung in die Überprüfung der Wirksamkeit / kausale Argumentation für einen korrekten Vergleich / Grundlagen der quantitativen Ökonomie" wurde (fast) in Python reproduziert.
• https://github.com/nekoumei/cibook-python
• In diesem Artikel ist hauptsächlich R eine Bibliothek! Ich werde den Teil erklären, der in Python nicht so ist

Einführung von Büchern

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

Die obige Tabelle ist bei Amazon gelistet, daher ist es meiner Meinung nach schnell zu sehen. .. Es ist ein sehr gutes Buch. Wie kann man Vorurteile loswerden, um genaue Entscheidungen zu treffen? Verschiedene kausale Inferenzmethoden (Propensity Score / DiD / RDD usw.) werden mit einer Implementierung durch R-Quellcode eingeführt. Wie können wir kausale Argumente verwenden, um die Wirksamkeit realer Probleme zu überprüfen? Es wurde aus der Perspektive von geschrieben, und ich fand es sehr praktisch.

Geschrieben in Python

https://github.com/nekoumei/cibook-python Ich habe ein Jupyter-Notizbuch erstellt, das dem ursprünglichen öffentlichen R-Quellcode entspricht (https://github.com/ghmagazine/cibook). Bei der Implementierung dieses Python wird außerdem mit Ausnahme einiger Teile plotly.express für die Grafikvisualisierungsbibliothek verwendet. Plotly-Diagramme können im Github Notebook-Rendering nicht angezeigt werden. Wenn Sie also online prüfen möchten, überprüfen Sie bitte die in der README-Liste aufgeführten Github-Seiten. (HTML unter Bildern anzeigen) In Bezug auf plotly.express ist der folgende Artikel auf Japanisch sehr hilfreich. Zusammenfassung der grundlegenden Zeichenmethode von Plotly Express, dem De-facto-Standard der Python-Zeichenbibliothek in der Reiwa-Ära

Wichtige Punkte beim Schreiben in Python

Regressionsanalyse

Ich verwende OLS und WLS von Statistikmodellen.

(Auszug aus ch2_regression.ipynb)

##Multiple Regression auf voreingenommene Daten
y = biased_df.spend
#In R lm werden kategoriale Variablen automatisch in Dummy-Variablen konvertiert, also reproduzieren Sie sie.
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

Hier gibt es zwei Unterschiede zur R lm-Funktion. (1) Es scheint, dass lm kategoriale Variablen automatisch in Dummy-Variablen konvertiert. Da sm.OLS dies nicht automatisch tun kann, wird "pd.get_dummies ()" verwendet. (2) In lm wird der Bias-Term automatisch addiert. Dies wird auch manuell hinzugefügt. (Referenz: Statistik: Versuchen Sie eine multiple Regressionsanalyse mit Python und R)

Lesen Sie den Datensatz im RData-Format

Einige der in der Implementierung verwendeten Datasets werden als R-Pakete wie experimentdatar oder im RData-Format veröffentlicht. es gibt. ~~ Es gibt keine Möglichkeit, diese nur in Python zu lesen. ~~ (Bitte lassen Sie mich wissen, wenn Sie eine haben) Im Kommentarbereich hat mich Upura-San unterrichtet! https://qiita.com/nekoumei/items/648726e89d05cba6f432#comment-0ea9751e3f01b27b0adb Die .rda-Datei ist ein Paket namens rdata, das ohne R gelesen werden kann. (Auszug aus ch2_voucher.ipynb)

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

Beim Lesen von Daten mit R über rpy2 und beim Konvertieren in Pandas DataFrame ist dies wie folgt. (Auszug aus ch2_voucher.ipynb)

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

Wie in ch2_voucher.ipynb beschrieben, wird das R-Paket im Voraus in der interaktiven R-Umgebung installiert. Außerdem ist der von ch3_lalonde.ipynb verwendete Datensatz eine .dta-Datei. Dies kann von pandas read_stata () gelesen werden. (Auszug aus 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')

Propensity Score Matching (nächster Matching)

Es gibt verschiedene Übereinstimmungsmethoden für die Neigungsbewertung, aber das Buch scheint die engste Übereinstimmung von MatchIt zu verwenden. Es schien keine gute Bibliothek für Python zu geben, also implementiere ich sie ehrlich. (Auszug aus ch3_pscore.ipynb)

def get_matched_dfs_using_propensity_score(X, y, random_state=0):
    #Berechnen Sie die Neigungsbewertung
    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('Es können nur 2 Gruppen abgeglichen werden. 2 Gruppen müssen sein[0, 1]Bitte mit ausdrücken.')
        raise ValueError
    # treatment ==1 bis Gruppe1, treatment ==Sei 0 Gruppe2. Da Gruppe2, die mit Gruppe1 übereinstimmt, extrahiert wird, sollte es sich um eine Schätzung der ATT handeln.
    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)

    #Standardabweichung der Gesamtneigungsbewertung* 0.2 ist die Schwelle
    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:
        #Finden und finden Sie einen nächsten Nachbarn in Nearest Neighbors
        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))
        #Duplikate entfernen
        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')
        #Löschen Sie Datensätze, die den Schwellenwert überschreiten
        distance_df = distance_df[distance_df.distance < threshold]
        if len(distance_df) == 0:
            break
        #Extrahieren und löschen Sie übereinstimmende Datensätze
        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)

    #Gibt einen übereinstimmenden Datensatz zurück
    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]

Das Abgleichen von 1 Punkt in der Behandlungsgruppe mit 1 Punkt in der Kontrollgruppe mit dem nächstgelegenen Neigungswert wird wiederholt. Zu diesem Zeitpunkt wird ein Schwellenwert festgelegt, so dass nur Paare extrahiert werden, deren Abstand innerhalb des 0,2-fachen des Standardwerts liegt. Die Details dieses Bereichs wurden in Propensity Score Concept and Its Practice detailliert beschrieben. Ich vergleiche mit dem Übereinstimmungsergebnis von MatchIt in ch3_pscore.ipynb, aber es gibt keine genaue Übereinstimmung. Die Schlussfolgerung ist dieselbe und das Gleichgewicht der Kovariaten ist gut (siehe unten), daher ist es im Allgemeinen gut. Auch hier gibt es keine praktische Bibliothek wie Rs Cobalt Love.plot (), also visualisiere ich sie selbst. (Aus images / ch3_plot1.html)

Inverse wahrscheinlichkeitsgewichtete Schätzung (IPW)

Das Gute an IPW ist, dass es einfach zu implementieren ist. Dies wird auch mit WeightIt verglichen, scheint aber im Allgemeinen korrekt zu sein. (Auszug aus ch3_pscore.ipynb)

def get_ipw(X, y, random_state=0):
    #Berechnen Sie die Neigungsbewertung
    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('Es können nur 2 Gruppen abgeglichen werden. 2 Gruppen müssen sein[0, 1]Bitte mit ausdrücken.')
        raise ValueError
    # treatment ==1 bis Gruppe1, treatment ==Sei 0 Gruppe2.
    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 Wie in ch4_did.ipynb beschrieben, werden die folgenden zwei Bibliotheken zum Vergleich verwendet.

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

Da die Ausgabe wunderschön ist, verwende ich normalerweise die kausale Auswirkung von dafiti, aber es war tcassous kausaler Einfluss, dass das Schätzergebnis nahe am Buch lag (Rs ursprüngliche kausale Auswirkung). Beide scheinen durch das Zustandsraummodell der Statistikmodelle geschätzt zu werden, aber ich habe den Unterschied in der Implementierung nicht verstanden. Lass es mich wissen, bitte. .. In beiden Fällen ist der Schätzfehler im Vergleich zur R-Implementierung extrem klein. Es ist vielleicht nicht so gut, dass es viele Kovariaten gibt.

dafiti / kausale Wirkungskurve

Handlung von tcassou / causal_impact

Regressionsdiskontinuierliches Design (RDD)

Mit R können Sie es schnell mit rddtools ausführen, mit Python jedoch nicht. Überprüfen Sie zunächst die in rdd_reg_lm von rddtools verwendete Regressionsgleichung. (Referenz: https://cran.r-project.org/web/packages/rddtools/rddtools.pdf P23)

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

Nebenbei dachte ich, dass RDD Regressionsmodelle links und rechts vom Grenzwert erstellen und die Differenz zwischen den geschätzten Werten beim Grenzwert, aber eine Regression, verwenden würde. Es wird durch einen Ausdruck ausgedrückt. Ist es in Bezug auf die Bedeutung dasselbe? Wobei D eine binäre Variable mit einem Wert von 1 ist, wenn X größer oder gleich dem Grenzwert ist, und 0, wenn es früher ist. Der Coef von D ist der Wert, den Sie als Wirkungsgrad überprüfen möchten. Außerdem ist c der Grenzwert. Basierend auf dem oben Gesagten werden wir es implementieren. Ich verweise auch auf den Quellcode von rddtools zur Implementierung. (https://github.com/MatthieuStigler/RDDtools/blob/master/RDDtools/R/model.matrix.RDD.R) (Auszug aus ch5_rdd.ipynb)

class RDDRegression:
#R-Paket rddtools rdd_reg_Lm reproduzieren
#Referenz: 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)

Polynommerkmale von scicit-learn werden als Vorverarbeitung für die Durchführung einer Polynomregression verwendet. (Aus images / ch5_plot2_3.html) Ich konnte es gut einschätzen. Wie Sie anhand des Notizbuchs sehen können, entspricht der geschätzte Effekt fast dem des Buches. War gut.

nonparametric RDD (RDestimate) Ich konnte nicht ... Es ist fast der Teil von "(fast) in Python reproduziert" am Anfang. RDestimate verwendet die Methode von Immunos und Kalyanaraman (2012). Optimale Bandbreitenauswahl für den Regressionsdiskontinuitätsschätzer, um die optimale Bandbreite auszuwählen, aber ich war mir nicht sicher, wie ich diese Bandbreite schätzen sollte. In ch5_rdd.ipynb habe ich versucht, MSE zu implementieren, um die optimale Bandbreite zu finden, wenn ich die Bandbreite grob ändere, aber es scheint nicht sehr gut zu funktionieren. (Aus ch5_plot4.html) Hmmm, nicht wahr? Es scheint unmöglich zu sein, außerhalb der Bandbreite vorherzusagen, ob dies durch reine Verengung der Bandbreite geschätzt wird, aber wie machen Sie das tatsächlich? Oder vielleicht bin ich nur daran interessiert, die abgeschnittene Nachbarschaft zu schätzen, damit ich mir keine Sorgen machen muss.

abschließend

Wenn Sie Probleme mit Ihrem Verständnis oder Ihrer Implementierung haben oder etwas nicht stimmt, teilen Sie uns dies bitte mit. Twitter