Come fare la Regressione Binomiale Negativa in Python
Iniziamo importando tutti i pacchetti necessari.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
In seguito, creiamo un DataFrame pandas per l’insieme di dati dei conteggi.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Aggiungiamo alcune variabili di regressione derivate alla matrice X.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Non useremo la variabile Date come regressore poiché contiene un valore assoluto di data ma non abbiamo bisogno di fare nulla di speciale per eliminare Date poiché è già consumata come indice del DataFrame di pandas. Quindi non sarà disponibile per noi nella matrice X.
Creiamo i set di dati di allenamento e di test.
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)))
Passo 1: Ora configureremo e adatteremo il modello di regressione Poisson sul set di dati di allenamento.
Impostare l’espressione di regressione in notazione patsy. Stiamo dicendo a patsy che BB_COUNT è la nostra variabile dipendente e dipende dalle variabili di regressione: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T e PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Impostare le matrici X e y per i set di dati di allenamento e di test. patsy rende questo molto semplice.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Utilizzando la classe GLM di statsmodels, addestrate il modello di regressione di Poisson sul set di dati di allenamento.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Questo finisce l’addestramento del modello di regressione di Poisson. Per vedere il risultato dell’addestramento, è possibile stampare il riepilogo dell’addestramento.
print(poisson_training_results.summary())
Questo stampa quanto segue:
Il nostro vero interesse è il vettore dei tassi montati λ prodotto dall’addestramento. Questo vettore di tassi è contenuto nel parametro poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Il seguente output mostra i primi e gli ultimi valori del vettore λ montato:
Questo completa lo STEP1: adattare il modello di regressione Poisson.
STEP 2: Ora adatteremo il modello ausiliario di regressione OLS sul set di dati e useremo il modello montato per ottenere il valore di α.
Importiamo il pacchetto api.
import statsmodels.formula.api as smf
Aggiungiamo il vettore λ come una nuova colonna chiamata ‘BB_LAMBDA’ al Data Frame del set di dati di allenamento. Ricordate che le dimensioni di λ sono (n x 1). Nel nostro esempio sarà (161 x 1). Ricordiamo anche che il vettore λ è disponibile in poisson_training_results.mu :
df_train = poisson_training_results.mu
In seguito, aggiungiamo una colonna derivata chiamata ‘AUX_OLS_DEP’ al Data Frame di pandas. Questa nuova colonna memorizzerà i valori della variabile dipendente della regressione OLS. È il lato sinistro dell’equazione di regressione OLS qui sotto:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
Nello snippet di codice qui sopra, la parte in grassetto è il lato sinistro dell’equazione aux OLSR qui sopra.
Utilizziamo patsy per formare la specifica del modello per l’OLSR. Vogliamo dire a patsy che AUX_OLS_DEP è la variabile dipendente ed è spiegata da BB_LAMBDA (che è il vettore dei tassi λ). Il ‘-1’ alla fine dell’espressione è la sintassi di patsy per dire: non usare un’intercetta di regressione; cioè basta adattare una linea retta passante per l’origine, come suggerito dai signori Cameron e Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Siamo ora pronti per adattare un modello OLSR.
Configura e adatta il modello OLSR:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Stampa i parametri di regressione:
print(aux_olsr_results.params)
Vedrai il seguente coefficiente singolo stampato corrispondente alla singola variabile di regressione BB_LAMBDA. Questo coefficiente è l’α che stavamo cercando:
BB_LAMBDA 0.037343
È α statisticamente significativo?
Ora dobbiamo rispondere a una domanda molto importante. Questo valore di α (0,037343) è statisticamente significativo? O può essere considerato zero a tutti gli effetti?
Perché è così importante scoprirlo? Ricordiamo che se α è zero, allora la seguente equazione:
…si riduce a Varianza = media. Questa è la funzione di varianza del modello di regressione Poisson.
Se il valore di α non è statisticamente significativo, allora il modello di regressione binomiale negativa non può fare un lavoro migliore di un modello di regressione Poisson per adattare il set di dati di allenamento.
L’oggetto OLSResults contiene il t-score del coefficiente di regressione α. Stampiamolo:
aux_olsr_results.tvalues
Questo stampa:
BB_LAMBDA 4.814096
Da un calcolatore di valori t, possiamo vedere che il valore t critico ad un livello di confidenza del 99% (coda destra), e gradi di libertà=(172 osservazioni) -( 6 variabili di regressione)- (1 parametro di dispersione α)=165 è 2.34916. Questo è comodamente inferiore alla statistica t di α che era 4,814096. Concludiamo che,
α=0,037343 è statisticamente significativo.
Questo completa lo STEP 2: La determinazione di α.
STEP 3: Forniamo il valore di alfa trovato nello STEP 2 nella classe statsmodels.genmod.families.family.NegativeBinomial, e addestriamo il modello NB2 sul set di dati di allenamento.
Questa è un’operazione di un solo passo in statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Come prima, stamperemo il riassunto dell’allenamento:
print(nb2_training_results.summary())
Che stampa il seguente riassunto:
PASSO 4: Facciamo qualche previsione usando il nostro modello NB2 allenato.
La previsione è di nuovo una procedura a passo singolo in statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Stampiamo le previsioni:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Seguono le prime righe dell’output:
Tracciamo anche i conteggi previsti rispetto a quelli reali per i dati di test.
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()
Ecco il risultato:
Non male! Il modello NB2 sembra seguire più o meno la tendenza dei conteggi delle biciclette. E proprio come per le prestazioni del modello di regressione di Poisson, in alcuni casi le sue previsioni sono molto lontane dai valori reali.