How to do Negative Binomial Regression in Python
We beginnen met het importeren van alle benodigde pakketten.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Volgende, maak een pandas DataFrame voor de tellingen dataset.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
We zullen een paar afgeleide regressie variabelen toevoegen aan de X matrix.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
We zullen de variabele Date niet als regressor gebruiken omdat deze een absolute datumwaarde bevat, maar we hoeven niets speciaals te doen om Date te laten vallen omdat deze al wordt gebruikt als de index van het pandas DataFrame. Het zal dus niet beschikbaar zijn voor ons in de X matrix.
Laten we de training en testing data sets creëren.
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)))
STAP 1: We zullen nu het Poisson regressie model configureren en fitten op de training data set.
Opzetten van de regressie expressie in patsy notatie. We vertellen patsy dat BB_COUNT onze afhankelijke variabele is en dat deze afhangt van de regressievariabelen: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T en PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Opzetten van de X en y matrices voor de training en testing data sets. patsy maakt dit heel eenvoudig.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Gebruik de statsmodels GLM klasse, train het Poisson regressie model op de training data set.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Dit is het einde van de training van het Poisson regressie model. Om het resultaat van de training te zien, kunt u de trainingssamenvatting afdrukken.
print(poisson_training_results.summary())
Dit drukt het volgende af:
Onze echte belangstelling gaat uit naar de vector van gepaste percentages λ die door de training wordt geproduceerd. Deze koersvector is opgenomen in de parameter poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
De volgende uitvoer toont de eerste paar en de laatste paar waarden van de gepaste λ-vector:
Dit voltooit STAP1: het passen van het Poisson-regressiemodel.
STAP 2: We passen nu het OLS-regressiemodel op de gegevensreeks toe en gebruiken het model om de waarde van α te verkrijgen.
Importeer het api-pakket.
import statsmodels.formula.api as smf
Voeg de λ-vector als een nieuwe kolom met de naam ‘BB_LAMBDA’ toe aan het gegevensframe van de trainingsgegevensreeks. Bedenk dat de afmetingen van λ (n x 1) zijn. In ons voorbeeld zal dat (161 x 1) zijn. Bedenk ook dat de λ-vector beschikbaar is in poisson_training_results.mu :
df_train = poisson_training_results.mu
Next, laten we een afgeleide kolom genaamd ‘AUX_OLS_DEP’ toevoegen aan het pandas Data Frame. Deze nieuwe kolom zal de waarden van de afhankelijke variabele van de OLS regressie opslaan. Het is het linkerdeel van onderstaande OLS-regressievergelijking:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
In het bovenstaande codefragment is het vetgedrukte gedeelte het linkerdeel van bovenstaande aux OLSR-vergelijking.
Laten we patsy gebruiken om de modelspecificatie voor de OLSR op te stellen. Wij willen patsy vertellen dat AUX_OLS_DEP de afhankelijke variabele is en dat deze wordt verklaard door BB_LAMBDA (die de koersvector λ is). De “-1” aan het eind van de uitdrukking is patsy’s syntaxis om te zeggen: gebruik geen regressie-intercept; pas dus gewoon een rechte lijn die door de oorsprong gaat, zoals voorgesteld door de heren Cameron en Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
We zijn nu klaar om een OLSR model te passen.
Configureer en fit het OLSR-model:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Print de regressieparameters:
print(aux_olsr_results.params)
U zult de volgende enkele coëfficiënt afgedrukt zien die correspondeert met de enkele regressievariabele BB_LAMBDA. Deze coëfficiënt is de α waarnaar we op zoek waren:
BB_LAMBDA 0.037343
Is α statistisch significant?
We moeten nu een zeer belangrijke vraag beantwoorden. Is deze waarde van α (0,037343) statistisch significant? Of kan zij voor alle praktische doeleinden als nul worden beschouwd?
Waarom is het zo belangrijk dit uit te zoeken? Bedenk dat als α nul is, de volgende vergelijking:
…reduceert tot Variantie = gemiddelde. Dit is de variantiefunctie van het Poisson-regressiemodel.
Als de waarde van α statistisch niet significant is, kan het Negative Binomial-regressiemodel de reeks opleidingsgegevens niet beter aanpassen dan een Poisson-regressiemodel.
Het OLSResults-object bevat de t-score van de regressiecoëfficiënt α. Laten we die afdrukken:
aux_olsr_results.tvalues
Dit drukt af:
BB_LAMBDA 4.814096
Uit een t-waardecalculator kunnen we afleiden dat de kritieke t-waarde bij een betrouwbaarheidsniveau van 99% (rechterstaart) en vrijheidsgraden=(172 waarnemingen) -( 6 regressievariabelen)- (1 dispersieparameter α)=165 2,34916 is. Dit is duidelijk minder dan de t-statistiek van α die 4,814096 bedroeg. We concluderen dat,
α=0,037343 statistisch significant is.
Dit voltooit STAP 2: De bepaling van α.
STAP 3: We voeren de in STAP 2 gevonden waarde van alpha in de statsmodels.genmod.families.family.NegativeBinomial-klasse in, en trainen het NB2-model op de reeks trainingsgegevens.
Dit is een bewerking in één stap in statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Zoals voorheen drukken we de trainingssamenvatting af:
print(nb2_training_results.summary())
Waardoor de volgende samenvatting wordt afgedrukt:
STAP 4: Laten we enkele voorspellingen doen met ons getrainde NB2-model.
Voorspellen is weer een procedure in één stap in statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Laten we de voorspellingen afdrukken:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Volgende zijn de eerste paar regels van de uitvoer:
Laten we ook de voorspelde tellingen tegenover de werkelijke tellingen voor de testgegevens uitzetten.
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 is de output:
Niet slecht! Het NB2-model lijkt de trend in de fietstellingen min of meer te volgen. En net als bij het Poisson regressiemodel wijken de voorspellingen in sommige gevallen ver af van de werkelijke waarden.