Como fazer Regressão Binomial Negativa em Python
Comecemos por importar todos os pacotes necessários.
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
Próximo, crie uma DataFrame pandas para o conjunto de dados de contagem.
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=, index_col=)
Adicionaremos algumas variáveis de regressão derivadas à matriz X.
ds = df.index.to_series()
df = ds.dt.month
df = ds.dt.dayofweek
df = ds.dt.day
Não usaremos a variável Data como um regressor já que ela contém um valor absoluto de data, mas não precisamos fazer nada de especial para deixar a Data como já está consumida como o índice da DataFrame pandas. Então ela não estará disponível para nós na matriz X.
Vamos criar os conjuntos de dados de treinamento e teste.
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: Vamos agora configurar e ajustar o modelo de regressão de Poisson no conjunto de dados de treinamento.
Configurar a expressão de regressão na notação de patsy. Estamos dizendo ao patsy que BB_COUNT é nossa variável dependente e depende das variáveis de regressão: DAY, DAY_OF_WEEK, MONTH, HIGH_T, LOW_T e PRECIP.
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
Configure as matrizes X e y para os conjuntos de dados de treino e teste. patsy torna isto muito simples.
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Usando a classe GLM dos modelos statsmodels, treine o modelo de regressão de Poisson no conjunto de dados de treinamento.
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
Acaba assim o treinamento do modelo de regressão de Poisson. Para ver o resultado do treinamento, você pode imprimir o resumo do treinamento.
print(poisson_training_results.summary())
Imprime o seguinte:
Nosso verdadeiro interesse está no vetor de taxas ajustadas λ produzido pelo treinamento. Este vector de taxas está contido no parâmetro poisson_training_results.mu.
print(poisson_training_results.mu)
print(len(poisson_training_results.mu))
A saída seguinte mostra os primeiros e últimos valores do vector ajustado λ vector:
Esta completa o PASSO 1: ajuste do modelo de regressão de Poisson.
STEP 2: Vamos agora ajustar o modelo de regressão auxiliar OLS no conjunto de dados e usar o modelo ajustado para obter o valor de α.
Importar o pacote api.
import statsmodels.formula.api as smf
Adicionar o vector λ como uma nova coluna chamada ‘BB_LAMBDA’ ao Quadro de Dados do conjunto de dados de treino. Recorde que as dimensões do λ são (n x 1). Em nosso exemplo será (161 x 1). Lembre-se também que o vetor λ está disponível em poisson_training_results.mu :
df_train = poisson_training_results.mu
Next, vamos adicionar uma coluna derivada chamada ‘AUX_OLS_DEP’ ao Quadro de Dados do pandas. Esta nova coluna irá armazenar os valores da variável dependente da regressão OLS. É o lado esquerdo da equação de regressão OLS abaixo:
df_train = df_train.apply(lambda x: ((x - x)**2 - x) / x, axis=1)
No trecho de código acima, a parte em negrito é o lado esquerdo da equação aux OLSR acima.
Vamos usar o bode expiatório para formar a especificação do modelo para o OLSR. Queremos dizer ao patsy que AUX_OLS_DEP é a variável dependente e é explicada por BB_LAMBDA (que é o vetor de taxa λ). O ‘-1’ no final da expressão é a sintaxe do patsy para dizer: não usar uma intercepção de regressão; isto é, apenas encaixar uma linha reta passando pela origem, como sugerido pelos senhores Cameron e Trivedi.
ols_expr = """AUX_OLS_DEP ~ BB_LAMBDA - 1"""
Estamos agora prontos para encaixar um modelo OLSR.
Configurar e encaixar o modelo OLSR:
aux_olsr_results = smf.ols(ols_expr, df_train).fit()
Imprimir os parâmetros de regressão:
print(aux_olsr_results.params)
Verá o seguinte coeficiente único a ser impresso correspondendo à variável de regressão única BB_LAMBDA. Este coeficiente é o α que estávamos procurando:
BB_LAMBDA 0.037343
Is α estatisticamente significante?
Agora precisamos responder a uma pergunta muito importante. Este valor de α (0,037343) é estatisticamente significante? Ou pode ser considerado zero para todos os fins práticos?
Por que é tão importante descobrir isto? Lembre-se que se α é zero, então a seguinte equação:
…reduz a Variância = média. Esta é a função de variância do modelo de regressão de Poisson.
Se o valor de α não for estatisticamente significativo, então o modelo de regressão Binomial Negativo não pode fazer um melhor ajuste do conjunto de dados de treinamento do que um modelo de regressão de Poisson.
O objeto OLSResults contém a pontuação t do coeficiente de regressão α. Vamos imprimir:
aux_olsr_results.tvalues
Esta impressão:
BB_LAMBDA 4.814096
De uma calculadora de valor t, podemos ver que o valor t crítico a um nível de confiança de 99% (cauda direita), e graus de liberdade=(172 observações) -( 6 variáveis de regressão)- (1 parâmetro de dispersão α)=165 é 2,34916. Isto é confortavelmente menor que a estatística t de α que era 4,814096. Concluímos que,
α=0,037343 é estatisticamente significativo.
>
Este conclui o PASSO 2: A determinação de α.
PASSO 3: Fornecemos o valor de alfa encontrado no PASSO 2 para a classe statsmodels.genmod.families.family.NegativeBinomial, e treinamos o modelo NB2 no conjunto de dados de treinamento.
Esta é uma operação de um passo em statsmodels:
nb2_training_results = sm.GLM(y_train, X_train,family=sm.families.NegativeBinomial(alpha=aux_olsr_results.params)).fit()
Como antes, vamos imprimir o resumo do treinamento:
print(nb2_training_results.summary())
Que imprime o seguinte resumo:
STEP 4: Vamos fazer algumas previsões usando o nosso modelo NB2 treinado.
Previsão é novamente um procedimento de um único passo em statsmodels:
nb2_predictions = nb2_training_results.get_prediction(X_test)
Imprimamos as previsões:
predictions_summary_frame = nb2_predictions.summary_frame()
print(predictions_summary_frame)
Seguir são as primeiras linhas da saída:
Vamos também plotar as contagens previstas versus as contagens reais para os dados do teste.
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()
Aqui está a saída:
Não é muito mau! O modelo NB2 parece estar mais ou menos seguindo a tendência na contagem de bicicletas. E tal como no modelo de regressão de Poisson, em alguns casos as suas previsões estão muito longe dos valores reais.