How to do Negative Binomial Regression in Python
Az összes szükséges csomag importálásával kezdjük.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Következő lépésként létrehozunk egy pandas DataFrame-et a számlálási adathalmazhoz.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Az X mátrixhoz hozzáadunk néhány származtatott regressziós változót.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
A Date változót nem fogjuk regresszornak használni, mivel abszolút dátumértéket tartalmaz, de a Date elhagyásához nem kell semmi különöset tennünk, mivel a pandas DataFrame indexeként már fel van használva. Így nem lesz elérhető számunkra az X mátrixban.
Létrehozzuk a gyakorló és teszt adathalmazt.
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)))
1. LÉPÉS: Most a Poisson regressziós modellt konfiguráljuk és illesztjük a gyakorló adathalmazon.
Elkészítjük a regressziós kifejezést patsy jelöléssel. Megmondjuk a patsy-nek, hogy a BB_COUNT a függő változónk, és függ a regressziós változóktól: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T és PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Beállítjuk az X és y mátrixokat a képzési és tesztelési adathalmazokhoz. patsy ezt nagyon egyszerűvé teszi.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
A statsmodels GLM osztály segítségével képezze ki a Poisson-regressziós modellt a gyakorló adathalmazon.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Ezzel befejeződik a Poisson-regressziós modell képzése. A képzés eredményének megtekintéséhez kinyomtathatja a képzési összefoglalót.
print(poisson_training_results.summary())
Ez a következőket nyomtatja ki:
Az igazi érdeklődésünk a képzés által előállított λ illesztett ráták vektorára irányul. Ezt az arányvektort a poisson_training_results.mu paraméter tartalmazza.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
A következő kimenet az illesztett λ vektor első és utolsó néhány értékét mutatja:
Ezzel befejeződött az 1. LÉPÉS: a Poisson-regressziós modell illesztése.
2. LÉPÉS: Most a segéd OLS regressziós modellt illesztjük az adathalmazra, és az illesztett modell segítségével megkapjuk α értékét.
Importáljuk az api csomagot.
import statsmodels.formula.api as smf
Adjuk hozzá a λ vektort egy új, ‘BB_LAMBDA’ nevű oszlopként a képzési adathalmaz adatkeretéhez. Emlékezzünk arra, hogy a λ dimenziói (n x 1). Példánkban ez (161 x 1) lesz. Emlékezzünk arra is, hogy a λ vektor a poisson_training_results.mu :
df_train = poisson_training_results.mu
A következő lépésben adjunk hozzá egy származtatott oszlopot ‘AUX_OLS_DEP’ néven a pandas Data Frame-hez. Ez az új oszlop az OLS-regresszió függő változójának értékeit fogja tárolni. Ez az alábbi OLS-regressziós egyenlet bal oldala:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
A fenti kódrészletben a félkövérrel szedett rész a fenti aux OLSR egyenlet bal oldala.
A patsy segítségével alakítsuk ki az OLSR modellspecifikációját. Azt akarjuk megmondani a patsy-nek, hogy az AUX_OLS_DEP a függő változó, és azt a BB_LAMBDA magyarázza (ami a λ sebességvektor). A ‘-1’ a kifejezés végén a patsy szintaxisa, amely azt mondja: ne használjunk regressziós metszetet; azaz csak egy, az origón áthaladó egyenest illesszünk, ahogy azt Cameron és Trivedi urak javasolják.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Most már készen állunk az OLSR modell illesztésére.
Konfiguráljuk és illesszük az OLSR modellt:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Kinyomtassuk a regressziós paramétereket:
print(aux_olsr_results.params)
A következő egyetlen együtthatót fogjuk látni, amely a BB_LAMBDA egyetlen regressziós változónak felel meg. Ez az együttható az α, amit kerestünk:
BB_LAMBDA 0.037343
Az α statisztikailag szignifikáns?
Most egy nagyon fontos kérdésre kell válaszolnunk. Ez az α érték (0,037343) statisztikailag szignifikáns? Vagy minden gyakorlati szempontból nullának tekinthető?
Miért olyan fontos ezt kideríteni? Emlékezzünk vissza, hogy ha α nulla, akkor a következő egyenlet:
…redukálódik a Variancia = átlagra. Ez a Poisson-regressziós modell varianciafüggvénye.
Ha α értéke statisztikailag nem szignifikáns, akkor a negatív binomiális regressziós modell nem tud jobban illeszkedni a gyakorló adathalmazra, mint a Poisson-regressziós modell.
Az OLSResults objektum tartalmazza az α regressziós együttható t-értékét. Nyomtassuk ki:
aux_olsr_results.tvalues
Ez kiírja:
BB_LAMBDA 4.814096
A t-érték kalkulátorból láthatjuk, hogy a kritikus t-érték 99%-os megbízhatósági szinten (jobb oldali farkú), és szabadságfokok=(172 megfigyelés)-( 6 regressziós változó)-(1 szórásparaméter α)=165 mellett 2,34916. Ez kényelmesen kisebb, mint az α t-statisztikája, amely 4,814096. Megállapítjuk,
α=0,037343 statisztikailag szignifikáns.
Ezzel befejeződik a 2. LÉPÉS: Az α meghatározása.
3. LÉPÉS: A 2. LÉPÉSBEN talált alfa értéket a statsmodels.genmod.families.family.NegativeBinomial osztályba adjuk, és az NB2 modellt a képzési adathalmazon betanítjuk.
Ez egylépéses művelet a statsmodelsben:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Az előzőekhez hasonlóan kiírjuk a képzési összegzést:
print(nb2_training_results.summary())
Az alábbi összegzést írja ki:
4. Lépés: Készítsünk néhány előrejelzést a betanított NB2 modellünkkel.
A jóslás ismét egylépéses eljárás a statsmodelsben:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Nyomtassuk ki a jóslatokat:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
A következőkben a kimenet első néhány sora következik:
A tesztadatok esetében ábrázoljuk a megjósolt számokat is a tényleges számokkal szemben.
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()
Itt a kimenet:
Nem rossz! Úgy tűnik, hogy az NB2 modell többé-kevésbé követi a kerékpárosok számának tendenciáját. És csakúgy, mint a Poisson-regressziós modell teljesítménye esetében, néhány esetben az előrejelzései messze eltérnek a tényleges értékektől
.