用于时间序列数据的泊松回归模型

泊松和类泊松回归模型常用于基于计数的数据集，即包含整数计数的数据。例如，每小时走进医院急诊室的人数就是一个这样的数据集。

基于普通最小二乘回归的线性模型或非线性模型（例如基于基于神经网络的回归技术的线性模型）不适用于此类数据集，因为它们可以预测负值。

如果数据集是计数的时间序列，则会产生额外的建模复杂性，因为时间序列数据通常是自相关的。以前的计数会影响将来计数的值。如果回归模型不能充分捕获这些相关性中包含的“信息”，“未解释的”信息将以自相关误差的形式泄漏到模型的残差中。在这种情况下，模型的拟合优度会很差。

解决这个问题的一般补救办法如下:

1. 在拟合回归模型之前，检查时间序列是否具有季节性，如果具有，则进行季节性调整。这样做，就算有季节性的自相关性，也可以解释得通。
2. 对所有t进行时间序列的第一次差分，即y_t - y_(t-1)，并对差分时间序列进行白噪声测试。如果差分时间序列是白噪声，则原始时间序列是随机游走。在这种情况下，不需要进一步建模。
3. 在季节性调整后的时间序列上拟合基于Poisson（或相关）计数的回归模型，但包括因变量y的滞后副本作为回归变量。

在本文中，我们将解释如何使用方法(3)在计数的时间序列上拟合泊松或类泊松模型。

MANUFACTURING STRIKES数据集

为了说明模型拟合过程，我们将使用以下在回归建模文献中广泛使用的开源数据集：

该数据是一个月度时间序列，显示了从1968年到1976年，美国制造业活动(偏离趋势线)与美国制造业合同罢工数量之间的关系。

这个数据可以使用statsmodels数据集包来获取。

因变量y是罢工数量。

我们将从导入所有必需的程序包开始：

 import statsmodels.api as sm
 import statsmodels.discrete.discrete_model as dm
 import numpy as np
 from patsy import dmatrices
 import statsmodels.graphics.tsaplots as tsa
 from matplotlib import pyplot as plt

让我们使用statsmodels将数据集加载到内存中：

 strikes_dataset = sm.datasets.get_rdataset(dataname='StrikeNb', package='Ecdat')

打印出数据集：

 print(strikes_dataset.data)

我们看到以下输出：

我们将前92个数据点视为训练集，其余16个数据点视为测试数据集

 strikes_data = strikes_dataset.data.copy()
 strikes_data_train = strikes_data.query('time<=92')
 strikes_data_test = strikes_data.query('time>92').reset_index().drop('index', axis=1)

打印出因变量的统计信息：

 print('Mean='+str(np.mean(strikes_data_train['strikes'])) + ' Variance='+str(np.var(strikes_data_train['strikes'])))

我们得到以下输出：

 Mean=5.5 Variance=14.728260869565217

我们看到y过于分散，因此它违反了Poisson模型的均值=方差假设。为了解决过度分散问题，我们将拟合具有以下NB2方差函数的负二项式回归模型：

这是我们的回归表达式。罢工是因变量，输出是我们的解释变量。假定存在回归截距：

 expr = 'strikes ~ output

我们将使用Patsy来计算X和y矩阵。Patsy会自动向X添加回归截距列：

 y_train, X_train = dmatrices(expr, strikes_data_train, return_type='dataframe')
 print(y_train)
 print(X_train)
 
 y_test, X_test = dmatrices(expr, strikes_data_test, return_type='dataframe')
 print(y_test)
 print(X_test)

使用NB2方差函数构建和训练负二项式回归模型：

 nb2_model = dm.NegativeBinomial(endog=y_train, exog=X_train, loglike_method='nb2')
 nb2_model_results = nb2_model.fit(maxiter=100)
 print(nb2_model_results.summary())

我们得到了拟合模型摘要的以下输出：

如系数的p值（0.034和0.000）所示，输出和色散参数alpha在95％置信度上都具有统计学意义。

贴合度

伪R平方仅为0.9％，说明训练数据集的拟合质量非常差。

对数似然比检验的p值为0.03589，表明该模型在95％的置信水平下比仅截距模型（又称为零模型）做得更好，但在99％或更高的置信水平下表现不佳。

让我们看一下拟合模型的残差的自相关图：

我们可以看到残差误差在时间滞后1、2和3时是自相关的，这表明因变量罢工中存在自相关，因为NB2模型无法完全解释导致泄漏到模型残差中的原因。。

总体而言，此模型的拟合优度非常差。

建立自回归泊松模型

为了解决残差自相关的情况，我们将引入y的滞后副本，具体为y(t-1)、y(t-2)和y(t-3)作为输出变量的回归变量。但我们不直接引入y(t-k)作为回归变量，而是使用ln[y(t-k)]来解决y(t-k)系数为正时的“模型爆炸”问题。

但是使用ln()变换产生了一个问题即如何处理0的y_t对于这个对数是没有定义的。所以我们使用Cameron和Trivedi在他们的书《Regression Analysis of Count Data》(见第7.5节:自回归模型)中概述的以下技巧来解决这个问题:

我们将为每个感兴趣的时延定义一个新的指标变量d_t，如下所示:

• 当y_t = 0时:d_t=1，我们将y_t设为0。
• 当y_t > 0: d_t=0时，我们将y_t保持原样。

让我们对dataframe进行这些更改。创建一个新列strikes_adj，如果strikes < 1，则设置为1，否则设置为strikes:

 strikes_data['strikes_adj'] = np.maximum(1, strikes_data['strikes'])

定义一个函数来设置指标变量d_t的值，如上所述。

 def indicator_func(x):
     if x == 0:
         return 1
     else:
         return 0

并使用此函数创建一个新的指标变量列：

 strikes_data['d'] = strikes_data['strikes'].apply(indicator_func)

现在为strikes_adj和d创建滞后变量：

 strikes_data['ln_strikes_adj_lag1'] = strikes_data['strikes_adj'].shift(1)
 strikes_data['ln_strikes_adj_lag2'] = strikes_data['strikes_adj'].shift(2)
 strikes_data['ln_strikes_adj_lag3'] = strikes_data['strikes_adj'].shift(3)
 
 strikes_data['d_lag1'] = strikes_data['d'].shift(1)
 strikes_data['d_lag2'] = strikes_data['d'].shift(2)
 strikes_data['d_lag3'] = strikes_data['d'].shift(3)

删除所有带有空单元格的行：

 strikes_data = strikes_data.dropna()

最后，取一个自然日志ln_strikes_adj_lag1。回想一下，我们想添加滞后变量y（t_1），y（t_2）和y_（t_3）的自然对数。

让我们看看我们的数据现在是什么样子：

 print(strikes_data)

让我们再次将数据分为训练和测试数据集：

 strikes_data_train=strikes_data.query('time<=92')
 strikes_data_test=strikes_data.query('time>92').reset_index().drop('index', axis=1)

我们的回归表达式也需要更新以包含滞后变量：

 expr = 'strikes ~ output + ln_strikes_adj_lag1 + ln_strikes_adj_lag2 + ln_strikes_adj_lag3 + d_lag1 + d_lag2 + d_lag3'

使用Patsy可以得出y和X矩阵：

 y_train, X_train = dmatrices(expr, strikes_data_train, return_type='dataframe')
 print(y_train)
 print(X_train)

最后，我们将在（y_train，X_train）上构建并拟合回归模型。这次，我们将使用直接的Poisson回归模型：

 poisson_model = dm.Poisson(endog=y_train, exog=X_train)
 poisson_model_results = poisson_model.fit(maxiter=100)
 print(poisson_model_results.summary()

我们看到以下结果：

贴合度

首先要注意的是，通过Pseudo-R-squared测量的拟合优度比早期的NB2模型提高了0.9％至15.69％。这是一个很大的进步。这次，LLR测试的p值也很小，为1.295e-15。这意味着我们可以以接近100％的置信度说滞后变量Poisson模型要比截距模型好。回想一下，我们可以说只有95％的置信度。

让我们看一下这个滞后变量泊松模型的残差的自相关图：

 tsa.plot_acf(poisson_model_results.resid, alpha=0.05)
 plt.show()

我们看到以下图：

除了在LAG 13处存在非常轻微的显着相关性之外，残差与所有其他滞后的相关性都很好地位于规定的alpha范围内。

我们在Poisson模型的回归变量中添加滞后罢工副本的策略似乎已经解释了很多罢工变量中的自相关。

变量的意义

我们从训练注意总结滞后变量的泊松模型的系数,输出,ln_strikes_adj_lag1和ln_strikes_adj_lag2是重要的在95%置信水平,第三个滞后ln_strikes_adj_lag3是重要的系数只有大约75%置信水平表示的假定值为0.237。此外，三个滞后指标变量d_lag1、d_lag2和d_lag3在95%置信水平下均不具有统计学意义。

预测

让我们使用拟合的滞后变量Poisson模型来预测我们先前预留的测试数据集的罢工次数。我们不应该寄希望于预测的质量过高。请记住，尽管此模型比以前的NB2模型拟合得好得多，但伪R平方仍然仅为16％。

我们将使用Patsy来查看(y_test, X_test):

 y_test, X_test = dmatrices(expr, strikes_data_test, return_type='dataframe')
 print(y_test)
 print(X_test)

对X_test进行预测:

 poisson_predictions = poisson_model_results.predict(X_test

绘制预测值和实际值:

 predicted_counts=poisson_predictions
 actual_counts = y_test['strikes']
 fig = plt.figure()
 fig.suptitle('Predicted versus actual strike counts')
 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=[predicted, actual])
 plt.show()

我们得到如下图:

下一步

我们可以尝试通过以下修改来改善滞后变量模型的拟合优度:

除输出外，还将输出变量的前三个时间滞后量作为回归变量。

将输出变量和罢工变量的时滞值作为回归变量。

使用负二项模型(使用NB1或NB2方差函数)代替泊松模型，并将上述类型的滞后变量作为回归变量。

论文和相关连接

Cameron A. Colin, Trivedi Pravin K., Regression Analysis of Count Data, Econometric Society Monograph №30, Cambridge University Press, 1998. ISBN: 0521635675

McCullagh P., Nelder John A., Generalized Linear Models, 2nd Ed., CRC Press, 1989, ISBN 0412317605, 9780412317606

Kennan J., The duration of contract strikes in U.S. manufacturing, Journal of Econometrics, Volume 28, Issue 1, 1985, Pages 5–28, ISSN 0304–4076, https://doi.org/10.1016/0304-4076(85)90064-8. PDF download link

Cameron C. A., Trivedi P. K., Regression-based tests for overdispersion in the Poisson model, Journal of Econometrics, Volume 46, Issue 3, 1990, Pages 347–364, ISSN 0304–4076, https://doi.org/10.1016/0304-4076(90)90014-K.

作者：Sachin Date

原文地址：https://towardsdatascience.com/poisson-regression-models-for-time-series-data-sets-54114e68c46d

deephub翻译组