How to do Negative Binomial Regression in Python
Aloitetaan tuomalla kaikki tarvittavat paketit.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Luotaan seuraavaksi pandas DataFrame laskentatietoaineistolle.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Lisätään X-matriisiin muutama johdettu regressiomuuttuja.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Emme käytä Date-muuttujaa regressorina, koska se sisältää absoluuttisen päivämäärän arvon, mutta meidän ei tarvitse tehdä mitään erityistä pudottaaksemme Date-muuttujan pois, koska se on jo kulutettu pandas DataFramen indeksinä. Se ei siis ole käytettävissämme X-matriisissa.
Luotaan harjoittelu- ja testausdatajoukot.
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)))
VAIHE 1: Määritämme nyt Poisson-regressiomallin ja sovitamme sen harjoitteludatajoukkoon.
Määritellään regressioilmaus patsy-merkintätavalla. Kerromme patsylle, että BB_COUNT on riippuvainen muuttujamme ja että se riippuu regressiomuuttujista: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T ja PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Asetetaan X- ja y-matriisit harjoittelu- ja testausdatajoukoille. patsy tekee tästä todella yksinkertaista.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Kouluta Poisson-regressiomalli statsmodels GLM-luokan avulla harjoitusdatajoukolla.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Tämä lopettaa Poisson-regressiomallin kouluttamisen. Jos haluat nähdä harjoittelun tuloksen, voit tulostaa harjoittelun yhteenvedon.
print(poisson_training_results.summary())
Tämä tulostaa seuraavan tulosteen:
>
Vaikutteinen kiinnostuksen kohteemme on harjoittelun tuloksena syntynyt sovitettujen kertoimien λ vektori. Tämä korkovektori sisältyy parametriin poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Seuraavassa tulosteessa näkyvät sovitetun λ-vektorin ensimmäiset ja viimeiset arvot:
Näillä on suoritettu VAIHE1: Poissonin regressiomallin sovittaminen.
VAIHE 2: Sovitamme nyt ylimääräisen OLS-regressiomallin aineistoon ja käytämme sovitettua mallia saadaksemme α:n arvon.
Importoi api-paketti.
import statsmodels.formula.api as smf
Lisää λ-vektori uudeksi sarakkeeksi nimeltä ’BB_LAMBDA’ harjoitusdatajoukon datakehykseen. Muista, että λ:n dimensiot ovat (n x 1). Esimerkissämme se on (161 x 1). Muistetaan myös, että λ-vektori on käytettävissä poisson_training_results.mu :
df_train = poisson_training_results.mu
Seuraavaksi lisätään johdettu sarake nimeltä ’AUX_OLS_DEP’ pandas Data Frameen. Tähän uuteen sarakkeeseen tallennetaan OLS-regression riippuvan muuttujan arvot. Se on alla olevan OLS-regressioyhtälön vasen puoli:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
Ylläolevassa koodin katkelmassa lihavoituna oleva osa on ylläolevan aux-OLSR-epätasapainoyhtälön vasen puoli.
Käytetään patsya OLSR:n mallispesifikaation muodostamiseen. Haluamme kertoa patsylle, että AUX_OLS_DEP on riippuvainen muuttuja ja sitä selittää BB_LAMBDA (joka on korkovektori λ). ’-1’ lausekkeen lopussa on patsyn syntaksi, jolla sanotaan: älä käytä regression leikkauspistettä; eli sovita vain origon kautta kulkeva suora, kuten herrat Cameron ja Trivedi ovat ehdottaneet.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Olemme nyt valmiita sovittamaan OLSR-mallin.
Konfiguroi ja sovita OLSR-malli:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Tulosta regressioparametrit:
print(aux_olsr_results.params)
Näet seuraavan yksittäisen kertoimen tulostuvan, joka vastaa yksittäistä regressiomuuttujaa BB_LAMBDA. Tämä kerroin on α, jota etsimme:
BB_LAMBDA 0.037343
Onko α tilastollisesti merkitsevä?
Meidän on nyt vastattava hyvin tärkeään kysymykseen. Onko tämä α:n arvo (0,037343) tilastollisesti merkitsevä? Vai voidaanko sitä pitää käytännössä nollana?
Miksi tämän selvittäminen on niin tärkeää? Muista, että jos α on nolla, seuraava yhtälö:
…redusoituu muotoon Varianssi = keskiarvo. Tämä on Poisson-regressiomallin varianssifunktio.
Jos α:n arvo ei ole tilastollisesti merkitsevä, negatiivisen binomialueen regressiomalli ei pysty sovittamaan harjoitusaineistoa paremmin kuin Poisson-regressiomalli.
OLSResults-olio sisältää regressiokertoimen α t-arvon. Tulostetaan se:
aux_olsr_results.tvalues
Tulostuu:
BB_LAMBDA 4.814096
Tulostetaan t-arvolaskimella, että kriittinen t-arvo 99 %:n luottamustasolla (oikeanpuoleinen) ja vapausasteilla=(172 havaintoa)-( 6 regressiomuuttujaa)- (1 hajontaparametri α)=165 on 2,34916. Tämä on huomattavasti pienempi kuin α:n t-statistiikka, joka oli 4,814096. Päättelemme, että,
α=0,037343 on tilastollisesti merkitsevä.
Tämä päättää VAIHEEN 2: α:n määrittäminen.
VAIHE 3: Toimitamme VAIHEESSA 2 löydetyn alfa-arvon statsmodels.genmod.families.family.NegativeBinomial-luokkaan ja harjoittelemme NB2-mallia harjoitteluaineistolla.
Tämä on yksivaiheinen operaatio statsmodelsissa:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Tulostamme kuten ennenkin koulutuksen yhteenvedon:
print(nb2_training_results.summary())
Joka tulostaa seuraavan yhteenvedon:
VAIHE 4: Tehdään joitakin ennusteita käyttäen koulutettua NB2-malliamme.
Ennustaminen on taas yksivaiheinen menettely statsmodelsissa:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Tulostetaan ennusteet:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Seuraavat ovat ensimmäiset rivit tulosteesta:
Piirretään myös ennustetut lukumäärät suhteessa todellisiin lukumääriin testidatan osalta.
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()
Tässä on tuloste:
Ei hassummin! NB2-malli näyttää enemmän tai vähemmän seuraavan polkupyörälaskentojen suuntausta. Ja aivan kuten Poisson-regressiomallin suorituskyvyn kohdalla, joissakin tapauksissa sen ennusteet ovat kaukana todellisista arvoista
.