Cum se face regresia binomială negativă în Python
Vom începe prin a importa toate pachetele necesare.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
În continuare, creați un pandas DataFrame pentru setul de date cu numărători.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Vom adăuga câteva variabile de regresie derivate la matricea X.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Nu vom folosi variabila Date ca regresor, deoarece conține o valoare absolută a datei, dar nu trebuie să facem nimic special pentru a renunța la Date, deoarece aceasta este deja consumată ca indice al pandas DataFrame. Așadar, nu ne va fi disponibilă în matricea X.
Să creăm seturile de date de antrenament și de testare.
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)))
PASUL 1: Acum vom configura și potrivi modelul de regresie Poisson pe setul de date de antrenament.
Configurați expresia de regresie în notația patsy. Îi spunem lui patsy că BB_COUNT este variabila noastră dependentă și că aceasta depinde de variabilele de regresie: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T și PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Configurați matricile X și y pentru seturile de date de instruire și de testare. patsy face acest lucru foarte simplu.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Utilizând clasa statsmodels GLM, antrenați modelul de regresie Poisson pe setul de date de antrenament.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Aceasta încheie antrenarea modelului de regresie Poisson. Pentru a vedea rezultatul antrenamentului, puteți tipări rezumatul antrenamentului.
print(poisson_training_results.summary())
Aceasta tipărește următoarele:
Interesul nostru real constă în vectorul de rate ajustate λ produs de instruire. Acest vector de rate este conținut în parametrul poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Vehiculul de ieșire de mai jos prezintă primele și ultimele valori ale vectorului λ ajustat:
Aceasta finalizează ETAPA 1: ajustarea modelului de regresie Poisson.
PASUL 2: Acum vom ajusta modelul auxiliar de regresie OLS pe setul de date și vom folosi modelul ajustat pentru a obține valoarea lui α.
Importați pachetul api.
import statsmodels.formula.api as smf
Adaugați vectorul λ ca o nouă coloană numită „BB_LAMBDA” în cadrul de date al setului de date de instruire. Reamintiți-vă că dimensiunile lui λ sunt (n x 1). În exemplul nostru, aceasta va fi (161 x 1). De asemenea, rețineți că vectorul λ este disponibil în poisson_training_results.mu :
df_train = poisson_training_results.mu
În continuare, să adăugăm o coloană derivată numită „AUX_OLS_DEP” la cadrul de date pandas. Această nouă coloană va stoca valorile variabilei dependente a regresiei OLS. Este partea stângă a ecuației de regresie OLS de mai jos:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
În fragmentul de cod de mai sus, partea în bold este partea stângă a ecuației aux OLSR de mai sus.
Să folosim patsy pentru a forma specificația modelului pentru OLSR. Dorim să îi spunem lui patsy că AUX_OLS_DEP este variabila dependentă și că aceasta este explicată de BB_LAMBDA (care este vectorul de rată λ). ‘-1’ de la sfârșitul expresiei este sintaxa lui patsy pentru a spune: nu folosiți o intercepție de regresie; adică potriviți doar o linie dreaptă care trece prin origine, așa cum au sugerat domnii Cameron și Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Suntem acum gata să potrivim un model OLSR.
Configurați și ajustați modelul OLSR:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Imprimați parametrii de regresie:
print(aux_olsr_results.params)
Vă veți vedea următorul coeficient unic tipărit corespunzător variabilei unice de regresie BB_LAMBDA. Acest coeficient este α pe care îl căutam:
BB_LAMBDA 0.037343
Este α semnificativ din punct de vedere statistic?
Acum trebuie să răspundem la o întrebare foarte importantă. Este această valoare a lui α (0,037343) semnificativă din punct de vedere statistic? Sau poate fi considerată ca fiind zero pentru toate scopurile practice?
De ce este atât de important să aflăm acest lucru? Reamintim că dacă α este zero, atunci următoarea ecuație:
…se reduce la Varianță = medie. Aceasta este funcția de varianță a modelului de regresie Poisson.
Dacă valoarea lui α nu este semnificativă din punct de vedere statistic, atunci modelul de regresie Binomial Negativ nu poate face o treabă mai bună de adaptare la setul de date de instruire decât un model de regresie Poisson.
Obiectul OLSResults conține scorul t al coeficientului de regresie α. Să o tipărim:
aux_olsr_results.tvalues
Aceasta se tipărește:
BB_LAMBDA 4.814096
Dintr-un calculator de valori t, putem vedea că valoarea t critică la un nivel de încredere de 99% (coadă dreaptă), și grade de libertate=(172 observații) -( 6 variabile de regresie)- (1 parametru de dispersie α)=165 este 2,34916. Aceasta este confortabil mai mică decât statistica t a lui α, care a fost de 4,814096. Concluzionăm că,
α=0,037343 este semnificativ din punct de vedere statistic.
Aceasta încheie ETAPA 2: Determinarea lui α.
ETAPA 3: Furnizăm valoarea lui alfa găsită în ETAPA 2 în clasa statsmodels.genmod.families.family.NegativeBinomial și antrenăm modelul NB2 pe setul de date de antrenament.
Aceasta este o operație într-o singură etapă în statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Ca și înainte, vom tipări rezumatul antrenamentului:
print(nb2_training_results.summary())
Care tipărește următorul rezumat:
ETAPA 4: Să facem câteva predicții folosind modelul nostru NB2 antrenat.
Predicția este din nou o procedură într-un singur pas în statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Să tipărim predicțiile:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
În cele ce urmează sunt primele câteva linii de ieșire:
Să reprezentăm, de asemenea, cifrele prezise față de cifrele reale pentru datele 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()
Iată rezultatul:
Nu prea rău! Modelul NB2 pare să urmărească mai mult sau mai puțin tendința în ceea ce privește numărul de bicicliști. Și, la fel ca în cazul performanței modelului de regresie Poisson, în unele cazuri predicțiile sale sunt departe de valorile reale.
.