How to do Negative Binomial Regression in Python
Zaczniemy od zaimportowania wszystkich wymaganych pakietów.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Następnie utwórz pandas DataFrame dla zbioru danych counts.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Do macierzy X dodamy kilka pochodnych zmiennych regresji.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Nie będziemy używać zmiennej Date jako regresora, ponieważ zawiera ona bezwzględną wartość daty, ale nie musimy robić nic specjalnego, aby upuścić Date, ponieważ jest już zużywana jako indeks pandas DataFrame. Nie będzie więc dla nas dostępna w macierzy X.
Utwórzmy zestawy danych treningowych i testowych.
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)))
KROK 1: Skonfigurujemy teraz i dopasujemy model regresji Poissona na zestawie danych treningowych.
Skonfigurujmy wyrażenie regresji w notacji patsy. Mówimy patsy, że BB_COUNT jest naszą zmienną zależną i zależy od zmiennych regresji: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T i PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Ustawiamy macierze X i y dla zbiorów danych treningowych i testowych. patsy sprawia, że jest to naprawdę proste.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Używając klasy statsmodels GLM, wytrenuj model regresji Poissona na zestawie danych treningowych.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
To kończy trening modelu regresji Poissona. Aby zobaczyć wynik treningu, można wydrukować podsumowanie treningu.
print(poisson_training_results.summary())
Wyświetla to następujące wyniki:
Naszym prawdziwym zainteresowaniem cieszy się wektor dopasowanych stawek λ otrzymany w wyniku szkolenia. Wektor ten jest zawarty w parametrze poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
Następujące dane wyjściowe pokazują kilka pierwszych i kilka ostatnich wartości dopasowanego wektora λ:
To kończy KROK 1: dopasowanie modelu regresji Poissona.
KROK 2: Dopasujemy teraz pomocniczy model regresji OLS na zbiorze danych i użyjemy dopasowanego modelu do uzyskania wartości α.
Importujemy pakiet api.
import statsmodels.formula.api as smf
Dodajemy wektor λ jako nową kolumnę o nazwie 'BB_LAMBDA’ do ramki danych zbioru danych treningowych. Przypomnijmy, że wymiar λ to (n x 1). W naszym przykładzie będzie to (161 x 1). Przypomnijmy również, że wektor λ jest dostępny w poisson_training_results.mu :
df_train = poisson_training_results.mu
Następnie dodajmy pochodną kolumnę o nazwie 'AUX_OLS_DEP’ do ramki danych pandas Data Frame. Ta nowa kolumna będzie przechowywać wartości zmiennej zależnej regresji OLS. Jest to lewa strona poniższego równania regresji OLS:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
W powyższym wycinku kodu część pogrubiona jest lewą stroną powyższego równania aux OLSR.
Użyjmy patsy do utworzenia specyfikacji modelu dla OLSR. Chcemy powiedzieć patsy, że AUX_OLS_DEP jest zmienną zależną i jest ona objaśniana przez BB_LAMBDA (który jest wektorem stopy λ). ’-1′ na końcu wyrażenia to składnia patsy mówiąca: nie używaj punktu przecięcia regresji; tzn. po prostu dopasuj linię prostą przechodzącą przez początek, jak sugerują panowie Cameron i Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Jesteśmy teraz gotowi do dopasowania modelu OLSR.
Konfiguracja i dopasowanie modelu OLSR:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Wydruk parametrów regresji:
print(aux_olsr_results.params)
Zobaczysz następujący pojedynczy współczynnik odpowiadający pojedynczej zmiennej regresji BB_LAMBDA. Ten współczynnik to α, którego szukaliśmy:
BB_LAMBDA 0.037343
Is α statistically significant?
Teraz musimy odpowiedzieć na bardzo ważne pytanie. Czy ta wartość α (0,037343) jest statystycznie istotna? Czy też można ją uznać za zero dla wszystkich praktycznych celów?
Dlaczego tak ważne jest, aby się tego dowiedzieć? Przypomnijmy, że jeśli α wynosi zero, to następujące równanie:
…sprowadza się do Wariancja = średnia. Jest to funkcja wariancji modelu regresji Poissona.
Jeśli wartość α jest statystycznie nieistotna, to model regresji ujemnej dwumianowej nie może lepiej dopasować zestawu danych treningowych niż model regresji Poissona.
Obiekt OLSResults zawiera t-score współczynnika regresji α. Wypiszmy go:
aux_olsr_results.tvalues
Wyświetla to:
BB_LAMBDA 4.814096
Z kalkulatora wartości t widzimy, że krytyczna wartość t przy 99% poziomie ufności (prawe ogony) i stopniach swobody=(172 obserwacje) -( 6 zmiennych regresji)- (1 parametr dyspersji α)=165 wynosi 2,34916. Jest to wartość wyraźnie mniejsza od t-statystyki α, która wyniosła 4,814096. Wnioskujemy, że,
α=0,037343 jest istotnie statystycznie.
To kończy KROK 2: Wyznaczenie α.
KROK 3: Dostarczamy wartość alfa znalezioną w KROKU 2 do klasy statsmodels.genmod.families.family.NegativeBinomial i trenujemy model NB2 na zbiorze danych treningowych.
Jest to operacja jednoetapowa w statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Tak jak poprzednio, wydrukujemy podsumowanie treningu:
print(nb2_training_results.summary())
Wydrukuje ono następujące podsumowanie:
KROK 4: Zróbmy kilka predykcji przy użyciu naszego wytrenowanego modelu NB2.
Predykcja jest ponownie procedurą jednoetapową w statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Wypisujmy predykcje:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Następnie kilka pierwszych linii danych wyjściowych:
Wykreślmy również przewidywane zliczenia w stosunku do rzeczywistych zliczeń dla danych testowych.
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()
Oto dane wyjściowe:
Nieźle! Model NB2 wydaje się w mniejszym lub większym stopniu podążać za trendem w liczbach rowerzystów. I podobnie jak w przypadku modelu regresji Poissona, w niektórych przypadkach jego przewidywania znacznie odbiegają od rzeczywistych wartości.
.