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

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

对照指数族分布的通式：

可得，

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

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

泊松回归是计数模型，因变量必须为自然数。一般针对具有稳定发生频数/率的事件进行建模。在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)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -4.6735  -0.9850  -0.1323   0.7891  21.2791
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.748568   0.001302    3648   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 9873.8  on 5113  degrees of freedom
## Residual deviance: 9873.8  on 5113  degrees of freedom
## AIC: 43525
##
## Number of Fisher Scoring iterations: 4

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

model.2 <- glm(death ~ temp, family = poisson(),data = chicagoNMMAPS)summary(model.2)
##
## Call:
## glm(formula = death ~ temp, family = poisson(), data = chicagoNMMAPS)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -4.3097  -0.8563  -0.0719   0.7530  22.5216
##
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)  4.7927697  0.0017346 2762.97   <2e-16 ***
## temp        -0.0044903  0.0001197  -37.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 9873.8  on 5113  degrees of freedom
## Residual deviance: 8471.8  on 5112  degrees of freedom
## AIC: 42125
##
## Number of Fisher Scoring iterations: 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)
##
## Deviance Residuals:
##     Min       1Q   Median       3Q      Max
## -3.8876  -0.8192  -0.0514   0.7029  22.6409
##
## Coefficients:
##                Estimate Std. Error  z value Pr(>|z|)
## (Intercept)   4.8028333  0.0017758 2704.577  < 2e-16 ***
## temp          0.0019664  0.0003720    5.286 1.25e-07 ***
## lag(temp, 1) -0.0011023  0.0005198   -2.121    0.034 *
## lag(temp, 2) -0.0021813  0.0005301   -4.115 3.87e-05 ***
## lag(temp, 3) -0.0005632  0.0005295   -1.064    0.287
## lag(temp, 4) -0.0008187  0.0005193   -1.577    0.115
## lag(temp, 5) -0.0028462  0.0003720   -7.651 1.99e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 9856.6  on 5108  degrees of freedom
## Residual deviance: 7768.2  on 5102  degrees of freedom
##   (5 observations deleted due to missingness)
## AIC: 41398
##
## Number of Fisher Scoring iterations: 4

偏移项

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

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

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

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

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

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

准泊松回归

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

## (Dispersion parameter for poisson family taken to be 1)

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

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))
##                 Estimate   Std. Error   z value      Pr(>|z|)
## (Intercept)  4.792769652 0.0017346450 2762.9687  0.000000e+00
## temp        -0.004490284 0.0001196991  -37.5131 5.633927e-308coef(summary(model.22))
##                 Estimate   Std. Error    t value      Pr(>|t|)
## (Intercept)  4.792769652 0.0023079821 2076.60605  0.000000e+00
## temp        -0.004490284 0.0001592622  -28.19428 1.102289e-162

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

相关阅读：

1. stats | 概率分布与随机数生成（一）——离散型分布

2. stats | 线性回归（一）——模型表达式和输出结果

3. stats | 广义线性模型（一）——广义线性模型的基本结构及与线性模型的比较

往期推荐阅读：

• 《数据处理通识》专辑-base | 使用apply族函数进行向量化运算

• 《制表与可视化》专辑-ggplot2 | ggplot2作图语法入门

• 《数学模型》专辑-car | 线性回归（三）——残差分析和异常点检验

• 《地理计算与分析》专辑-spdep | 如何在R语言中计算空间自相关指数