在Python中做广义泊松回归的步骤教程
常规的泊松回归模型通常是基于计数的数据集的首选模型。泊松回归模型的主要假设是：计数的方差与平均值相同，也就是说，数据是等距分散的。不幸的是，现实世界的数据很少是等距分散的，这就促使统计学家采用其他的计数模型，如：
-
负二项式（NB）回归模型 -
广义泊松回归模型
这两个模型都没有对数据集进行 方差等于均值的假设。这使得它们成为对基于计数的数据进行建模的更好的实际选择。在以前的文章中，我详细介绍了NB模型。在这篇文章中，我们将介绍通用泊松回归模型。
文章大纲
在这篇文章中，我们将涵盖以下主题：
-
什么是基于计数的数据集？ -
什么是等离散、欠离散、过离散？我们将介绍泊松模型对低分散和高分散数据集的限制。 -
介绍Consul的广义泊松回归(GP-1)模型和Famoye的受限广义泊松回归(GP-2)模型。 -
基于Python的教程，用于建立和训练GP-1和GP-2模型，并比较它们与标准泊松回归模型的性能。
什么是基于计数的数据？
基于计数的数据集是指因变量 代表某些事件的发生次数。下面是一些例子：
-
每个月网站的访问PV。 -
每年年来太阳黑子的数量。 -
每小时走进急诊室的人数。
在上述每个数据集中，因变量 代表观察到的计数，而解释变量矩阵 的选择，即被认为可以解释 y 的变量集，主要留给了统计建模者的判断。
等离散、欠离散、过离散和泊松模型的限制
人们可以首先尝试将因果计数变量建模为泊松过程。不幸的是，在许多现实世界的数据集中，泊松过程无法令人满意地解释观察到的计数的变化性。这主要是由于泊松回归模型对数据所作的 变量=平均值 的假设。换句话说，泊松模型不正确地假设计数是等距分散的。
泊松模型无法摆脱这一假设，因为该模型是基于因变量 是遵循泊松概率分布的随机变量这一假设，可以证明，泊松分布的随机变量的方差等于其平均值。
下面的公式表示泊松分布随机变量的概率分布函数(也称为Probability Mass F函数)。可以证明：
单位时间内发生的事件数。
在单位时间内发生λ个事件的情况下，看到 个事件的概率。
下面是一些看起来很有趣的概率分布图，由上述公式生成的不同的λ值：
不同的λ值下看到 个事件的概率
顺便说一下，在回归模型中，谈论无条件方差和无条件平均数是不常见的。相反，人们更喜欢把方差说成是[[条件概率]]
条件方差可以表示为Variance( )。同理，条件平均数也是如此Mean( )
如前所述，许多现实世界的数据集不是等散的，所以 的条件方差不等于 的条件均值。其中有些数据集的方差大于均值，这种现象被称为过度分散，而在其他数据集中，方差小于均值，又称欠分散。
一般来说：
当你的因变量 y 的方差大于你的模型所假设的方差时，那么从你的模型的角度来看，你的数据集是**过度分散的。
当你的因变量**y**的方差小于你的模型所假设的方差时，从你的模型的角度看，你的数据集是欠分散的。
过度（欠）分散的影响是，你的模型将不能充分地解释观察到的计数的变异性。相应地，其预测的质量也会很差。
对这个问题的一个常见的解决方法是假设方差是平均值的某个一般函数，即
f(.)的常用形式之一是如下：
一个常用的方差函数
其中 p=0,1,2…
α 被称为分散参数，它代表了 的额外变异性，由一些未知的变量集引入，导致y的总体变异性与你的回归模型所预期的不同。
请注意，当 α=0 时，方差=平均值 ，你会得到标准的泊松回归模型。
当 α>0 时，我们考虑两种有趣的情况，即当 p=1和p=2时 在这两种情况下，我们得到的是所谓的负二项式（NB）回归模型。
在NB回归模型中，我们假设观察到的计数 y 是一个泊松分布的随机变量，事件率为 λ ，λ 本身是一个伽马分布的随机变量。
负二项回归模型（通常被称为泊松-伽马混合物）被证明是许多现实世界中计数数据集的一个优秀模型。我在下面的文章中更详细地介绍了它：
[[负二项式回归]]
还有许多其他的方法来概括泊松回归模型，以增加其对过度分散和欠分散的数据的影响。在这篇文章中，我们将介绍两个这样的模型。
Consul的广义泊松回归模型
1989年，Prem C. Consul在他的书通用泊松分布：属性与应用中，提出了一种修改泊松分布的概率分布的方法，使其既能处理过度分散的数据，又能处理分散不足的数据。这个模型后来被称为GP-1（Generalized Poisson-1）模型。GP-1模型假定因变量 y 是一个具有以下概率分布的随机变量：
GP-1模型的概率分布和方差及平均函数
如果在上述方程中把色散参数 α 设为0，则GP-1的PMF、均值和方差基本上会减少到标准泊松分布的那些参数。
Famoye的受限广义泊松分布
1993年，Felix Famoye引入了他所说的受限广义泊松回归模型，作为一种扩展标准泊松模型的方法，以处理过度分散和不足分散的数据集。这个模型后来被称为GP-2（Generalized Poisson-2）模型。
GP-2模型假设因变量 y 是一个随机变量，其概率分布如下：
如前所述，当分散参数 α 被设置为0时，GP-2的PMF、平均数和方差函数基本上减少到标准泊松分布的那些。
分散参数α是用以下公式估计的：
N: 训练样本数量 k: 回归变量的数量 : 第i个观测值 : 对应第i个 泊松率 对应的预测值 p： 1 或者 2 GP1 或者GP2 模型
在 statsmodels 库中，α 是作为GP-1和GP-2模型拟合过程的估计的一部分。
在Python和Statsmodels中实现GP-1和GP-2
statsmodels库包含了通过statsmodels.discrete.discrete_model.GeneralizedPoisson类对GP-1和GP-2模型的实现。
在本节中，我们将展示如何使用GP-1和GP-2对以下现实世界的计数数据集进行建模。
一个真实世界的计数数据集
下表包含了骑车人在纽约市各座桥梁上行驶的车辆计数。时间范围在2017年4月1日至2017年10月31日。
我们将重点分析每天通过 布鲁克林大桥 的自行车数量。下面是在布鲁克林桥上看到的自行车人数的时间序列图：
我们的回归目标
我们的回归目标是预测在任何一天穿越布鲁克林大桥的骑自行车的人数。
我们的回归策略
我们将使用一组来自数据集的回归变量，即日、周（来自日期）、月（来自日期）、当日最高高温、当日最低温和降水，来 "解释"布鲁克林大桥上观察到的数量的差异。
回归矩阵和观察到的骑车人数量向量
我们的回归模型的训练算法将把观察到的计数 y 与回归矩阵 **X** 相匹配。
一旦模型训练完成，我们将在模型训练期间完全没有看到的测试数据集上测试其模型表现。
回归策略--一步一步来
-
我们将首先在这个数据集上训练标准泊松回归模型。标准泊松模型也将作为 "对照 "模型，用于测试广义泊松模型的功效。 -
如果数据集的方差与平均值大致相同，那么除了检查泊松模型的拟合度之外，就没有什么可做的了。 -
如果方差大于或小于泊松平均值，那么我们将在数据集上依次训练GP-1和GP-2模型，并将其性能与泊松模型进行比较。
我们将从导入所有需要的包：
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('泊松回归_自行车数量.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
让我们打印出我们的数据集的前几行，看看它看起来是什么样子：
print(df.head(10))
我们来创建训练和测试数据集：
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是我们的因变量 y ，它取决于回归变量 X ： day、DAY_OF_WEEK、month，HIGH_T、LOW_T和PRECIP。
expr = 'BB_COUNT ~ DAY + DAY_OF_WEEK + MONTH + HIGH_T + LOW_T + PRECIP'
让我们用Patsy来刻画训练和测试数据集的 X 和 y 矩阵：
y_train, X_train = dmatrices(expr, df_train, return_type='dataframe')
y_test, X_test = dmatrices(expr, df_test, return_type='dataframe')
使用statsmodels GLM类，在训练数据集上训练泊松回归模型。
poisson_training_results = sm.GLM(y_train, X_train, family=sm.families.Poisson()).fit()
打印训练模型总结：
print(poisson_training_results.summary())
这就打印出了以下内容：
我在输出中强调了重要的部分。
我们可以看到，所有的回归系数（在回归术语中称为 β 向量）在95%的置信度下都具有统计学意义，因为它们的p值<0.05。
最大对数似然值（-11872）将在以后用于比较该模型与GP-1和GP-2模型的模型表现。
让我们打印出数据集的方差和平均值：
print('variance='+str(df['BB_COUNT'].var()))
print('mean='+str(df['BB_COUNT'].mean()))
这就打印出了以下内容：
variance=730530.6601948135
mean=2680.042056074766
方差显然比平均值大得多。数据严重过度分散，泊松模型的主要假设并不成立。
因此，我们接下来将建立并训练GP-1和GP-2模型，看看它们的表现是否更好。
Statsmodels让你在3行代码中就能完成这个任务!
建立Consul的广义泊松回归模型，也就是GP-1：
gen_poisson_gp1 = sm.GeneralizedPoisson(y_train, X_train, p=1)
拟合(训练)模型：
gen_poisson_gp1_results = gen_poisson_gp1.fit()
打印结果：
print(gen_poisson_gp1_results.summary())
这将打印出以下内容（我在输出中强调了几个有趣的地方）：
请注意，除了DAY变量的系数，其他所有的回归系数都在95%的置信水平上具有统计学意义，即它们的P值>0.05。请记住，在泊松模型中，所有的回归系数都在95%的置信水平上具有统计学意义。
GP-1模型的拟合度
GP-1模型的最大似然估计值为-1350.6，大于无效模型的最大似然估计值-1475.9。空模型是一个简单的仅有截距的模型，即一条通过Y截距的水平线。但是，这种差异在统计上有意义吗？概率比（LR）测试的p值显示为3.12e-51，一个极其微小的数字。因此，是的，GP-1模型实际上确实比简单的仅有截距的模型做了更好的数据建模。
但GP-1是否比普通的泊松模型做得更好呢？
请记住，常规泊松模型的最大似然估计值是-11872。与GP-1的最大似然估计值-1350.6相比。显然，GP-1的拟合度更高（它的最大似然估计值大得多）。如果我们愿意，我们可以用这两个似然估计值做一个 LR 检验，然后用 Chi-squared 检验来评估其显著性。在这种情况下，这是没有必要的，因为GP-1的最大似然估计值要比常规泊松模型的大得多。
预测
让我们用GP-1来预测骑自行车的人的数量，使用模型在训练期间没有看到的测试数据集：
gen_poisson_gp1_predictions = gen_poisson_gp1_results.predict(X_test)
gen_poisson_gp1_predictions是一个pandas系列对象，包含X/test矩阵中每一行的预测自行车人数。记住，y/test包含实际观察到的人数。
让我们绘制预测数和实际数来直观地评估预测的质量：
predicted_counts=gen_poisson_gp1_prediction
sactual_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()
我们得到了以下预测与实际骑车人数的关系图：
来自GP-1模型的预测与实际骑行者数量的关系图
正如你所看到的，除了少数计数外，GP-1模型在预测骑自行车的人数方面做得相当好。
如果需要，我们可以比较GP-1和普通泊松的预测质量，方法是比较GP-1的预测的均方根误差（RMSE）和普通泊松模型对相同测试数据集的预测。我将把这作为一个练习。
最后，我们也来试试法莫耶的限制性广义泊松回归模型，即GP-2：
我们将使用同样的3步方法来建立和训练模型。 注意参数p=2：
#Build Famoye's Restricted Generalized Poison regression model, know as GP-2
gen_poisson_gp2 = sm.GeneralizedPoisson(y_train, X_train, p=2)
#Fit the model
gen_poisson_gp2_results = gen_poisson_gp2.fit()
#print the results
print(gen_poisson_gp2_results.summary())
我们看到以下训练输出：
GP-2模型的训练输出
分析和比较GP-2的拟合度和预测质量的方法与GP-1相同。然而请注意，在这种情况下，该模型的训练算法无法收敛。因此，我们将不再研究这个模型，而是选择GP-1而不是标准的泊松模型来为自行车人数数据建模。
摘要
-
标准泊松模型假定基于计数的数据的方差与平均值相同。这一假设经常被现实世界的数据集所违反，这些数据集要么过度分散，要么分散不足。 -
因此，我们需要对基于计数的数据使用其他模型，如负二项模型或广义泊松回归模型，这些模型并不假设数据是等距分散的。这类模型假定方差是平均值的某个函数。 -
Consul的广义泊松回归模型（称为GP-1）和Famoye的限制性广义泊松回归模型（GP-2）是两个这样的GP模型，可以用来为现实世界中基于计数的数据集建模。 -
Python库Statsmodels刚好对建立和训练GP-1和GP-2模型有很好的支持。