Comment faire de la régression binomiale négative en Python
Nous commencerons par importer tous les paquets requis.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Puis, créez un DataFrame pandas pour l’ensemble de données de comptage.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Nous ajouterons quelques variables de régression dérivées à la matrice X.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Nous n’utiliserons pas la variable Date comme régresseur puisqu’elle contient une valeur absolue de date mais nous n’avons pas besoin de faire quoi que ce soit de spécial pour abandonner Date car elle est déjà consommée comme indice du DataFrame de pandas. Elle ne sera donc pas disponible pour nous dans la matrice X.
Créons les ensembles de données de formation et de 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)))
STEP 1 : Nous allons maintenant configurer et ajuster le modèle de régression de Poisson sur l’ensemble de données de formation.
Configurer l’expression de régression en notation patsy. Nous disons à patsy que BB_COUNT est notre variable dépendante et qu’elle dépend des variables de régression : DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T et PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Constituer les matrices X et y pour les ensembles de données de formation et de test. patsy rend cela vraiment simple.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
En utilisant la classe statsmodels GLM, entraînez le modèle de régression de Poisson sur l’ensemble de données d’entraînement.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Cela termine l’entraînement du modèle de régression de Poisson. Pour voir le résultat de l’entraînement, vous pouvez imprimer le résumé de l’entraînement.
print(poisson_training_results.summary())
Cela imprime ce qui suit:
Notre intérêt réel réside dans le vecteur de taux ajustés λ produit par l’entraînement. Ce vecteur de taux est contenu dans le paramètre poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
La sortie suivante montre les quelques premières et les quelques dernières valeurs du vecteur λ ajusté :
Cela termine l’ÉTAPE1 : ajustement du modèle de régression de Poisson.
Étape 2 : Nous allons maintenant ajuster le modèle de régression MCO auxiliaire sur l’ensemble de données et utiliser le modèle ajusté pour obtenir la valeur de α.
Importer le paquet api.
import statsmodels.formula.api as smf
Ajouter le vecteur λ comme une nouvelle colonne appelée ‘BB_LAMBDA’ au cadre de données de l’ensemble de données d’entraînement. Rappelez-vous que les dimensions de λ sont (n x 1). Dans notre exemple, elle sera de (161 x 1). Rappelez-vous également que le vecteur λ est disponible dans poisson_training_results.mu :
df_train = poisson_training_results.mu
Puis, ajoutons une colonne dérivée appelée ‘AUX_OLS_DEP’ au Data Frame de pandas. Cette nouvelle colonne va stocker les valeurs de la variable dépendante de la régression MCO. C’est le côté gauche de l’équation de régression MCO ci-dessous :
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
Dans l’extrait de code ci-dessus, la partie en gras est le côté gauche de l’équation MCO auxiliaire ci-dessus.
Utilisons patsy pour former la spécification du modèle pour le OLSR. Nous voulons dire à patsy que AUX_OLS_DEP est la variable dépendante et qu’elle est expliquée par BB_LAMBDA (qui est le vecteur de taux λ). Le ‘-1’ à la fin de l’expression est la syntaxe de patsy pour dire : ne pas utiliser d’intercept de régression ; c’est-à-dire ajuster simplement une ligne droite passant par l’origine, comme suggéré par MM. Cameron et Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Nous sommes maintenant prêts à ajuster un modèle OLSR.
Configurer et ajuster le modèle OLSR :
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Imprimer les paramètres de régression :
print(aux_olsr_results.params)
Vous verrez s’imprimer le coefficient unique suivant correspondant à la variable unique de régression BB_LAMBDA. Ce coefficient est le α que nous recherchions:
BB_LAMBDA 0.037343
Est-ce que α est statistiquement significatif?
Nous devons maintenant répondre à une question très importante. Cette valeur de α (0,037343) est-elle statistiquement significative ? Ou peut-elle être considérée comme nulle à toutes fins utiles ?
Pourquoi est-il si important de le savoir ? Rappelons que si α est nul, alors l’équation suivante :
…se réduit à Variance = moyenne. Il s’agit de la fonction de variance du modèle de régression de Poisson.
Si la valeur de α n’est pas statistiquement significative, alors le modèle de régression binomiale négative ne peut pas faire un meilleur travail d’ajustement de l’ensemble de données d’entraînement qu’un modèle de régression de Poisson.
L’objet OLSResults contient le t-score du coefficient de régression α. Imprimons-le :
aux_olsr_results.tvalues
Ceci s’imprime :
BB_LAMBDA 4.814096
À partir d’un calculateur de valeur t, nous pouvons voir que la valeur t critique à un niveau de confiance de 99 % (queue droite), et des degrés de liberté=(172 observations) -( 6 variables de régression)- (1 paramètre de dispersion α)=165 est 2,34916. Ce chiffre est confortablement inférieur à la statistique t de α qui était de 4,814096. Nous concluons que,
α=0,037343 est statistiquement significatif.
Ceci termine l’ÉTAPE 2 : La détermination de α.
Étape 3 : Nous fournissons la valeur de alpha trouvée dans l’ÉTAPE 2 dans la classe statsmodels.genmod.families.family.NegativeBinomial, et entraînons le modèle NB2 sur l’ensemble de données d’entraînement.
C’est une opération en une étape dans statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Comme précédemment, nous allons imprimer le résumé de l’entraînement:
print(nb2_training_results.summary())
Qui imprime le résumé suivant :
STEP 4 : Faisons quelques prédictions en utilisant notre modèle NB2 entraîné.
La prédiction est à nouveau une procédure en une seule étape dans statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Imprimons les prédictions:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Voici les premières lignes de la sortie :
Traçons également les effectifs prédits en fonction des effectifs réels pour les données de 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()
Voici le résultat:
Pas trop mal ! Le modèle NB2 semble plus ou moins suivre la tendance des comptages de bicyclettes. Et tout comme pour les performances du modèle de régression de Poisson, dans certains cas, ses prédictions sont très éloignées des valeurs réelles.
.