Cómo hacer Regresión Binomial Negativa en Python
Comenzaremos importando todos los paquetes necesarios.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
A continuación, crearemos un DataFrame de pandas para el conjunto de datos de recuentos.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Añadiremos unas cuantas variables de regresión derivadas a la matriz X.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
No utilizaremos la variable Fecha como regresor ya que contiene un valor de fecha absoluto pero no necesitamos hacer nada especial para dejar caer Fecha ya que ya se consume como índice del DataFrame de pandas. Así que no estará disponible para nosotros en la matriz X.
Creemos los conjuntos de datos de entrenamiento y de prueba.
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)))
PASO 1: Ahora configuraremos y ajustaremos el modelo de regresión de Poisson en el conjunto de datos de entrenamiento.
Configuramos la expresión de regresión en notación patsy. Le estamos diciendo a patsy que BB_COUNT es nuestra variable dependiente y que depende de las variables de regresión: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T y PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Establezca las matrices X e y para los conjuntos de datos de entrenamiento y prueba. patsy hace que esto sea realmente sencillo.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Usando la clase GLM de statsmodels, entrena el modelo de regresión de Poisson en el conjunto de datos de entrenamiento.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Esto termina el entrenamiento del modelo de regresión de Poisson. Para ver el resultado del entrenamiento, puede imprimir el resumen del entrenamiento.
print(poisson_training_results.summary())
Esto imprime lo siguiente:
Nuestro verdadero interés radica en el vector de tasas ajustadas λ producido por el entrenamiento. Este vector de tasas está contenido en el parámetro poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
La siguiente salida muestra los primeros y últimos valores del vector λ ajustado:
Esto completa el PASO1: ajuste del modelo de regresión de Poisson.
PASO 2: Ahora ajustaremos el modelo de regresión OLS auxiliar en el conjunto de datos y utilizaremos el modelo ajustado para obtener el valor de α.
Importar el paquete api.
import statsmodels.formula.api as smf
Añadir el vector λ como una nueva columna llamada ‘BB_LAMBDA’ al Marco de datos del conjunto de datos de entrenamiento. Recuerde que las dimensiones de λ son (n x 1). En nuestro ejemplo será (161 x 1). Recordemos también que el vector λ está disponible en poisson_training_results.mu :
df_train = poisson_training_results.mu
A continuación, vamos a añadir una columna derivada llamada ‘AUX_OLS_DEP’ al Data Frame de pandas. Esta nueva columna almacenará los valores de la variable dependiente de la regresión OLS. Es el lado izquierdo de la ecuación de regresión OLS a continuación:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
En el fragmento de código anterior, la parte en negrita es el lado izquierdo de la ecuación OLSR auxiliar anterior.
Usemos patsy para formar la especificación del modelo para el OLSR. Queremos decirle a patsy que AUX_OLS_DEP es la variable dependiente y que se explica por BB_LAMBDA (que es el vector de tasas λ). El ‘-1’ al final de la expresión es la sintaxis de patsy para decir: no usar un intercepto de regresión; es decir, sólo ajustar una línea recta que pase por el origen, como sugieren los señores Cameron y Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Ahora estamos listos para ajustar un modelo OLSR.
Configure y ajuste el modelo OLSR:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Imprima los parámetros de regresión:
print(aux_olsr_results.params)
Verá que se imprime el siguiente coeficiente único correspondiente a la única variable de regresión BB_LAMBDA. Este coeficiente es el α que estábamos buscando:
BB_LAMBDA 0.037343
¿Es α estadísticamente significativo?
Ahora tenemos que responder a una pregunta muy importante. ¿Es este valor de α (0,037343) estadísticamente significativo? O puede considerarse nulo a efectos prácticos?
¿Por qué es tan importante averiguarlo? Recordemos que si α es cero, entonces la siguiente ecuación:
…se reduce a Varianza = media. Esta es la función de varianza del modelo de regresión de Poisson.
Si el valor de α no es estadísticamente significativo, entonces el modelo de regresión binomial negativa no puede hacer un mejor trabajo de ajuste del conjunto de datos de entrenamiento que un modelo de regresión de Poisson.
El objeto OLSResults contiene la puntuación t del coeficiente de regresión α. Imprimámoslo:
aux_olsr_results.tvalues
Esto imprime:
BB_LAMBDA 4.814096
Desde una calculadora de valores t, podemos ver que el valor t crítico a un nivel de confianza del 99% (cola derecha), y grados de libertad=(172 observaciones) -( 6 variables de regresión)- (1 parámetro de dispersión α)=165 es 2,34916. Esto es cómodamente menos que el estadístico t de α, que era de 4,814096. Concluimos que,
α=0,037343 es estadísticamente significativo.
Esto completa el PASO 2: La determinación de α.
PASO 3: Suministramos el valor de alfa encontrado en el PASO 2 en la clase statsmodels.genmod.families.family.NegativeBinomial, y entrenamos el modelo NB2 en el conjunto de datos de entrenamiento.
Esta es una operación de un solo paso en statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Como antes, imprimiremos el resumen de entrenamiento:
print(nb2_training_results.summary())
Que imprime el siguiente resumen:
PASO 4: Hagamos algunas predicciones utilizando nuestro modelo NB2 entrenado.
La predicción es de nuevo un procedimiento de un solo paso en statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Imprimamos las predicciones:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Las siguientes son las primeras líneas de la salida:
También vamos a trazar los recuentos predichos frente a los recuentos reales para los datos de prueba.
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()
Aquí está el resultado:
¡No está nada mal! El modelo NB2 parece seguir más o menos la tendencia de los recuentos de bicicletas. Y al igual que con el rendimiento del modelo de regresión de Poisson, en algunos casos sus predicciones están muy lejos de los valores reales.