data18_1 <- haven::read_sav("datasets/例18-01.sav",encoding = "GBK")
data18_1 <- haven::as_factor(data18_1)
str(data18_1)
## tibble [8 × 7] (S3: tbl_df/tbl/data.frame)
## $ X1 : Factor w/ 2 levels "有暴露","无暴露": 1 2 1 2 1 2 1 2
## ..- attr(*, "label")= chr "砷暴露情况"
## $ YEARGRP: Factor w/ 4 levels "40-49岁组","50-59岁组",..: 1 1 2 2 3 3 4 4
## ..- attr(*, "label")= chr "年龄分组"
## $ X2 : num [1:8] 1 1 0 0 1 1 1 1
## ..- attr(*, "label")= chr "50-59岁组"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ X3 : num [1:8] 1 1 1 1 0 0 1 1
## ..- attr(*, "label")= chr "60-69岁组"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ X4 : num [1:8] 1 1 1 1 1 1 0 0
## ..- attr(*, "label")= chr ">=70岁组"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ Y : num [1:8] 7 14 42 38 59 58 17 41
## ..- attr(*, "label")= chr "死亡数"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ N : num [1:8] 11026 38337 10792 31019 6898 ...
## ..- attr(*, "label")= chr "观察单位数"
## ..- attr(*, "format.spss")= chr "F8.2"
data18_1
## # A tibble: 8 × 7
## X1 YEARGRP X2 X3 X4 Y N
## <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 有暴露 40-49岁组 1 1 1 7 11026.
## 2 无暴露 40-49岁组 1 1 1 14 38337.
## 3 有暴露 50-59岁组 0 1 1 42 10792.
## 4 无暴露 50-59岁组 0 1 1 38 31019.
## 5 有暴露 60-69岁组 1 0 1 59 6898.
## 6 无暴露 60-69岁组 1 0 1 58 17496.
## 7 有暴露 >=70岁组 1 1 0 17 2581.
## 8 无暴露 >=70岁组 1 1 0 41 6842.23 泊松回归和负二项回归
23.1 泊松回归简介
人类稀有疾病或一些卫生事件，如恶性肿瘤、非遗传性先天性疾病、癫痫患者在两周内癫痫病发作次数、某地在一个月内因交通事故死亡人数、某病患者在一年内住院次数、1ml水中大肠杆菌数、1L空气中粉尘粒子数和放射性物质在一定时间内放射性质点数等计数资料（count data），具有发病率低或者不像二项分布资料有分母能计算比例（proportion）等特点；因此，这些事件数的多少除了取决于事件的实际发生数外，还取决于计数时研究者所观察的范围，即观察多长时间、多大人群、多少体积或面积等。使用发病密度（incidence density，ID）等密度指标描述这些事件的群体特征比较合适。对于上述罕见事件的发生，如果事件之间彼此相互独立，观察样本含量较大时，则具有平均计数等于方差的特点；这类事件的发生次数往往服从Poisson分布。
Poisson回归（Poisson regression）主要用于单位时间、单位面积、单位空间内某事件发生数的影响因素分析。
23.2 泊松回归应用
孙振球《医学统计学》第4版例18-1。研究者为检查某冶炼厂的砷暴露与因呼吸道疾病死亡之间的关系，对该厂1978一2009年的职工进行了回顾性队列研究，其结果如下。请对该资料进行分析。
Y是因变量，X1和YEARGRP是自变量，X2/X3/X4是YEARGRP的哑变量编码形式，建模时用不到，因为软件会自动进行这一步，详情可见：分类变量进行回归分析时的编码方式。
X1的因子水平需要修改一下，让无暴露作为reference：
# 查看下因子水平
levels(data18_1$X1)
## [1] "有暴露" "无暴露"
levels(data18_1$YEARGRP)
## [1] "40-49岁组" "50-59岁组" "60-69岁组" ">=70岁组"
# 修改下X1的因子水平
data18_1$X1 <- factor(data18_1$X1,levels = c("无暴露","有暴露"))下面就是建立泊松回归模型即可。本例资料的观察单位为人年数，事件数（因变量）为因呼吸道疾病死亡人数；影响因素有两个，砷暴露情况和年龄。由于观察单位不同，所以在建模时需要设置偏移量，也就是offset：
f <- glm(Y ~ X1 + YEARGRP,data = data18_1,family = poisson(),offset = log(N))
# 或者也可以写成
# glm(Y ~ X1 + YEARGRP+offset(log(N)), data = data18_1, family = poisson())
#f
summary(f)
##
## Call:
## glm(formula = Y ~ X1 + YEARGRP, family = poisson(), data = data18_1,
## offset = log(N))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.0086 0.2233 -35.859 < 2e-16 ***
## X1有暴露 0.8109 0.1210 6.699 2.09e-11 ***
## YEARGRP50-59岁组 1.4702 0.2453 5.994 2.04e-09 ***
## YEARGRP60-69岁组 2.3661 0.2372 9.976 < 2e-16 ***
## YEARGRP>=70岁组 2.6238 0.2548 10.297 < 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: 260.9304 on 7 degrees of freedom
## Residual deviance: 9.9303 on 3 degrees of freedom
## AIC: 61.342
##
## Number of Fisher Scoring iterations: 4常数项和各变量的回归系数和书中是一样的，但是标准误有所不同(导致可信区间也不同)，推测是由于各种计算过程中的细节差异导致，并不影响结果。z值是系数除以标准误得到的。
# 单独查看结果以及其他结果
coef(f) # 查看回归系数
## (Intercept) X1有暴露 YEARGRP50-59岁组 YEARGRP60-69岁组
## -8.0086495 0.8108698 1.4701505 2.3661111
## YEARGRP>=70岁组
## 2.6237532
confint(f) # 回归系数的标准误
## 2.5 % 97.5 %
## (Intercept) -8.4780952 -7.598323
## X1有暴露 0.5724818 1.047494
## YEARGRP50-59岁组 1.0094838 1.976026
## YEARGRP60-69岁组 1.9239329 2.858522
## YEARGRP>=70岁组 2.1411250 3.145495
exp(coef(f)) # 相对危险度RR值（或者叫发病率比值，IRR）
## (Intercept) X1有暴露 YEARGRP50-59岁组 YEARGRP60-69岁组
## 3.325736e-04 2.249864e+00 4.349890e+00 1.065587e+01
## YEARGRP>=70岁组
## 1.378737e+01
exp(confint(f)) # RR值的可信区间
## 2.5 % 97.5 %
## (Intercept) 0.0002079745 5.012913e-04
## X1有暴露 1.7726609281 2.850500e+00
## YEARGRP50-59岁组 2.7441841472 7.214017e+00
## YEARGRP60-69岁组 6.8478373075 1.743573e+01
## YEARGRP>=70岁组 8.5090045641 2.323118e+01结果解释：在控制年龄因素后，砷暴露组因呼吸道疾病死亡的风险是非暴露组的exp(0.8109)=2.25倍（95%CI：1.77～2.85）。在控制砷暴露因素后，因呼吸道疾病死亡风险随着年龄增加越来越大，其中“50～59岁”年龄组因呼吸道疾病死亡风险是“40～49岁”年龄组的exp(1.4702)=4.35倍（95%CI：2.747.21)，“60～69岁”年龄组因呼吸道疾病死亡风险是“40～49岁”年龄组的exp(2.366)=10.66倍（95%CI：6.85～17.44），70岁以上年龄组因呼吸道疾病死亡风险是“40～49”岁年龄组的exp(2.6238)=13.79倍（95%CI:8.5123.23)。
基线发病密度（即所有自变量均为0时的发病密度）或40～49岁无砷暴露史人群因呼吸道疾病的死亡密度为exp(-8.0086)=33.260/10万人年。40～49岁有砷暴露史的人群因呼吸道疾病的死亡密度为exp(-8.0086+0.8109×1+1.4702×0+2.3661x0+2.6238x0)=74.831/10万人年；其他各组的死亡密度计算与此类似。
下面对模型进行拟合优度检验，采用pearson卡方统计量和残差偏差统计量进行评价，参考资料见（https://rstudio-pubs-static.s3.amazonaws.com/1047952_9306ae04c1de4543812af559d777dd72.html）。
首先是残差偏差法：
# 计算P值
1 - pchisq(deviance(f), df = f$df.residual)
## [1] 0.01916785结果小于0.05，说明拟合不好。
再看pearson卡方的P值：
Pearson <- sum((data18_1$Y - f$fitted.values)^2 / f$fitted.values)
1 - pchisq(Pearson, df = f$df.residual)
## [1] 0.02137005结果也是小于0.05，说明拟合不好。书中的方法与此不同，其结果都是P值大于0.25，但其实个人感觉这个模型确实不太好。
下面我们对这个模型进行过度离散（overdispersion）检验，如果存在过度离散，说明使用该方法是不太合适的。
library(performance)
check_overdispersion(f)
## # Overdispersion test
##
## dispersion ratio = 3.231
## Pearson's Chi-Squared = 9.692
## p-value = 0.021结果表明该数据确实存在过度离散，因此使用poisson()是不太合适的，我们可以使用quasipoisson()试试（或者用下一节介绍的负二项回归）。
泊松回归假设方差等于均值，也就是所谓的等离散（equidispersion）。而实际数据中经常出现过度离散（overdispersion），也就是方差大于均值的情况。这时候，使用poisson族可能不合适，标准误会被低估，导致p值过小，可能得出错误的显著性结论。quasipoisson就是用来处理这种情况的，它引入了一个离散参数来调整方差和标准误，但不改变系数估计值。
f1 <- glm(Y ~ X1 + YEARGRP,data = data18_1,family = quasipoisson(),offset = log(N))
summary(f1)
##
## Call:
## glm(formula = Y ~ X1 + YEARGRP, family = quasipoisson(), data = data18_1,
## offset = log(N))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -8.0086 0.4014 -19.950 0.000275 ***
## X1有暴露 0.8109 0.2176 3.727 0.033640 *
## YEARGRP50-59岁组 1.4702 0.4408 3.335 0.044555 *
## YEARGRP60-69岁组 2.3661 0.4263 5.550 0.011534 *
## YEARGRP>=70岁组 2.6238 0.4580 5.729 0.010558 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasipoisson family taken to be 3.23081)
##
## Null deviance: 260.9304 on 7 degrees of freedom
## Residual deviance: 9.9303 on 3 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 4系数没变，P值和标准误变大了，结果解读同上，不再重复。
23.3 负二项回归简介
Poisson回归的应用条件之一是计数资料服从Poisson分布，即满足计数值的平均值等于方差。但是，有许多事件的发生是非独立的（如传染性疾病、遗传性疾病、地方病、致病生物的分布和一些原因不明疾病的空间聚集现象等），它们的计数资料会发生方差远远大于平均值，即存在过度离散现象（over-dispersion）；若用Poisson回归来分析这些事件的影响因素，会导致模型参数估计值的标准误偏小，参数检验的假阳性率增加。此时，宜选用负二项回归模型（negativebinomialregression）来分析这些资料。
负二项回归模型是基于计数资料服从负二项分布的。负二项分布实际上是当Poisson分布中强度参数λ服从γ分布（Gamma distribution）时所得到的复合分布。在Poisson分布中，λ是一个常数：在负二项分布中，λ是一个随机变量，并服从γ分布；因此，负二项分布又称γ-Poisson分布（Gamma-Poisson distribution）。
23.4 负二项回归应用
孙振球《医学统计学》第4版例18-2。某学者为了检查居住地类型与蚊虫幼虫滋生的关系，对299个不同居住地的家庭进行调查，结果如下。请对该资料选择合适的统计方法进行分析。
data18_2 <- haven::read_sav("datasets/例18-02.sav",encoding = "GBK")
data18_2 <- haven::as_factor(data18_2)
str(data18_2)
## tibble [24 × 3] (S3: tbl_df/tbl/data.frame)
## $ x : Factor w/ 3 levels "农村","城市贫民区",..: 1 1 1 1 1 1 1 1 2 2 ...
## ..- attr(*, "label")= chr "居住地情况"
## $ y : num [1:24] 0 1 2 3 4 5 6 11 0 1 ...
## ..- attr(*, "label")= chr "受蚊子幼虫滋生的容器数"
## ..- attr(*, "format.spss")= chr "F8.0"
## $ count: num [1:24] 136 23 10 5 2 1 1 1 38 8 ...
## ..- attr(*, "label")= chr "频数"
## ..- attr(*, "format.spss")= chr "F8.0"
head(data18_2)
## # A tibble: 6 × 3
## x y count
## <fct> <dbl> <dbl>
## 1 农村 0 136
## 2 农村 1 23
## 3 农村 2 10
## 4 农村 3 5
## 5 农村 4 2
## 6 农村 5 1这是频数资料的格式，不是原始数据格式，我们先把这个数据变成原始数据的格式：
x1 <- rep(data18_2$x,data18_2$count)
y1 <- rep(data18_2$y,data18_2$count)
df <- data.frame(x1,y1)
str(df) # 原始数据的结构
## 'data.frame': 299 obs. of 2 variables:
## $ x1: Factor w/ 3 levels "农村","城市贫民区",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ y1: num 0 0 0 0 0 0 0 0 0 0 ...
xtabs(~x1+y1,data = df) # 检查下和原来的数据是否一样
## y1
## x1 0 1 2 3 4 5 6 11
## 农村 136 23 10 5 2 1 1 1
## 城市贫民区 38 8 2 0 0 0 0 0
## 城市 67 5 0 0 0 0 0 0使用glm.nb()函数（来自MASS包）来拟合负二项回归模型：
library(MASS)
f <- glm.nb(y1~x1, data = df)
summary(f)
##
## Call:
## glm.nb(formula = y1 ~ x1, data = df, init.theta = 0.3002652205,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.7100 0.1731 -4.102 4.1e-05 ***
## x1城市贫民区 -0.6762 0.4274 -1.582 0.113612
## x1城市 -1.9572 0.5256 -3.724 0.000196 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(0.3003) family taken to be 1)
##
## Null deviance: 174.95 on 298 degrees of freedom
## Residual deviance: 156.37 on 296 degrees of freedom
## AIC: 426.23
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 0.3003
## Std. Err.: 0.0764
##
## 2 x log-likelihood: -418.2280结果与书中一致。输出结果和泊松回归的结果基本一样，其中有一个Theta：0.3003是负二项分布的过度离散参数，表示数据的过度离散程度。如果theta接近1，表明数据接近泊松分布。
结果显示城市家庭和贫民区家庭滋生蚊虫幼虫机会都低于农村家庭，但是只有城市家庭才有统计学意义。或者，农村家庭滋生蚊虫幼虫机会是城市家庭的exp（1.9572）=7.08倍（95%CI：2.527～19.834），而与贫民区家庭之间没有差别。
# 演示下这个95%CI的计算方法。
# 先计算系数的可信区间：均值±1.96*标准误
cbind(-1.9572+1.96*0.5256,-1.9572-1.96*0.5256) # 系数的可信区间
## [,1] [,2]
## [1,] -0.927024 -2.987376
cbind(exp(0.9270),exp(2.9874)) # RR值的可信区间
## [,1] [,2]
## [1,] 2.526917 19.8340523.5 参考资料
- When to use an offset in a Poisson regression?
- https://rstudio-pubs-static.s3.amazonaws.com/1047952_9306ae04c1de4543812af559d777dd72.html