R语言统计分析——泊松回归
R语言统计分析——泊松回归
泊松回归是一种用于预测计数型结果变量的统计分析方法，特别适用于通过一系列连续型和/或类别型预测变量来预测计数型结果变量的场景。本文将通过一个具体的案例（Breslow癫痫数据）来阐述泊松回归模型的拟合过程，并探讨一些可能出现的问题。
数据准备
我们使用robust包中的Breslow癫痫数据集，该数据集记录了遭受轻微或严重间歇性癫痫的病人的年龄和癫痫发病数。数据包含病人被随机分配到药物组或者安慰剂组前八周和随机分配后八周两种情况。响应变量为sumY（随机化后八周内癫痫发病数），预测变量为治疗条件（Trt）、年龄（Age）和前八周内的基础癫痫发病数。
首先，加载数据集并查看其统计汇总信息：
# 加载数据集
data(breslow.dat, package = "robust")
# 查看数据列
names(breslow.dat)
# 查看统计汇总数据
summary(breslow.dat[c(6, 7, 8, 10)])
数据可视化
接下来，对数据进行可视化，以直观了解数据的分布特征：
# 对数据进行可视化
opar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
attach(breslow.dat)
hist(sumY, breaks = 20, xlab = "Seizure Count", main = "Distribution of Seizure")
boxplot(sumY ~ Trt, xlab = "Treatment", main = "Group comparisons")
par(opar)
detach(breslow.dat)
从直方图可以看到因变量的偏倚特性以及可能的离群点。从箱线图可以看到，药物治疗下癫痫发病数似乎变小了，且方差也变小了（泊松分布中，较小的方差伴随着较小的均值）。与标准最小二乘回归（OLS）不同，泊松回归并不关注方差异质性。
泊松回归拟合
使用glm函数拟合泊松回归模型：
# 拟合泊松回归
fit <- glm(sumY ~ Base + Age + Trt, family = poisson(link = "log"), data = breslow.dat)
# 查看拟合结果
summary(fit)
输出结果列出了偏差、回归参数、标准误和参数为0的检验。注意，此处预测变量在p<0.05的水平下都非常显著。
解释模型参数
使用coef()函数可以获取模型系数，或者调用summary()函数的输出结果中的Coefficients表格。
在泊松回归中，因变量以条件均值的对数形式来建模。年龄的回归参数为0.02274，表明保持其他预测变量不变，年龄增加一岁，癫痫病发病的对数均值响应的0.03。截距项即当预测变量都为0时，癫痫发病的对数均值。由于不可能为0岁，且调查对象的基础癫痫发病数均不为0，因此本例中的截距项没有实际意义。
通常在因变量的初始尺度（癫痫病发病数，而非发病数的对数）上解释回归系数比较容易，因此，可以进行指数化系数：
现在可以看到，保持其他变量不变，年龄增加一岁，期望的癫痫发病数将乘以1.023。这意味着年龄的增加与较高的癫痫发病数相关联。更为重要的是，一单位Trt的变化（即从安慰剂到治疗组），期望的发病数将乘以0.86，也就是说，保持基础癫痫发病数和年龄不变，服药组相对于安慰剂组癫痫病发病数降低了20%。
过度离势
泊松分布的方差和均值相等。当响应变量观测的方差比依据泊松分布预测的方差大时，泊松回归可能发生过度离势。处理计数型数据时经常发生过度离势，且过度历史会对结果的可解释性造成负面影响。
发生过度离势可能以下几个原因：
- 遗漏了某个重要的预测变量。
- 可能因为事件相关。在泊松分布的观测中，计数中每次事件都被认为是独立发生的。以癫痫数据为例，这意味着对于任何病人，每次癫痫发病的概率与其他癫痫发病的概率相互独立。但是这个假设通常都无法满足。对于某个病人，在已知他已经发生了39次癫痫时，第一次发生癫痫的概率不可能与第40次发生癫痫的概率相同。
- 在纵向数据分析中，重复测量的数据由于内在群聚特性可导致过度离势。此处并不讨论纵向泊松模型。
如果存在过度离势，在模型中我们无法解释，那么可能会得到很小的标准误和置信区间，并且显著性检验也过于宽松。
与Logistic回归类似，泊松回归如果残差偏差与残差自由度的比例远大于1，那么表明存在过度离势。对于癫痫数据，它的比例是：
显然比例明显大于1。qcc包提供了一个对泊松模型过度离势的检验方法，代码如下：
# 安装qcc包
install.packages("qcc")
# 加载qcc包
library(qcc)
# 泊松回归过度离势检验
qcc.overdispersion.test(breslow.dat$sumY, type = "poisson")
检验结果显示，p小于0.05，进一步表明确实存在过度离势。
通过family = quasipoisson替换family = poisson，我们仍然可以使用glm()函数对该数据进行拟合：
fit.od <- glm(sumY ~ Base + Age + Trt, family = quasipoisson(link = "log"), data = breslow.dat)
# 查看结果
summary(fit.od)
注意，使用类泊松（quasi-Poisson）方法所得的参数估计与泊松方法相同，单标准误变大了许多。此处，标准误越大将会导致Trt和Age的p值变大。当考虑过度离势，并控制技术癫痫数和年龄时，并没有充足的证据表明药物质量相对于使用安慰剂能显著降低癫痫发病次数。