Negative Binomialregression in Python
Zunächst importieren wir alle erforderlichen Pakete.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Nachfolgend erstellen wir einen Pandas DataFrame für den Zähldatensatz.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Wir fügen der X-Matrix einige abgeleitete Regressionsvariablen hinzu.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Wir werden die Variable Datum nicht als Regressor verwenden, da sie einen absoluten Datumswert enthält, aber wir brauchen nichts Besonderes zu tun, um das Datum wegzulassen, da es bereits als Index des Pandas DataFrame verwendet wird. Es wird uns also nicht in der X-Matrix zur Verfügung stehen.
Lassen Sie uns die Trainings- und Testdatensätze erstellen.
mask = np.random.rand(len(df)) < 0.8
df_train = df
df_test = df
print('Training data set length='+str(len(df_train)))
print('Testing data set length='+str(len(df_test)))
STEP 1: Wir werden nun das Poisson-Regressionsmodell auf dem Trainingsdatensatz konfigurieren und anpassen.
Richten Sie den Regressionsausdruck in der Patsy-Notation ein. Wir teilen patsy mit, dass BB_COUNT unsere abhängige Variable ist und dass sie von den Regressionsvariablen abhängt: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T und PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Richten Sie die X- und y-Matrizen für die Trainings- und Testdatensätze ein. patsy macht dies wirklich einfach.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Verwenden Sie die GLM-Klasse von statsmodels und trainieren Sie das Poisson-Regressionsmodell auf dem Trainingsdatensatz.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Damit ist das Training des Poisson-Regressionsmodells abgeschlossen. Um das Ergebnis des Trainings zu sehen, können Sie die Trainingszusammenfassung ausdrucken.
print(poisson_training_results.summary())
Das gibt Folgendes aus:
Unser eigentliches Interesse gilt dem Vektor der angepassten Raten λ, der durch das Training erzeugt wurde. Dieser Ratenvektor ist im Parameter poisson_training_results.mu enthalten.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Die folgende Ausgabe zeigt die ersten und letzten Werte des angepassten λ-Vektors:
Damit ist Schritt 1 abgeschlossen: die Anpassung des Poisson-Regressionsmodells.
SCHRITT 2: Wir passen nun das OLS-Hilfsregressionsmodell an den Datensatz an und verwenden das angepasste Modell, um den Wert von α zu erhalten.
Importieren Sie das api-Paket.
import statsmodels.formula.api as smf
Fügen Sie den λ-Vektor als neue Spalte mit dem Namen „BB_LAMBDA“ zum Datenrahmen des Trainingsdatensatzes hinzu. Erinnern Sie sich, dass die Dimensionen von λ (n x 1) sind. In unserem Beispiel ist es (161 x 1). Denken Sie auch daran, dass der λ-Vektor in poisson_training_results.mu :
df_train = poisson_training_results.mu
verfügbar ist. Als Nächstes fügen wir dem Pandas Data Frame eine abgeleitete Spalte mit dem Namen ‚AUX_OLS_DEP‘ hinzu. In dieser neuen Spalte werden die Werte der abhängigen Variablen der OLS-Regression gespeichert. Es handelt sich um die linke Seite der folgenden OLS-Regressionsgleichung:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
Im obigen Codeschnipsel ist der fettgedruckte Teil die linke Seite der obigen aux OLSR Gleichung.
Lassen Sie uns patsy verwenden, um die Modellspezifikation für OLSR zu bilden. Wir wollen patsy mitteilen, dass AUX_OLS_DEP die abhängige Variable ist und dass sie durch BB_LAMBDA (den Ratenvektor λ) erklärt wird. Das ‚-1‘ am Ende des Ausdrucks ist die Patsy-Syntax, die besagt, dass kein Regressionsabschnitt verwendet werden soll, d.h. es soll einfach eine gerade Linie durch den Ursprung angepasst werden, wie von den Herren Cameron und Trivedi vorgeschlagen.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Wir sind nun bereit, ein OLSR-Modell anzupassen.
Konfigurieren Sie das OLSR-Modell und passen Sie es an:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Drucken Sie die Regressionsparameter:
print(aux_olsr_results.params)
Sie werden sehen, dass der folgende einzelne Koeffizient ausgedruckt wird, der der einzelnen Regressionsvariablen BB_LAMBDA entspricht. Dieser Koeffizient ist das α, nach dem wir gesucht haben:
BB_LAMBDA 0.037343
Ist α statistisch signifikant?
Wir müssen nun eine sehr wichtige Frage beantworten. Ist dieser Wert von α (0,037343) statistisch signifikant? Oder kann er für alle praktischen Zwecke als Null betrachtet werden?
Warum ist es so wichtig, dies herauszufinden? Erinnern Sie sich, dass, wenn α Null ist, die folgende Gleichung:
…reduziert sich auf Varianz = Mittelwert. Dies ist die Varianzfunktion des Poisson-Regressionsmodells.
Wenn der Wert von α statistisch nicht signifikant ist, dann kann das Negativ-Binomial-Regressionsmodell den Trainingsdatensatz nicht besser anpassen als ein Poisson-Regressionsmodell.
Das OLSResults-Objekt enthält den t-Score des Regressionskoeffizienten α. Drucken wir ihn aus:
aux_olsr_results.tvalues
Das Ergebnis lautet:
BB_LAMBDA 4.814096
Aus einem t-Wert-Rechner können wir ersehen, dass der kritische t-Wert bei einem Konfidenzniveau von 99 % (rechtsseitig) und Freiheitsgraden=(172 Beobachtungen) -( 6 Regressionsvariablen)- (1 Dispersionsparameter α)=165 2,34916 beträgt. Dies ist deutlich weniger als die t-Statistik von α, die 4,814096 beträgt. Wir schlussfolgern, dass,
α=0,037343 statistisch signifikant ist.
Damit ist Schritt 2 abgeschlossen: Die Bestimmung von α.
Schritt 3: Wir geben den in Schritt 2 gefundenen Wert von alpha in die Klasse statsmodels.genmod.families.family.NegativeBinomial ein und trainieren das NB2-Modell auf dem Trainingsdatensatz.
Dies ist ein Ein-Schritt-Vorgang in statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Wie zuvor drucken wir die Trainingszusammenfassung:
print(nb2_training_results.summary())
Die folgende Zusammenfassung wird gedruckt:
SCHRITT 4: Lassen Sie uns einige Vorhersagen mit unserem trainierten NB2-Modell machen.
Die Vorhersage ist wieder ein einstufiges Verfahren in statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Lassen Sie uns die Vorhersagen ausdrucken:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Nachfolgend sind die ersten Zeilen der Ausgabe:
Lassen Sie uns auch die vorhergesagten Zählungen gegenüber den tatsächlichen Zählungen für die Testdaten darstellen.
predicted_counts=predictions_summary_frameactual_counts = y_testfig = plt.figure()fig.suptitle('Predicted versus actual bicyclist counts on the Brooklyn bridge')predicted, = plt.plot(X_test.index, predicted_counts, 'go-', label='Predicted counts')actual, = plt.plot(X_test.index, actual_counts, 'ro-', label='Actual counts')plt.legend(handles=)plt.show()
Hier ist die Ausgabe:
Nicht so schlecht! Das NB2-Modell scheint mehr oder weniger dem Trend bei den Fahrradzählungen zu folgen. Und genau wie bei der Leistung des Poisson-Regressionsmodells liegen seine Vorhersagen in einigen Fällen weit von den tatsächlichen Werten entfernt.