Lernen Sie APPELPY in der 10 Minutes-Python Applied Metric Economics Library

** In 10 Minuten lernen APPEL PY **

Dieser Beitrag ist eine japanische Übersetzung von diesem Notizbuch.

import pandas as pd
import numpy as np
import statsmodels.api as sm  # for datasets
import matplotlib.pyplot as plt
%matplotlib inline

import warnings
warnings.filterwarnings('ignore')

** Einführung in die Funktionen von Appelpy. ** ** **

• Explorative Datenanalyse
• Modellschätzung
• Modelldiagnose

Hinweis: Dieses Notizbuch wurde mit der neuesten Version von Appelpy v0.4.2 und seinen Abhängigkeiten ausgeführt.

Modelldatensatzeinstellungen

Lassen Sie uns nun die Funktionen von Appelpy anhand eines der umfangreichsten Datensätze in der Metrikökonomie untersuchen. Der California Test Score-Datensatz (** Caschool **).

Laden Sie den Datensatz mit Pandas und rufen Sie die Daten mit der Statsmodels-Methode ab.

# Load data
df = sm.datasets.get_rdataset('Caschool', 'Ecdat').data
# Add together the two test scores (total score)
df['scrtot'] = df['readscr'] + df['mathscr']

df.head()

Überprüfen Sie die Daten

Nachdem wir die Daten haben, schauen wir uns die Verteilung einiger Schlüsselvariablen vor der Modellierung genauer an.

Appelpy hat ein ** eda Modul ** (explorative Datenanalyse).

from appelpy.eda import statistical_moments

Holen Sie sich alle statistischen Momente auf einmal.

1. Durchschnitt
2. Verteilt
3. Verzerrung (Standard ist Fisher)
4. Schärfe

statistical_moments(df)

Modellschätzung

Im Zentrum von Appelpy steht die Modellierung der Vorgänge im Datensatz.

Das Modul ** linear_model ** enthält Klassen zum Modellieren von Daten. Beispiel: OLS (Ordinary Least Squares)

from appelpy.linear_model import OLS

** Hier ist ein Rezept zur Schätzung eines linearen Modells. ** ** **

1. Erstellen Sie zwei Listen, "y_list" und "X_list", für die abhängige Variable (y) und die unabhängige Variable (X).

y_list = ['scrtot']
X_list = ['avginc', 'str', 'enrltot', 'expnstu']

2. Passen Sie das Modell in eine Zeile an, indem Sie den Modelldatenrahmen und zwei Listen als Eingabe für das OLS-Objekt übergeben.

model_nonrobust = OLS(df, y_list, X_list).fit()

3. Ermitteln Sie das Ergebnis des Schlüsselmodells mit dem Attribut model_selection_stats, einschließlich des quadratischen mittleren Quadratwurzelfehlers (Root MSE).

model_nonrobust.model_selection_stats

{'root_mse': 25.85923055850994,
 'r_squared': 0.5438972040191434,
 'r_squared_adj': 0.5395010324916171,
 'aic': 3929.1191674301117,
 'bic': 3949.320440986499}

Holen Sie sich eine Zusammenfassung der Statistikmodelle mit dem Attribut results_output ...

model_nonrobust.results_output

... Jetzt ist es endlich einfach, die Ergebnisse eines standardisierten Modells mit dem Attribut results_output_standardized zu erhalten.

model_nonrobust.results_output_standardized
# This is a Styler object; add .data to access the underlying Pandas dataframe

Die Ausgabe ähnelt dem Befehl listcoef von Stata.

Modelldiagnose

Rufen Sie die Methode diagnostic_plot des Modellobjekts auf, um viele visuelle ** Regressionsdiagnosen anzuzeigen: **.

--pp_plot: P-P-Plot -- qq_plot: Q-Q-Plot --rvp_plot: Rest- und Vorhersagewertdiagramm --rvf_plot: Darstellung der Rest- und Messwerte

Appelpys ** Diagnose ** -Modul verfügt über eine Reihe von Klassen und Funktionen für die Regressionsdiagnose.

• Wenn die visuelle Diagnose nicht ausreicht, rufen Sie eine andere Methode wie "heteroskedasticity_test" auf, um einen formalen Test durchzuführen.

• Gibt es Beobachtungen, die das Modell verschmutzen? Übergeben Sie das Modellobjekt an die Klasse "BadApples", um die Einflusspunkte, Ausreißer und die hohe Hebelwirkung zu ermitteln.

Zeichnung

** Hier ist ein Rezept für einen 2x2-Satz von Diagnoseplots. ** ** **

fig, axarr = plt.subplots(2, 2, figsize=(10, 10))
model_nonrobust.diagnostic_plot('pp_plot', ax=axarr[0][0])
model_nonrobust.diagnostic_plot('qq_plot', ax=axarr[1][0])
model_nonrobust.diagnostic_plot('rvf_plot', ax=axarr[1][1])
axarr[0, 1].axis('off')
plt.tight_layout()

Guter Kummer. Das angegebene Modell hat ein spezielles Problem mit Testergebnissen über 1350.

Es scheint eine nicht konstante Varianz im RVF-Diagramm zu geben.

Es ist wahrscheinlich gleichmäßig verteilt, aber lassen Sie uns einen genaueren Test durchführen.

Gleiche Dispersion

from appelpy.diagnostics import heteroskedasticity_test

Rufen Sie einen von mehreren ** Isodispersitätstests ** auf. Beispiel:

• • Breusch-Pagan * (Standardeinstellung "Hettest" des Statas)
• • Bresuch-Pagan studentisiert * (Standard für Rs Befehl "bptest")
• • Weiß * (Statas "heißester, weißer" Befehl)

bp_stats = heteroskedasticity_test('breusch_pagan', model_nonrobust)
print('Breusch-Pagan test :: {}'.format(bp_stats['distribution'] + '({})'.format(bp_stats['nu'])))
print('Test statistic: {:.4f}'.format(bp_stats['test_stat']))
print('Test p-value: {:.4f}'.format(bp_stats['p_value']))

Breusch-Pagan test :: chi2(1)
Test statistic: 2.5579
Test p-value: 0.1097

bps_stats = heteroskedasticity_test('breusch_pagan_studentized', model_nonrobust)
print('Breusch-Pagan test (studentized) :: {}'.format(bps_stats['distribution'] + '({})'.format(bps_stats['nu'])))
print('Test statistic: {:.4f}'.format(bps_stats['test_stat']))
print('Test p-value: {:.4f}'.format(bps_stats['p_value']))

Breusch-Pagan test (studentized) :: chi2(1)
Test statistic: 11.3877
Test p-value: 0.0225

Der studentisierte Breusch-Pagan-Test gilt als robuster und zeigt signifikante Unterschiede zur ursprünglichen Hypodispersitätshypothese. Passen wir ihn daher an ein Modell mit dem heteroskedastisch-robusten Standardfehler (HC1-Fehler) an. ..

model_hc1 = OLS(df, y_list, X_list, cov_type='HC1').fit()

model_hc1.results_output

Was sind jetzt die wesentlichen Vergeltungsmaßnahmen (unabhängige Variablen)?

model_hc1.significant_regressors(0.05)

['avginc', 'enrltot', 'expnstu']

Die Rückkehrer von "str" (Schüler-Lehrer-Verhältnis) und "expnstu" sind besonders relevant für Investitionen in Schüler. Sie können einen ** Wald-Test ** durchführen, um zu testen, ob die beiden Vergeltungsmaßnahmen zusammen sind und sich erheblich von 0 unterscheiden.

from appelpy.diagnostics import wald_test

wald_stats = wald_test(model_hc1, ['str', 'expnstu'])
print('Wald test :: {}'.format(wald_stats['distribution'] + '({})'.format(wald_stats['nu'])))
print('Test statistic: {:.4f}'.format(wald_stats['test_stat']))
print('Test p-value: {:.4f}'.format(wald_stats['p_value']))

Wald test :: chi2(2)
Test statistic: 4.7352
Test p-value: 0.0937

Es überrascht nicht, dass der Standardfehler einiger Variablen beim HC1-Fehler höher ist als beim nicht robusten Fehler.

pd.concat([pd.Series(model_nonrobust.results.bse, name='std_error_nonrobust'),
           pd.Series(model_hc1.results.bse, name='std_error_hc1')], axis='columns')

Anomale Beobachtungen: Hebelwirkung, Ausreißer und Einfluss

Lassen Sie uns die beobachteten Ergebnisse untersuchen.

• Hoher Einfluss
• Hohe Hebelwirkung
• Früher

from appelpy.diagnostics import BadApples

Richten Sie eine Instanz von "BadApples" mithilfe eines Modellobjekts ein.

bad_apples = BadApples(model_hc1).fit()

Zeichnen wir den Hebelwert gegen das Quadrat des normalisierten Residuums.

Die Beobachtungen im oberen rechten Quadranten weisen überdurchschnittlich hohe Hebel- und Restwerte auf (dh die einflussreichsten Beobachtungen im Modell).

fig, ax = plt.subplots(figsize=(11,7))
bad_apples.plot_leverage_vs_residuals_squared(annotate=True)
plt.show()

Hebelwirkung

Der Hebelwert wird im Attribut Measures_leverage gespeichert.

bad_apples.measures_leverage.sort_values(ascending=False).head()

48     0.113247
34     0.085633
404    0.076762
413    0.065741
159    0.061353
Name: leverage, dtype: float64

Werfen wir einen Blick auf die fünf am häufigsten verwendeten Beobachtungen.

bad_apples.show_extreme_observations().loc[[48,34,404,413,159]]

model_hc1.X_standardized.loc[[48,34,404,413,159]]

Diese hoch verschuldeten Punkte haben Werte, bei denen mindestens eine unabhängige Variable mindestens 4 Standardabweichungen vom Mittelwert aufweist.

Ausreißer

Das BadApples-Objekt berechnet eine Heuristik, um festzustellen, ob der Skalierungswert "extrem" ist.

Extreme Indizes werden in Attributen gespeichert, beispielsweise werden Ausreißer in "indices_outliers" gespeichert.

bad_apples.indices_outliers

[5, 6, 7, 8, 9, 10, 13, 14, 37, 38, 47, 134, 370, 392, 403, 404, 415, 417]

Ausreißerheuristik: Wenn der Absolutwert des standardisierten Residuums oder der Absolutwert des studentisierten Residuums größer als 2 ist, ist der beobachtete Wert ein Ausreißer.

print(f"Proportion of observations as outliers: {len(bad_apples.indices_outliers) / len(df):.4%}")

Proportion of observations as outliers: 4.2857%

Aufgrund seiner Heuristik wird erwartet, dass weniger als etwa 5% der beobachteten Werte als Ausreißer gelten.

Beeinflussen

Berücksichtigen Sie auch verschiedene Einflussmessungen wie Kochabstand, Dffits und Dfbeta.

bad_apples.measures_influence.sort_values('cooks_d', ascending=False).head()

einpacken

--Laden Sie den Datensatz und verarbeiten Sie ihn mit Pandas vor. --Appelpy ** eda **: Überprüfen Sie den Modelldatensatz mit nützlichen Funktionen wie statistischen Momenten. --Appelpy ** Modellierung **. Linear_model enthält Klassen wie OLS und WLS zur Modellschätzung. Verwenden Sie für jedes Modell ein kurzes Rezept, um auf Modellschätzungen zuzugreifen. --Appelpy ** Diagnose **: Überprüft, ob Modellschätzungen mit Diagnosediagrammen, statistischen Tests und Diagnosen wie BadApples für die Auswirkungsanalyse gültig sind.

Lassen Sie uns auch andere Funktionen überprüfen.

--Logistische Regression (unter discrete_model) --DummyEncoder- und InteractionEncoder-Klassen