Sådan laver man negativ binomialregression i Python
Vi starter med at importere alle de nødvendige pakker.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Næst skal vi oprette et pandas DataFrame for datasættet counts.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Vi tilføjer et par afledte regressionsvariabler til X-matrixen.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Vi vil ikke bruge variablen Date som en regressor, da den indeholder en absolut datoværdi, men vi behøver ikke at gøre noget særligt for at droppe Date, da den allerede forbruges som indeks i pandas DataFrame. Så den vil ikke være tilgængelig for os i X-matrixen.
Lad os oprette trænings- og testdatasættene.
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)))
STEG 1: Vi vil nu konfigurere og tilpasse Poisson-regressionsmodellen på træningsdatasættet.
Sæt regressionsudtrykket op i patsy notation. Vi fortæller patsy, at BB_COUNT er vores afhængige variabel, og at den afhænger af regressionsvariablerne: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T og PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Opret X- og y-matricer for trænings- og testdatasættene. patsy gør dette meget enkelt.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Træn Poisson-regressionsmodellen på træningsdatasættet ved hjælp af statsmodels GLM-klassen.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Dermed er træningen af Poisson-regressionsmodellen afsluttet. Hvis du vil se resultatet af træningen, kan du udskrive træningsoversigten.
print(poisson_training_results.summary())
Dette udskriver følgende:
Vores egentlige interesse ligger i vektoren af tilpassede satser λ, der er produceret ved træningen. Denne satsvektor er indeholdt i parameteren poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Det følgende output viser de første få og de sidste få værdier af den tilpassede λ-vektor:
Dette afslutter STEP1: tilpasning af Poisson-regressionsmodellen.
STEG 2: Vi tilpasser nu den supplerende OLS-regressionsmodel på datasættet og bruger den tilpassede model til at få værdien af α.
Importer api-pakken.
import statsmodels.formula.api as smf
Føj λ-vektoren som en ny kolonne kaldet “BB_LAMBDA” til datarammen for træningsdatasættet. Husk på, at λ’s dimensioner er (n x 1). I vores eksempel vil den være (161 x 1). Husk også på, at λ-vektoren er tilgængelig i poisson_training_results.mu :
df_train = poisson_training_results.mu
Næst skal vi tilføje en afledt kolonne kaldet “AUX_OLS_DEP” til pandas Data Frame. Denne nye kolonne vil gemme værdierne for den afhængige variabel i OLS-regressionen. Det er venstre side af OLS-regressionsligningen nedenfor:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
I ovenstående kodestump er den del med fed skrift den venstre side af ovenstående aux OLSR-ligning.
Lad os bruge patsy til at danne modelspecifikationen for OLSR. Vi ønsker at fortælle patsy, at AUX_OLS_DEP er den afhængige variabel, og at den forklares af BB_LAMBDA (som er satsvektoren λ). “-1” i slutningen af udtrykket er patsy-syntaks for at sige: brug ikke et regressionsintercept; dvs. tilpas blot en ret linje, der går gennem oprindelsen, som foreslået af de herrer Cameron og Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Vi er nu klar til at tilpasse en OLSR-model.
Konfigurer og tilpas OLSR-modellen:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Udskriv regressionsparametrene:
print(aux_olsr_results.params)
Du vil se følgende enkeltkoefficient blive udskrevet, der svarer til den enkelte regressionsvariabel BB_LAMBDA. Denne koefficient er den α, som vi søgte:
BB_LAMBDA 0.037343
Er α statistisk signifikant?
Vi skal nu besvare et meget vigtigt spørgsmål. Er denne værdi af α (0,037343) statistisk signifikant? Eller kan den for alle praktiske formål betragtes som nul?
Hvorfor er det så vigtigt at finde ud af dette? Husk på, at hvis α er nul, så er følgende ligning:
…reduceres til Varians = middelværdi. Dette er variansfunktionen for Poisson-regressionsmodellen.
Hvis værdien af α ikke er statistisk signifikant, kan den negative binomiale regressionsmodel ikke passe bedre til træningsdatasættet end en Poisson-regressionsmodel.
Objektet OLSResults indeholder t-score for regressionskoefficienten α. Lad os udskrive den:
aux_olsr_results.tvalues
Dette udskriver:
BB_LAMBDA 4.814096
Fra en t-værdiberegner kan vi se, at den kritiske t-værdi ved et konfidensniveau på 99 % (højrehalet) og frihedsgrader=(172 observationer) -( 6 regressionsvariabler)- (1 dispersionsparameter α)=165 er 2,34916. Dette er komfortabelt mindre end t-statistikken for α, som var 4,814096. Vi konkluderer, at,
α=0,037343 er statistisk signifikant.
Dette afslutter TRIN 2: Bestemmelse af α.
TRIN 3: Vi leverer værdien af alpha fundet i TRIN 2 til statsmodels.genmod.families.family.NegativeBinomial-klassen og træner NB2-modellen på træningsdatasættet.
Dette er en operation i ét trin i statsmodeller:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Som før udskriver vi træningsoversigten:
print(nb2_training_results.summary())
Hvilket udskriver følgende oversigt:
STEP 4: Lad os lave nogle forudsigelser ved hjælp af vores trænede NB2-model.
Forudsigelse er igen en enkelt-trin procedure i statsmodeller:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Lad os udskrive forudsigelserne:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Følgende er de første par linjer i output:
Lad os også plotte de forudsagte tællinger i forhold til de faktiske tællinger for testdataene.
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()
Her er outputtet:
Ingen særlig dårlig! NB2-modellen ser ud til mere eller mindre at følge tendensen i cykeltællingerne. Og ligesom med Poisson-regressionsmodellens resultater er dens forudsigelser i nogle tilfælde langt fra de faktiske værdier.