stats | 广义线性模型（二）—

stats | 广义线性模型（二）——泊松回归

本篇介绍一种常见的广义线性模型：泊松回归。泊松分布是离散型分布，它的概率分布函数如下：

写成指数族分布的形式如下：

对照指数族分布的通式：

可得，

广义线性模型假设与解释变量存在线性关系，即

又因为泊松分布的数学期望就是，因此泊松回归可写成如下形式：

泊松回归是计数模型，因变量必须为自然数。一般针对具有稳定发生频数/率的事件进行建模。在glm函数中，泊松回归的family参数设置为poisson(link = "log")或poisson()。

示例

示例数据是来自dlnm工具包的chicagoNMMAPS数据集，它记录的是芝加哥逐日的死亡人口数和温度、环境污染物浓度等信息。

library(dlnm)names(chicagoNMMAPS)##[1]"date""time""year""month""doy""dow""death""cvd""resp"##[10]"temp""dptp""rhum""pm10""o3"

在一段时间内，自然状态下地区的死亡人口可以看作是等频数发生的，因此可以使用泊松回归进行建模。

model<-glm(death~1,family=poisson(),data=chicagoNMMAPS)summary(model)####Call:##glm(formula=death~1,family=poisson(),data=chicagoNMMAPS)####DevianceResiduals:##Min1QMedian3QMax##-4.6735-0.9850-0.13230.789121.2791####Coefficients:##EstimateStd.ErrorzvaluePr(>|z|)##(Intercept)4.7485680.0013023648<2e-16***##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1####(Dispersionparameterforpoissonfamilytakentobe1)####Nulldeviance:9873.8on5113degreesoffreedom##Residualdeviance:9873.8on5113degreesoffreedom##AIC:43525####NumberofFisherScoringiterations:4

在模型中添加温度变量temp作为解释变量：

model.2<-glm(death~temp,family=poisson(),data=chicagoNMMAPS)summary(model.2)####Call:##glm(formula=death~temp,family=poisson(),data=chicagoNMMAPS)####DevianceResiduals:##Min1QMedian3QMax##-4.3097-0.8563-0.07190.753022.5216####Coefficients:##EstimateStd.ErrorzvaluePr(>|z|)##(Intercept)4.79276970.00173462762.97<2e-16***##temp-0.00449030.0001197-37.51<2e-16***##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1####(Dispersionparameterforpoissonfamilytakentobe1)####Nulldeviance:9873.8on5113degreesoffreedom##Residualdeviance:8471.8on5112degreesoffreedom##AIC:42125####NumberofFisherScoringiterations:4

在模型中添加温度的时间滞后项作为解释变量：

library(dplyr)model.3<-glm(death~temp+lag(temp,1)+lag(temp,2)+lag(temp,3)+lag(temp,4)+lag(temp,5),family=poisson(),data=chicagoNMMAPS)summary(model.3)####Call:##glm(formula=death~temp+lag(temp,1)+lag(temp,2)+lag(temp,##3)+lag(temp,4)+lag(temp,5),family=poisson(),data=chicagoNMMAPS)####DevianceResiduals:##Min1QMedian3QMax##-3.8876-0.8192-0.05140.702922.6409####Coefficients:##EstimateStd.ErrorzvaluePr(>|z|)##(Intercept)4.80283330.00177582704.577<2e-16***##temp0.00196640.00037205.2861.25e-07***##lag(temp,1)-0.00110230.0005198-2.1210.034*##lag(temp,2)-0.00218130.0005301-4.1153.87e-05***##lag(temp,3)-0.00056320.0005295-1.0640.287##lag(temp,4)-0.00081870.0005193-1.5770.115##lag(temp,5)-0.00284620.0003720-7.6511.99e-14***##---##Signif.codes:0'***'0.001'**'0.01'*'0.05'.'0.1''1####(Dispersionparameterforpoissonfamilytakentobe1)####Nulldeviance:9856.6on5108degreesoffreedom##Residualdeviance:7768.2on5102degreesoffreedom##(5observationsdeletedduetomissingness)##AIC:41398####NumberofFisherScoringiterations:4

偏移项

在很多情况下，保持恒定的并非计数变量（如示例中的死亡人口），而是比率变量。在上述示例中，死亡人口受到总人口或老年人口（如60岁以上人口）数量的影响：

泊松回归要求因变量必须为非负整数，比率变量显然不能满足这一条件。

利用对数的运算法则，有如下变形：

在上式中，log(pop)作为系数恒为1的解释变量，即偏移项。因此可以设置如下形式的模型：

model.4<-glm(death~temp,offset=log(pop),family=poisson(),data=chicagoNMMAPS)

chicagoNMMAPS数据集中无pop变量，上式只作为演示，不能运行。

准泊松回归

泊松分布的均值和方差相等，均为：

##(Dispersionparameterforpoissonfamilytakentobe1)

但很多情况下该条件并满足：

deviance(model.2)/df.residual(model.2)##[1]1.657238

这种情况可以使用准泊松分布族quasipoisson()：

model.22<-glm(death~temp,family=quasipoisson(),data=chicagoNMMAPS)summary(model.22)

结果对比：

coef(summary(model.2))##EstimateStd.ErrorzvaluePr(>|z|)##(Intercept)4.7927696520.00173464502762.96870.000000e+00##temp-0.0044902840.0001196991-37.51315.633927e-308coef(summary(model.22))##EstimateStd.ErrortvaluePr(>|t|)##(Intercept)4.7927696520.00230798212076.606050.000000e+00##temp-0.0044902840.0001592622-28.194281.102289e-162

泊松回归和准泊松回归的系数估计值是一样的，但后者的标准误（Std. Error）较大。