使用 Python 进行泊松回归的教程
在本文中，我们将讨论以下主题：
-
基于计数的数据的特点：计数数据集是世界上最常见的数据。我们将探讨基于计数的数据的不同之处。 -
预测计数的回归模型：我们将详细介绍泊松回归模型。负二项（NB）回归模型是另一种常用的计数数据模型。我将在以后的文章中介绍。 -
关于泊松回归的 Python 教程：我将一步一步地教你如何使用statsmodels的GLM 类在 Python 中创建泊松回归模型，以及如何在真实世界的数据集上对其进行训练。
1. 在我们开始之前，有几个要点...
-
关于泊松回归的Python教程，请向下滚动**到本文的最后几节。 -
本文使用的现实世界自行车骑行者计数数据集是在此。 -
有关随机变量、泊松过程的入门知识，以及模拟泊松过程的 Python 程序，请点击此处： [[泊松回归]]
2. 什么是计数数据？
基于计数的数据包含以一定速率发生的事件。发生率可能会随着时间或观察结果的不同而变化。下面是一些基于计数的数据示例：
-
每小时通过十字路口的车辆数 -
每月到医生办公室就诊的人数 -
每月发现的类地行星数量
计数数据集具有以下特征：
-
整数数据：数据由非负整数组成： [0......∞] 普通最小二乘回归等回归技术可能不适合对此类数据建模，因为 OLSR 最适用于实数 -
偏斜分布：数据可能包含大量数据点，但只有几个值，因此频率分布相当偏斜。参见上面的直方图。 -
稀疏性：数据可能反映了罕见事件的发生，如伽马射线暴，从而使数据变得稀疏。 -
事件发生率：为了建立模型，可以假设有一定的事件发生率 来驱动此类数据的生成。随着时间的推移，事件发生率可能会发生偏移。
3. 现实世界的计数数据集
下表包含自行车骑行者通过纽约市各座桥梁的计数。这些计数是在 2017 年 4 月 1 日至 2017 年 10 月 31 日期间每天测量的。
下面是桥上自行车骑行者数量的时序图：
4. 计数回归模型
泊松回归模型和负二项回归模型是建立计数回归模型的两种常用技术。其他还有 有序 Logit、有序 Probit 和 非线性最小二乘法 模型。
5. 回归策略
好的做法是从泊松回归模型开始，并将其作为更复杂或限制较少的模型的 “对照”。在他们的著作 Regression Analysis of Count Data中，Cameron 和 Trivedi 这样描述：
“合理的做法是同时估计泊松模型和负二项模型”。
在本文中，我们将使用泊松回归模型对桥上观察到的骑自行车人数进行回归。
6. 介绍泊松模型
泊松分布具有以下概率质量函数：
泊松分布的期望值（平均值）为 λ.，因此，在没有其他信息的情况下，我们应该期望在任何单位时间间隔（如 1 小时、1 天等）内看到 λ 事件。对于任何时间间隔 t，人们都会期望看到 λt 事件。
7. 适用于 λ 的泊松回归模型
如果事件发生率 λ 是常数，我们就可以简单地使用修正的平均值模型来预测未来的事件计数。在这种情况下，我们可以将所有预测的计数值设为这个常数 λ。
下图说明了恒定 λ 的情况：
下面的 Python 代码使用 λ=5 的泊松过程生成蓝点（过去时间步骤中的实际计数）。黑色点（预测值）全部设置为相同的值 5。
8. 一个 非常数 λ 的泊松回归模型
现在我们进入有趣的部分。让我们来研究一种更常见的情况，即λ会随着观察结果的变化而变化。在这种情况下，我们假设λ的值受一个变量的向量影响，也称为预测变量、回归变量或回归变量。我们称这个回归变量矩阵为 。
回归模型的工作是将观测到的计数 与回归值矩阵 拟合。
在自行车骑行者计数数据集中，回归变量包括 日期 、星期、高温、低温 和 降水。我们还可以引入额外的回归变量，如从 日期 派生出来的 月份 和 每月的哪一天 ，我们也可以放弃现有的回归变量，如 日期 。
将 与 拟合是通过固定回归系数向量 的值来实现的。
在泊松回归模型中，事件计数 被假定为泊松分布，这意味着观测到 的概率是事件率向量 的函数。
泊松回归模型的工作是通过一个链接函数将观测到的计数y与回归矩阵X拟合，该链接函数将率向量λ表示为以下各项的函数
回归系数 回归矩阵
下图说明了泊松回归模型的结构。
从左到右扫描： 泊松回归模型的结构
连接 和 的 是一个怎样的链接函数呢？事实证明，下面的指数链接函数非常有效：
泊松回归模型的指数链接函数
即使回归系数 或回归系数 为负值，该链接函数也能保持 非负值。这是对基于计数的数据的要求。
一般来说，我们有:
泊松回归模型的正式说明
基于计数数据的泊松回归模型的完整说明如下：
对于数据集中的 ith 个观测值（以 表示，与回归变量 的行相对应），观测到计数 的概率是泊松分布的，如下所示：
给定 （根据泊松 PMF 公式），观察到计数 的概率
其中，第 个样本的平均率 由前面所示的指数链接函数给出。我们在此将其复制：
泊松回归模型的指数联系函数
一旦模型在数据集上得到充分训练，回归系数 已知，模型就可以进行预测了。要预测与已观察到的输入行回归系数 相对应的事件计数 ，可以使用下面的公式：
泊松回归模型的预测方程
所有这一切都取决于我们能否成功地训练模型，使回归系数向量 为已知。
让我们来看看如何进行训练。
训练泊松回归模型
训练泊松回归模型需要找到回归系数 的值，使观测到的计数向量 最有可能出现。
确定系数 的技术称为 Maximum Likelihood Estimation (MLE)。
让我们来了解一下 MLE 方法。
了解最大似然估计（MLE）
我将使用自行车骑行者计数数据集来说明MLE技术。请看数据集的前几行：
大桥上自行车骑行者的日计数
我们的假设是，红框中显示的自行车骑行者人数来自泊松过程。因此，我们可以说它们出现的概率是由泊松 PMF 给出的。以下是前 4 次出现的概率：
根据相应的回归向量，观察到前几次出现的骑自行车人数的概率
同样，我们也可以计算出在训练集中观察到的所有 n 计数的概率。
请注意，在上述公式中， 是使用链接函数计算的，如下所示：
与最初几天的计数相对应的事件发生率
其中， 是回归矩阵的前 4 行。
训练集中整个 n 个计数 的出现概率是各个计数的联合出现概率。
计数 y 是泊松分布， 是独立随机变量，相应地给定 。因此， 的联合概率可以表示为单个概率的简单相乘。下面是整个训练集的联合概率：
似然函数 L(β)表示为联合概率质量函数。
让我们回忆一下， 是通过回归系数 与回归向量 联系在一起的。
什么值的 会使给定的观察计数集 最有可能？就是上式所示的联合概率达到最大值 值。换句话说，它是联合概率函数随 变化的速率为 0 的 值。
微分联合概率方程的对数比微分原方程更容易。对数方程的解与 的最优值相同。
这个对数方程称为对数似然函数。对于泊松回归，对数似然函数由以下方程给出：
泊松回归模型的对数似然函数
将 ) 代替 后，对前面所示的联合概率函数两边取自然对数，即可得到上式。
如前所述，我们将这个对数概率方程与 β 进行微分，并将其设为零。通过这一运算，我们可以得到下面的方程：
的泊松 MLE 就是这个方程的解
对回归系数 求解这个方程，就会得到 的 最大似然估计值 (MLE) 。
要求解上述方程，需要使用迭代法，迭代加权最小二乘法（IRLS）。在实践中，我们不会手工求解这个方程。取而代之的是使用 Python statsmodels软件包等统计软件，在数据集上训练泊松回归模型时，这些软件会为您完成所有计算。
泊松回归的步骤摘要
总之，以下是对基于计数的数据集进行泊松回归的步骤：
-
首先，确保您的数据集包含计数。一种判断方法是，数据集只包含非负整数值，这些值代表某个时间间隔内某个事件的发生次数。在自行车骑行者计数数据集中，它是每天穿过布大桥的自行车骑行者计数。 -
找出（或猜测）会影响观察到的计数的回归变量。在自行车骑行者计数数据集中，回归变量是 _星期、最低温度、最高温度、降水量等。 -
挑选出一个训练数据集和一个测试数据集，前者用于训练回归模型，后者则放在一边。不要在测试数据上训练模型。 -
使用合适的统计软件（如 Python statsmodels 软件包）在训练数据集上配置和拟合泊松回归模型。 -
在测试数据集上运行该模型，以生成预测计数，从而测试模型的性能。将它们与测试数据集中的实际计数进行比较。 -
使用拟合优度来确定模型在训练数据集上的训练效果。
如何用 Python 进行泊松回归
让我们学以致用。Python statmodels软件包支持泊松回归。
我们的目标是为观测到的自行车骑行者人数 建立一个泊松回归模型。完成训练后，我们将使用训练好的模型来预测模型在训练过程中未来的大桥上每天的自行车骑行者人数。
大桥上自行车骑行者的每日计数
我们先导入所有需要的软件包。
import pandas as pd
from patsy import dmatrices
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
为计数数据集创建一个 pandas DataFrame
df = pd.read_csv('nyc_bb_bicyclist_counts.csv', header=0, infer_datetime_format=True, parse_dates=[0], index_col=[0])
我们将在 X 矩阵中添加一些衍生回归变量。
ds = df.index.to_series()
df['MONTH'] = ds.dt.month
df['DAY_OF_WEEK'] = ds.dt.dayofweek
df['DAY'] = ds.dt.day
我们不会使用 Date 变量作为回归变量，因为它包含一个绝对日期值，但我们不需要做任何特殊处理来删除 Date，因为它已经作为 pandas DataFrame 的索引被消耗掉了。因此，我们在 矩阵中将无法使用它。
让我们创建训练数据集和测试数据集。
mask = np.random.rand(len(df)) < 0.8
df_train = df[mask]
df_test = df[~mask]
print('Training data set length='+str(len(df_train)))
print('Testing data set length='+str(len(df_test)))
用 patsy 符号设置回归表达式。我们要告诉 patsy，BB/COUNT 是因变量，它取决于回归变量： 日（DAY）、周（DAY_OF_WEEK）、月（MONTH）、高（HIGH_T）、低（LOW_T）和前（PRECIP）。
expr = """BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP"""
为训练数据集和测试数据集设置 和 矩阵。
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
Using the statsmodels GLM class, train the Poisson regression model on the training data set.
poisson_training_results = sm.GLM(y_train, X_train,family=sm.families.Poisson()).fit()
打印培训摘要
print(poisson_training_results.summary())
打印输出如下:
泊松回归模型的训练总结
我们的模型效果如何？让我们对测试数据集进行一些预测。
poisson_predictions = poisson_training_results.get_prediction(X_test)
#.summary_frame() returns a pandas
DataFramepredictions_summary_frame = poisson_predictions.summary_frame()
print(predictions_summary_frame)
下面是输出结果的前几行：
-
mean：预测的 -
mean_se: 预测的 的标准误 -
mean_ci_lower: 置信区间下限 -
mean_ci_upper:置信区间上限 来自 poisson_predictions.summary_frame() 的前几行输出
让我们绘制测试数据的预测计数与实际计数的对比图。
predicted_counts=predictions_summary_frame['mean']
actual_counts = y_test['BB_COUNT']
fig = 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=[predicted, actual])plt.show()
输出结果如下：
大桥上自行车骑行者数量的预测与实际对比
该模型似乎或多或少地跟踪了实际人数的趋势，尽管**在许多情况下，其预测值与实际值相差甚远。
我们还可以绘制实际骑车人数与预测骑车人数的对比图。
plt.clf()
fig = plt.figure()
fig.suptitle('Scatter plot of Actual versus Predicted counts')
plt.scatter(x=predicted_counts, y=actual_counts, marker='.')
plt.xlabel('Predicted counts')
plt.ylabel('Actual counts')
plt.show()
情节是这样的
实际计数与预测计数的散点图
泊松回归模型的拟合优度
请记住，泊松分布的期望值（即均值）和方差都是 λ。 大多数实际数据都违反了这一相当严格的条件。
泊松回归模型失败的一个常见原因是数据不满足泊松分布所规定的 均值 = 方差 准则。
statsmodels GLMResult s 类上的 summary()方法显示了几个有用的拟合优度统计量，可以帮助您评估泊松回归模型是否能够成功拟合训练数据。让我们来看看它们的值：
泊松回归模型的训练摘要
所报告的偏离度和皮尔逊卡方值非常大。从这些值来看，几乎不可能拟合得很好。为了在一定置信度（例如 95%，p=0.05）下定量确定拟合优度，我们在 χ2 表中查找 p=0.05 和残差自由度=172 的值。(DF 残差 = 观察数 减 模型DF ）。我们将这个卡方值与观察到的统计量进行比较，在本例中，观察到的统计量是偏差或 GLMResults 中报告的皮尔逊卡方值。我们发现，在 p=0.05 和 DF 残差=172 时，标准卡方表 中的卡方值为 203，远小于报告的统计量 25091。因此，根据这项测试，尽管泊松回归模型对测试数据集的拟合效果 “还可以”，但对训练数据的拟合效果却相当差。
结论和下一步
对于基于计数的数据，泊松回归模型是一个有用的起点。然后，我们可以将其性能与其他流行的基于计数的模型进行比较，例如：
-
如果怀疑数据中包含过多的零，即比普通泊松模型所能解释的零多很多，则可以使用 零膨胀泊松模型。 -
负二项回归模型，它不对数据做出均值 = 方差的假设。 -
广义泊松回归 模型同样适用于过度分散或分散不足的数据集。