第 37 章 计数型因变量

计数型变量在临床医学/流行病学研究中也十分常见，下面是一些例子：

  1. 某个呼吸科诊所的患者中，每个人在过去一个月中哮喘发作的次数；
  2. 癫痫患者在过去一年中癫痫发作次数；
  3. 接受脑部 CT 扫描的患者中，每个人被诊断出颅内肿瘤个数。

最早的泊松模型可以追溯到普鲁士骑兵连中被马蹄踢死士兵的人数模型。

37.1 泊松 GLM

一个计数型的随机变量，只能取大于等于零的正整数，\(0,1,\cdots\)。泊松随机变量可以理解为产生于发生在一段时间内的事件次数。泊松模型可以用于计数型数据的回归模型的构建：

\[ \begin{aligned} Y &\sim \text{Po}(\mu) \\ \text{P} (Y = y) & = \frac{\mu^y e^{-\mu}}{y!} \end{aligned} \]

所以，一个泊松回归，默认的前提是因变量 \(Y\) 服从一个以预测变量 \(x_1, \cdots, x_p\) 为条件的泊松分布。其标准链接方程是 \(\theta=\text{log}(\mu)\)。

\[ \begin{aligned} Y_i & \sim \text{Po}(\mu_i) \\ \text{log}(\mu_i) & = \alpha + \beta_1 x_{i1} + \cdots + \beta_p x_{ip} \end{aligned} \]

观测对象1，用模型中全部的预测变量\(\mathbf{x_1}=(x_{11},\cdots,x_{1p})\) 计算获得的拟合值，和另一个观测对象0 的拟合值之比为：

\[ \begin{aligned} & \frac{\text{exp}(\alpha + \beta_1 x_{11} + \cdots + \beta_p x_{1p})}{\text{exp}(\alpha + \beta_1 x_{01} + \cdots + \beta_p x_{0p})} \\ = & \exp(\beta_1(x_{11}-x_{01}) + \cdots + \beta_p(x_{1p} - x_{0p})) \end{aligned} \]

其中，

  • 线性预测方程linear predictor 中的截距\(\alpha\) 的含义是，当所有的预测变量均等于零\(\mathbf{x_1} = 0\) 时，因变量\(Y\) 的均值之对数。
  • \(\beta_1\) 的含义是，其余预测变量保持不变时，预测变量 \(x_1\) 每增加一个单位时，因变量变化量的对数。
  • 回归系数的指数 (自然底数) 大小，可以被理解为是率比 (rate ratio) (详见下一章率的 GLM)。

37.2 泊松回归实例

下列数据来自 UCLA 的统计学网站。数据内容是某高中全部学生，获奖的次数。预测变量包括，1) 获奖种类 “一般 General”，“学术类 Academic”，“技能类 Vocational”；和所有学生期末数学考试分数。

p <- read_csv("../Datas/poisson_sim.csv")
p <- within(p, {
  prog <- factor(prog, levels=1:3, labels=c("General", "Academic",
                                                     "Vocational"))
  id <- factor(id)
})
summary(p)
##        id        num_awards           prog          math      
##  1      :  1   Min.   :0.00   General   : 45   Min.   :33.00  
##  2      :  1   1st Qu.:0.00   Academic  :105   1st Qu.:45.00  
##  3      :  1   Median :0.00   Vocational: 50   Median :52.00  
##  4      :  1   Mean   :0.63                    Mean   :52.65  
##  5      :  1   3rd Qu.:1.00                    3rd Qu.:59.00  
##  6      :  1   Max.   :6.00                    Max.   :75.00  
##  (Other):194

下面的代码拟合因变量为获奖次数，预测变量为获奖种类(分类) 和数学成绩(连续) 的泊松分布，泊松分布默认的链接方程就是\(\text{log}\)，所以你可以像第一行那样把链接方程部分省略。结果也是一样的。

m1 <- glm(num_awards ~ prog, family="poisson", data=p)
m2 <- glm(num_awards ~ prog, family=poisson(link = log), data=p)
summary(m1); summary(m2)
## 
## Call:
## glm(formula = num_awards ~ prog, family = "poisson", data = p)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.6094     0.3333  -4.828 1.38e-06 ***
## progAcademic     1.6094     0.3473   4.634 3.59e-06 ***
## progVocational   0.1823     0.4410   0.413    0.679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 234.46  on 197  degrees of freedom
## AIC: 416.51
## 
## Number of Fisher Scoring iterations: 6
## 
## Call:
## glm(formula = num_awards ~ prog, family = poisson(link = log), 
##     data = p)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.6094     0.3333  -4.828 1.38e-06 ***
## progAcademic     1.6094     0.3473   4.634 3.59e-06 ***
## progVocational   0.1823     0.4410   0.413    0.679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 234.46  on 197  degrees of freedom
## AIC: 416.51
## 
## Number of Fisher Scoring iterations: 6

输出结果的回归系数部分，

  • 该学校学生获得学术类奖项的平均次数和获得一般奖项的平均次数的比值是\(\text{exp}(1.6094) = 4.999\)，所以获得的学术类奖平均次数要高于一般奖次数\(390\ %\)；
  • 获得技能类奖的平均次数和一般奖平均次数的比值是 \(\text{exp}(0.1823) = 1.199\)，也就是高出了 \(19.9\%\)；
  • 该校学生获得一般类奖的次数平均每人是 \(\text{exp}(-1.6094) = 0.20\) 次；
  • 该校学生获得学术奖的次数平均每人是 \(\text{exp}(-1.6094 + 1.6094) = 1\) 次；(一人一次够流弊)
  • 该校学生获得技能类奖的次数平均每人是 \(\text{exp}(-1.6094 + 0.182) = 0.24\) 次。

看来该校师生很重视学术。

当然也可以用下面定义的函数来帮助我们计算上面这些数值，及其置信区间。

glm.RR <- function(GLM.RESULT, digits = 2) {

    if (GLM.RESULT$family$family == "binomial") {
        LABEL <- "OR"
    } else if (GLM.RESULT$family$family == "poisson") {
        LABEL <- "RR"
    } else {
        stop("Not logistic or Poisson model")
    }

    COEF      <- stats::coef(GLM.RESULT)
    CONFINT   <- stats::confint(GLM.RESULT)
    TABLE     <- cbind(coef=COEF, CONFINT)
    TABLE.EXP <- round(exp(TABLE), digits)

    colnames(TABLE.EXP)[1] <- LABEL

    TABLE.EXP
}
glm.RR(m1)
##                 RR 2.5 % 97.5 %
## (Intercept)    0.2  0.10   0.36
## progAcademic   5.0  2.68  10.63
## progVocational 1.2  0.51   2.94

37.3 过度离散 overdispersion

泊松分布的前提条件之一是，方差和均值相等。这是一个非常强的假设，很多计数型数据其实是无法满足这个条件的。许多时候 (包括上面的例子也是) 方差要大于或者小于均值：

epiDisplay::summ(p$num_awards[p$prog == "Academic"], graph = FALSE)
##  obs. mean   median  s.d.   min.   max.  
##  105  1      1       1.279  0      6
epiDisplay::summ(p$num_awards[p$prog == "General"], graph = FALSE)
##  obs. mean   median  s.d.   min.   max.  
##  45   0.2    0       0.405  0      1
epiDisplay::summ(p$num_awards[p$prog == "Vocational"], graph = FALSE)
##  obs. mean   median  s.d.   min.   max.  
##  50   0.24   0       0.517  0      2

试想一下，实际的数据中其实是经常出现这样的违反泊松分布前提的计数型数据的。例如某两个观测对象，如果他们二者的线性预测方程给出相等的结果(他们各自的预测变量可以完全不同)，会被认为服从相同均值，相同方差的泊松分布，这显然是不合理的。例如本章用到的学校学生获奖的例子，有的学生成绩好，那么获得学术类奖的平均次数(及其方差) 自然和成绩排在后面的学生不同，强制这样的两个学生服从相同均值，相同方差的泊松分布显然是不合情理的。手工好的学生，可能更倾向于获得更多得技能类奖。实际情况下，还有许许多多其他的未知因素会影响学生获奖的次数，例如家庭教育背景的不同，有些学生钢琴获奖多，因为他每天都去练习弹钢琴等等，这些都是没有被收集到的数据。

真实情况应该是这样的，当有其他的我们不知道的因素存在时，这些因素会导致某些人的均值高于其他人。如果对象\(i\) 的因变量\(Y_i\) 服从均值为\(\mu_i\) 的泊松分布，那么对于所有的\(\mu_i\)，其均值(overall mean) 是\(\mu\)，方差(overall variance) 是\(\sigma^2\)。这是一个典型的随机效应模型random effect model，我们会在后面的hierarchical data analysis 再深入讨论，但是这里的重点是，每个观测对象自己的均值\(\mu_i\)，是我们在普通泊松回归中忽略掉的随机共变量(the effects of omitted covariates)。

所以样本数据来自的人群如果共同均值 (或者叫边际效应均值，marginal mean) 为 \(\mu\)：

\[ E(Y_i) = E(E(Y_i | \mu_i)) = E(\mu_i) = \mu \]

和共同方差 (边际效应方差) ，需要用到 总体方差法则 (Law of total variance) 概念：

\[ \begin{aligned} \text{Var}(Y_i) & = E(\text{Var}(Y_i | \mu_i)) + \text{Var}(E(Y_i | \mu_i)) \\ & = E(\mu_i) + \text{Var}(\mu_i) \\ & = \mu + \sigma^2 \end{aligned} \]

37.3.1 过度离散怎么查？

如果，我们的泊松回归模型中的共变量全部都是分类型变量，我们可以把观测值 \(Y\) 对每一个分类变量分别作简单的数据总结，观察其均值和方差是否可以认为大致相同。但是许多时候模型中不会只有分类型变量。

R 输出的结果中的 模型偏差 deviance，可以用来初步判断整体模型的拟合优度。如果模型偏差除以残差获得的残差偏差(residual deviance) 足够小，说明拟合的模型跟数据本身比较接近，也就是模型和数据拟合程度较好，反之则提示模型本身具有较高的过度离散overdispersion。另外，模型偏差由于在个人数据 (individual data) 情况下不适用 (因为模型偏差值就不再服从卡方分布了)，下面的检验结果仅仅只能作为极为微弱的参考证据。此时应该推荐使用 Pearson 的模型拟合检验。如果 Pearson 统计量，除以残差的自由度获得的值远大于 1，就提示存在过度离散。

with(m1, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance  df          p
## [1,]       234.46 197 0.03496117

Goodness of fit 检验结果 提示本模型可能存在过度离散，数据拟合度不理想。值得注意的是如果样本很大时，模型偏差的检验统计量将不再服从卡方分布，应用的时候一定要慎重。

37.3.2 负二项式分布模型 negative binomial model

如果普通泊松回归模型拟合数据时，发现数据本身有过度离散的嫌疑，那么建议使用负二项式分布模型来重新拟合数据。负二项式分布模型其实是泊松分布的扩展版本，即考虑了个体的方差和均值的随机效应 subject-specific random effect。如果设每个观测对象的随机效应部分为\(a_i\)，预测变量为向量\(\mathbf{x_i} = (x_{i1}, \cdots, x_{ip})\)，那么因变量\(Y_i\) 服从均值为\(\text{exp}(\beta^T\mathbf{x_i}+a_i)\) 泊松分布。在负二项式分布中，个体的随机效应部分的自然底数的指数\(e^{a_i}\) 其实是服从均值为1， 方差为\(\alpha\) 的伽马分布(gamma distribution)。 \(\alpha\) 越大，说明过度离散越明显。

接下来用相同的数据，使用负二项式分布模型在 R 里作模型的拟合，你就会看到差别：

R 里拟合负二项式分布模型的函数 glm.nb 在基本包 MASS 里。

m1 <- MASS::glm.nb(num_awards ~ prog, data = p)
m2 <- glm(num_awards ~ prog, family=poisson(link = log), data=p)
summary(m1)
## 
## Call:
## MASS::glm.nb(formula = num_awards ~ prog, data = p, init.theta = 1.72267107, 
##     link = log)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.6094     0.3522  -4.570 4.87e-06 ***
## progAcademic     1.6094     0.3729   4.316 1.59e-05 ***
## progVocational   0.1823     0.4679   0.390    0.697    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(1.7227) family taken to be 1)
## 
##     Null deviance: 211.26  on 199  degrees of freedom
## Residual deviance: 171.07  on 197  degrees of freedom
## AIC: 406.53
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  1.723 
##           Std. Err.:  0.717 
## 
##  2 x log-likelihood:  -398.532
summary(m2)
## 
## Call:
## glm(formula = num_awards ~ prog, family = poisson(link = log), 
##     data = p)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -1.6094     0.3333  -4.828 1.38e-06 ***
## progAcademic     1.6094     0.3473   4.634 3.59e-06 ***
## progVocational   0.1823     0.4410   0.413    0.679    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 287.67  on 199  degrees of freedom
## Residual deviance: 234.46  on 197  degrees of freedom
## AIC: 416.51
## 
## Number of Fisher Scoring iterations: 6

仔细比较普通泊松分布回归和负二项式分布回归的输出结果，你会发现

  1. 回归系数的计算是完全相同的 (由于我们只放了一个简单的分类型变量作为预测变量，一般来说泊松回归和负二项式分布回归计算的回归系数会有些许不同)；
  2. 另外一个变化是标准误的估计量在负二项式分布模型中明显变大了，这就是我们放宽了前提条件，允许模型考虑个体的随机效应的体现。如果泊松模型被数据本身的过度离散影响显著，那么泊松回归计算获得的参数标准无是偏低的；
  3. 负二项式分布回归的结果最底下出现的 Theta: 1.723 部分，它的倒数是前面提到的个体的随机效应部分 \(a_i\) 服从的伽马分布的方差 \(\alpha\)。它是关键的离散程度参数 (dispersion parameter)。 在STATA 里，如果用nbreg 拟合负二项式分布回归的模型，输出的结果最底下会有\(\alpha\) 值的报告，注意它和R 输出的Theta 结果互为倒数。另外，STATA 的输出结果还会对 \(\alpha = 0\) 直接进行检验。在 R 里面则需要给两个模型分别进行拟合优度检验，多数情况下你会发现负二项式分布回归的模型更加拟合数据：
with(m1, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance  df         p
## [1,]     171.0661 197 0.9090187
with(m2, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance  df          p
## [1,]       234.46 197 0.03496117

另一种获取没有被低估的回归系数的标准误的方法来自稳健统计学手段。在 R 里，拟合完普通泊松回归以后，用 sandwich 包里的 vcovHC() 命令进行稳健的参数误差估计 (具体说是夹心方差矩阵估计 sandwich estimator of variance)：

m2 <- glm(num_awards ~ prog, family=poisson(link = log), data=p)
cov.m2 <- sandwich::vcovHC(m2, type = "HC0")
std.err <- sqrt(diag(cov.m2))
robust.est <- cbind(Estimate= coef(m2), "Robust SE" = std.err,
"Pr(>|z|)" = 2 * pnorm(abs(coef(m2)/std.err), lower.tail=FALSE),
LL = coef(m1) - 1.96 * std.err,
UL = coef(m1) + 1.96 * std.err)
robust.est
##                  Estimate Robust SE     Pr(>|z|)         LL        UL
## (Intercept)    -1.6094379 0.2981424 6.730573e-08 -2.1937970 -1.025079
## progAcademic    1.6094379 0.3229681 6.251792e-07  0.9764204  2.242455
## progVocational  0.1823216 0.4242641 6.673877e-01 -0.6492360  1.013879

37.4 GLM例05

在这次练习中，我们来探索几个不同的计数型数据的模型，进一步探讨如何处理过度离散的方法。数据来自Stata的网站，记录的是美国亚利桑那州医院的住院时长数据。使用如下代码下载该数据:

medpar <- haven::read_dta("../Datas/medpar.dta")
head(medpar)
# A tibble: 6 × 10
  provnum  died white   hmo   los age80   age type1 type2 type3
  <chr>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 030001      0     1     0     4     0     4     1     0     0
2 030001      0     1     1     9     0     4     1     0     0
3 030001      1     1     1     3     1     7     1     0     0
4 030001      0     1     0     9     0     6     1     0     0
5 030001      1     1     0     1     1     7     1     0     0
6 030001      1     1     0     4     0     5     1     0     0

我们主要使用的数据是下面这几列：

Variable Description
los length of hospital stay, in days
age Age group
type1 Binary variable indicating elective admission
type2 Binary variable indicating urgent admission
type3 Binary variable indicating emergency admission
  1. 分析住院时间长短和年龄及其他共变量之间关系之前，先了解一下 los 本身的特征。首先，计算 los 的平均值，及其 Wald 法的 95% 置信区间。在怎样的前提条件下，这个置信区间可以被认为有效 (valid)？你认为这些前提条件在这里得到满足了吗？
psych::describe(medpar$los)
   vars    n mean   sd median trimmed  mad min max range skew kurtosis   se
X1    1 1495 9.85 8.83      8    8.61 5.93   1 116   115 3.65    26.16 0.23
t.test(medpar$los)$conf.int
[1]  9.406072 10.302289
attr(,"conf.level")
[1] 0.95

这里默认的前提条件是，住院时长 (days) 服从正态分布。即使住院时长这一数据可能并不100% 服从正态分布，但是如果样本量足够大，那么该 95% 置信区间依然可以被认为近似可以涵盖95%的可重复实验的次数。这一结论依据的是中心极限定理。在这里，住院时长的数据其实分布的非常的偏，并非正态分布。但是我们可以认为由于样本量接近1500，可以认为计算获得的95%置信区间是渐进有效的。

  1. 接下来我们使用 glm 命令来估计 los 的边际均值 (marginal mean)，不加任何预测变量。根据你拟合的泊松回归模型获得的结果，请思考 los 的模型估计 95% 置信区间是多少。和前一步简单的估计相比较，他们是否相似或者有怎样的不同，原因是什么。你认为哪种估计更加有意义？
mA <- glm(los ~ 1, family=poisson(link = log), data=medpar)
jtools::summ(mA, digits = 6, confint = TRUE)
Observations 1495
Dependent variable los
Type Generalized linear model
Family poisson
Link log
χ²(0) -0.000000
Pseudo-R² (Cragg-Uhler) 0.000000
Pseudo-R² (McFadden) 0.000000
AIC 14618.283648
BIC 14623.593529
Est. 2.5% 97.5% z val. p
(Intercept) 2.287896 2.271748 2.304044 277.694425 0.000000
Standard errors: MLE

根据这个模型的结果，住院时长的均值可以被估计为 \(\exp(2.287896) = 9.854\)。但是其95%置信区间的估计是：

下限为，\(\exp(2.287896 - 1.96\times0.008239) = 9.696\)

上限为，\(\exp(2.287896 + 1.96\times0.008239) = 10.014\)

我们发现，均值的点估计，和第一步简单归纳时的结果一致，都是 9.854 天。但是模型估计的95%置信区间 (9.696, 10.014) 相比 Wald 法的 95% 置信区间 (9.406, 10.302) 更加精确 (范围更窄)。当然可以理解为，当数据本身分布较偏 (skew) 时，使用泊松模型分析获得的结果更加可靠且更加有效 (more efficient)。在这个数据中，模型估计的置信区间和wald法置信区间之间的差别更加可能是由于住院时长数据本身的过度离散问题导致的。在R里我们获得的结果只有残差离差量 (residual deviance): Residual deviance: 8901.1 on 1494 degrees of freedom。如果你用的是 Stata，还可以获得 Pearson 统计量，以及他们除以自由度以后获得的数字都显著地大于1。这说明其实病人住院时长这样的数据很大程度上有相当大的差异，因为每个病人各自住院的时间长度更加取决于他们本身患病的程度。 这样的数据不会是通过泊松分布可以简单拟合的，因为泊松分布的均值和方差是严格相等的。 

Iteration 0:   log likelihood = -7354.5568  
Iteration 1:   log likelihood = -7308.2078  
Iteration 2:   log likelihood = -7308.1418  
Iteration 3:   log likelihood = -7308.1418  

Generalized linear models                          No. of obs      =      1495
Optimization     : ML                              Residual df     =      1494
                                                   Scale parameter =         1
Deviance         =  8901.134077                    (1/df) Deviance =  5.957921
Pearson          =  11828.70662                    (1/df) Pearson  =  7.917474

Variance function: V(u) = u                        [Poisson]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  9.778116
Log likelihood   = -7308.141824                    BIC             = -2019.829

------------------------------------------------------------------------------
             |                 OIM
         los |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   2.287896   .0082389   277.69   0.000     2.271748    2.304044
------------------------------------------------------------------------------

也就是下面这两行所提示的内容。

## Deviance         =  8901.134077                   (1/df) Deviance =   5.957921
## Pearson          =  11828.70662                   (1/df) Pearson  =   7.917474
  1. 下面我们在泊松回归模型中加入不同的入院类型的哑变量，看他们是否和患者住院时长有关。尝试用医学文献的文章写法描述这个模型的结果。
mB <- glm(los ~ type2 + type3, family=poisson(link = log), data=medpar)
jtools::summ(mB, digits = 6, confint = TRUE, exp = TRUE)
Observations 1495
Dependent variable los
Type Generalized linear model
Family poisson
Link log
χ²(2) 717.506464
Pseudo-R² (Cragg-Uhler) 0.381200
Pseudo-R² (McFadden) 0.049090
AIC 13904.777184
BIC 13920.706828
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 8.830688 8.659413 9.005350 217.975720 0.000000
type2 1.267877 1.216985 1.320897 11.354985 0.000000
type3 2.065477 1.963233 2.173046 28.003160 0.000000
Standard errors: MLE

记得模型中省略掉了 type1 因为它被当作参考组 (reference group)。

下面的文献描述可以用于参考：

There is strong evidence (p < 0.0001) that the length of hospital admission is related to the type of admission. The mean length of stay for elective admission is 8.83 days (95% CI 8.66 to 9.01 days). Urgent admissions are associated with stays that are on average 26.8% (95% CI 21.7% to 32.1%) longer than those resulting from elective admissions. Emergency admissions result in stays that are on average 2.06 (95% CI 1.96 to 2.17) times as long as those from elective admissions .

值得注意的是，这些结果所依据的泊松回归模型的前提条件很可能因为数据本身存在过度离散 (overdispersion) 的问题而无法得到满足。

  1. 在上述模型中如果加入年龄作为解释变量，试分析年龄是否可以认为是住院时长的独立解释变量 (独立于住院形态 type of admission)。
mC <- glm(los ~ as.factor(age) + type2 + type3, 
          family=poisson(link = log), data=medpar)
jtools::summ(mC, digits = 6, confint = TRUE, exp = TRUE)
Observations 1495
Dependent variable los
Type Generalized linear model
Family poisson
Link log
χ²(10) 735.948666
Pseudo-R² (Cragg-Uhler) 0.388787
Pseudo-R² (McFadden) 0.050351
AIC 13902.334982
BIC 13960.743678
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 11.329126 9.015810 14.236002 20.830249 0.000000
as.factor(age)2 0.775566 0.608964 0.987746 -2.059888 0.039409
as.factor(age)3 0.817641 0.647619 1.032300 -1.692696 0.090513
as.factor(age)4 0.791049 0.628003 0.996425 -1.990375 0.046550
as.factor(age)5 0.765953 0.608122 0.964748 -2.264798 0.023525
as.factor(age)6 0.795170 0.631446 1.001345 -1.948541 0.051350
as.factor(age)7 0.759360 0.601769 0.958220 -2.319574 0.020364
as.factor(age)8 0.723916 0.570628 0.918383 -2.661289 0.007784
as.factor(age)9 0.731655 0.571365 0.936912 -2.476470 0.013269
type2 1.265918 1.214954 1.319021 11.246930 0.000000
type3 2.064878 1.962417 2.172690 27.922756 0.000000
Standard errors: MLE 
lmtest::lrtest(mC, mB)
Likelihood ratio test

Model 1: los ~ as.factor(age) + type2 + type3
Model 2: los ~ type2 + type3
  #Df  LogLik Df  Chisq Pr(>Chisq)  
1  11 -6940.2                       
2   3 -6949.4 -8 18.442    0.01814 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

加入和年龄组作为预测变量的模型结果如上所示。和没有加年龄的模型相比较的似然比检验(likelihood ratio test) 结果显示，如果泊松回归模型前提得到满足，那么有证据证明(p = 0.018)，在调整了住院形态之后，年龄依然是住院时长独立的预测变量。

  1. 重新对加入年龄的泊松回归模型加入稳健统计标准误 (robust standard error) 的估计，获得新的稳健置信区间估计。与非稳健置信区间相比较，你能得出怎样的结论？
jtools::summ(mC, digits = 6, confint = TRUE, exp = TRUE, 
             robust = "HC1")
Observations 1495
Dependent variable los
Type Generalized linear model
Family poisson
Link log
χ²(10) 735.948666
Pseudo-R² (Cragg-Uhler) 0.388787
Pseudo-R² (McFadden) 0.050351
AIC 13902.334982
BIC 13960.743678
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 11.329126 7.637401 16.805336 12.065268 0.000000
as.factor(age)2 0.775566 0.492271 1.221892 -1.095887 0.273128
as.factor(age)3 0.817641 0.539154 1.239974 -0.947602 0.343332
as.factor(age)4 0.791049 0.528283 1.184513 -1.137915 0.255156
as.factor(age)5 0.765953 0.511209 1.147641 -1.292452 0.196201
as.factor(age)6 0.795170 0.529183 1.194851 -1.103144 0.269965
as.factor(age)7 0.759360 0.503126 1.146089 -1.310721 0.189952
as.factor(age)8 0.723916 0.476998 1.098653 -1.517924 0.129034
as.factor(age)9 0.731655 0.469288 1.140705 -1.378954 0.167909
type2 1.265918 1.140543 1.405076 4.431306 0.000009
type3 2.064878 1.640359 2.599262 6.174562 0.000000
Standard errors: Robust, type = HC1

可以看到加入了 robust 选项之后，并不会改变每个变量的回归系数的点估计 (point estimates)。但是，可以发现每个变量的回归系数对应的标准误发生了较大的变化 – 置信区间的范围都无一例外地变大了。由于使用 robust = "HC1" 选项，这里的标准误估计不再依据泊松模型的前提条件 – 泊松分布的特征。在 Stata 中拟合相同的模型，我们可以获得是否有过度离散的指标型数据结果，也就是残差离差量 (residual deviance) 和 Pearson 统计量，以及他们对各自的自由度之比：

Iteration 0:   log pseudolikelihood = -7023.6853  
Iteration 1:   log pseudolikelihood = -6940.3577  
Iteration 2:   log pseudolikelihood = -6940.1675  
Iteration 3:   log pseudolikelihood = -6940.1675  

Generalized linear models                          No. of obs      =      1495
Optimization     : ML                              Residual df     =      1484
                                                   Scale parameter =         1
Deviance         =   8165.18541                    (1/df) Deviance =  5.502147
Pearson          =  9346.752373                    (1/df) Pearson  =  6.298351

Variance function: V(u) = u                        [Poisson]
Link function    : g(u) = ln(u)                    [Log]

                                                   AIC             =  9.299221
Log pseudolikelihood = -6940.167491                BIC             = -2682.679

------------------------------------------------------------------------------
             |               Robust
         los |        IRR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       type2 |   1.265918    .067136     4.45   0.000     1.140942    1.404585
       type3 |   2.064878   .2416633     6.20   0.000     1.641625    2.597257
             |
         age |
          2  |   .7755658   .1792692    -1.10   0.272     .4930223    1.220031
          3  |   .8176414    .173137    -0.95   0.342     .5399076    1.238244
          4  |   .7910486   .1623995    -1.14   0.254     .5289985     1.18291
          5  |   .7659533   .1574871    -1.30   0.195     .5119027    1.146086
          6  |   .7951697   .1646582    -1.11   0.268     .5299061    1.193221
          7  |   .7593596   .1589473    -1.32   0.188     .5038207    1.144508
          8  |   .7239164   .1535639    -1.52   0.128     .4776652    1.097118
          9  |   .7316548   .1652243    -1.38   0.166     .4699868    1.139008
             |
       _cons |   11.32913   2.271633    12.11   0.000     7.647503    16.78314
------------------------------------------------------------------------------

可以看到残差离差量 (residual deviance) 和 Pearson 统计量与各自的自由度之比均显著大于1。提示该数据有相当程度的过度离散，也就是泊松回归模型的前提泊松分布无法得到满足。而使用了 robust = "HC1" (in R) 或者 robust (in Stata) 的选项之后，就不再需要这一前提假设。我们也看到稳健置信区间相较于没有使用该选项时要不那么精确，也就是区间范围都变宽了。

  1. 当考虑了过度离散的模型被采纳后 robust (in Stata)，我们无法再使用似然比检验法检验年龄是否是住院时长的独立预测变量。但是我们可以使用 Stata 里特有的模型拟合之后的 test 命令来实施联合 Wald 检验 (joint Wald test) 年龄是否在稳健泊松模型下依然是住院时长的独立预测变量。试解读该检验结果和之前未考虑过度离散现象时使用的模型的年龄的独立预测变量检验之间有何不同，原因是什么呢？
 ( 1)  [los]2.age - [los]3.age = 0
 ( 2)  [los]2.age - [los]4.age = 0
 ( 3)  [los]2.age - [los]5.age = 0
 ( 4)  [los]2.age - [los]6.age = 0
 ( 5)  [los]2.age - [los]7.age = 0
 ( 6)  [los]2.age - [los]8.age = 0
 ( 7)  [los]2.age - [los]9.age = 0
 ( 8)  [los]2.age = 0

           chi2(  8) =    4.09
         Prob > chi2 =    0.8493

这个联合 Wald 检验的结果是 p = 0.8493。这和之前未考虑数据过度离散时使用的模型时进行的对年龄这一变量的似然比模型检验结果大相径庭 (p = 0.018)。之所以结果相差如此之大，我们相信主要是因为之前忽略数据过度离散问题的模型其实是错误的。而考虑了数据过度离散特征之后，可以认为使用了稳健泊松回归模型之后的联合 Wald 检验结果才是真的值得相信的。

  1. 请使用负二项回归模型 (negative binomial regression model) 来拟合上述模型，先拟合一个不加任何预测变量的负二项回归模型。请解释模型结果的下半部分出现的似然比检验的意义，和无预测变量的泊松回归模型结果做一下对比。

在 Stata 里：

Fitting Poisson model:

Iteration 0:   log likelihood = -7308.1418  
Iteration 1:   log likelihood = -7308.1418  

Fitting constant-only model:

Iteration 0:   log likelihood = -4988.8172  
Iteration 1:   log likelihood = -4857.3372  
Iteration 2:   log likelihood =  -4856.494  
Iteration 3:   log likelihood =  -4856.494  

Fitting full model:

Iteration 0:   log likelihood =  -4856.494  
Iteration 1:   log likelihood =  -4856.494  

Negative binomial regression                      Number of obs   =       1495
                                                  LR chi2(0)      =       0.00
Dispersion     = mean                             Prob > chi2     =          .
Log likelihood = -4856.494                        Pseudo R2       =     0.0000

------------------------------------------------------------------------------
         los |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   2.287896   .0198946   115.00   0.000     2.248903    2.326888
-------------+----------------------------------------------------------------
    /lnalpha |  -.7128727   .0431289                     -.7974038   -.6283417
-------------+----------------------------------------------------------------
       alpha |   .4902339   .0211432                       .450497    .5334758
------------------------------------------------------------------------------
Likelihood-ratio test of alpha=0:  chibar2(01) = 4903.30 Prob>=chibar2 = 0.000

在 R 里：

mD <- MASS::glm.nb(los ~ 1, data=medpar)
jtools::summ(mD, digits = 6, confint = TRUE, exp = TRUE)
Error in glm.control(...) : 
  unused argument (family = list("Negative Binomial(2.0398)", "log", function (mu) 
log(mu), function (eta) 
pmax(exp(eta), .Machine$double.eps), function (mu) 
mu + mu^2/.Theta, function (y, mu, wt) 
2 * wt * (y * log(pmax(1, y)/mu) - (y + .Theta) * log((y + .Theta)/(mu + .Theta))), function (y, n, mu, wt, dev) 
{
    term <- (y + .Theta) * log(mu + .Theta) - y * log(mu) + lgamma(y + 1) - .Theta * log(.Theta) + lgamma(.Theta) - lgamma(.Theta + y)
    2 * sum(term * wt)
}, function (eta) 
pmax(exp(eta), .Machine$double.eps), expression({
    if (any(y < 0)) stop("negative values not allowed for the negative binomial family")
    n <- rep(1, nobs)
    mustart <- y + (y == 0)/6
}), function (mu) 
all(mu > 0), function (eta) 
TRUE, function (object, nsim) 
{
    ftd <- fitted(object)
    rnegbin(nsim * length(ftd), ftd, .Theta)
}))
Warning: Something went wrong when calculating the pseudo R-squared. Returning NA
instead.
Observations 1495
Dependent variable los
Type Generalized linear model
Family Negative Binomial(2.0398)
Link log
χ²(NA) NA
Pseudo-R² (Cragg-Uhler) NA
Pseudo-R² (McFadden) NA
AIC 9716.988008
BIC 9727.607771
exp(Est.) 2.5% 97.5% z val. p
(Intercept) 9.854181 9.477334 10.246011 115.000849 0.000000
Standard errors: MLE 
summary(mD)

Call:
MASS::glm.nb(formula = los ~ 1, data = medpar, init.theta = 2.03984274, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.28790    0.01989     115   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(2.0398) family taken to be 1)

    Null deviance: 1571.5  on 1494  degrees of freedom
Residual deviance: 1571.5  on 1494  degrees of freedom
AIC: 9717

Number of Fisher Scoring iterations: 1

              Theta:  2.0398 
          Std. Err.:  0.0880 

 2 x log-likelihood:  -9712.9880

可以看到在 R 里最下方出现的是 Theta 也就是评价过度离散程度的指标，他和 Stata 的输出报告中的 alpha 互为倒数。不同的是 Stata 的报告中还对 alpha = 0 做了检验。检验结果提示住院时长这个数据并不服从泊松分布。也就是实际的住院时长的数据比泊松分布时的结果的方差要大的多。 (there is more variability than would be expected from a Poisson variable) 其实，简单对 los 分析一下就知道它的样本均值和样本方差分别是：（他们相差巨大，不符合泊松分布的特征）

mean(medpar$los) 
[1] 9.854181
var(medpar$los)
[1] 78.02022
  1. 如果你有兴趣，请拟合一个只有住院形态一个预测变量的负二项回归模型。使用其输出的结果来计算一下不同住院形态下的平均住院时长。比较模型预测的平均住院时长和观测到的不同住院形态下的实际住院时长的平均值之间的差别。

在 Stata 里：

Fitting Poisson model:

Iteration 0:   log likelihood = -6949.6915  
Iteration 1:   log likelihood = -6949.3886  
Iteration 2:   log likelihood = -6949.3886  

Fitting constant-only model:

Iteration 0:   log likelihood = -4988.8172  
Iteration 1:   log likelihood = -4857.3372  
Iteration 2:   log likelihood =  -4856.494  
Iteration 3:   log likelihood =  -4856.494  

Fitting full model:

Iteration 0:   log likelihood = -4802.8473  
Iteration 1:   log likelihood = -4800.2234  
Iteration 2:   log likelihood = -4800.2189  
Iteration 3:   log likelihood = -4800.2189  

Negative binomial regression                      Number of obs   =       1495
                                                  LR chi2(2)      =     112.55
Dispersion     = mean                             Prob > chi2     =     0.0000
Log likelihood = -4800.2189                       Pseudo R2       =     0.0116

------------------------------------------------------------------------------
         los |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       type2 |   .2373439   .0502166     4.73   0.000     .1389211    .3357666
       type3 |   .7253612   .0757014     9.58   0.000     .5769893    .8737332
       _cons |   2.178233   .0222433    97.93   0.000     2.134637    2.221829
-------------+----------------------------------------------------------------
    /lnalpha |  -.8033559   .0443828                     -.8903446   -.7163671
-------------+----------------------------------------------------------------
       alpha |   .4478236   .0198757                      .4105142    .4885238
------------------------------------------------------------------------------
Likelihood-ratio test of alpha=0:  chibar2(01) = 4298.34 Prob>=chibar2 = 0.000

在 R 里：

mE <- MASS::glm.nb(los ~ type2 + type3, data=medpar)
jtools::summ(mE, digits = 6, confint = TRUE)
Error in glm.control(...) : 
  unused argument (family = list("Negative Binomial(2.233)", "log", function (mu) 
log(mu), function (eta) 
pmax(exp(eta), .Machine$double.eps), function (mu) 
mu + mu^2/.Theta, function (y, mu, wt) 
2 * wt * (y * log(pmax(1, y)/mu) - (y + .Theta) * log((y + .Theta)/(mu + .Theta))), function (y, n, mu, wt, dev) 
{
    term <- (y + .Theta) * log(mu + .Theta) - y * log(mu) + lgamma(y + 1) - .Theta * log(.Theta) + lgamma(.Theta) - lgamma(.Theta + y)
    2 * sum(term * wt)
}, function (eta) 
pmax(exp(eta), .Machine$double.eps), expression({
    if (any(y < 0)) stop("negative values not allowed for the negative binomial family")
    n <- rep(1, nobs)
    mustart <- y + (y == 0)/6
}), function (mu) 
all(mu > 0), function (eta) 
TRUE, function (object, nsim) 
{
    ftd <- fitted(object)
    rnegbin(nsim * length(ftd), ftd, .Theta)
}))
Warning: Something went wrong when calculating the pseudo R-squared. Returning NA
instead.
Observations 1495
Dependent variable los
Type Generalized linear model
Family Negative Binomial(2.233)
Link log
χ²(NA) NA
Pseudo-R² (Cragg-Uhler) NA
Pseudo-R² (McFadden) NA
AIC 9608.437779
BIC 9629.677305
Est. 2.5% 97.5% z val. p
(Intercept) 2.178233 2.134637 2.221829 97.927402 0.000000
type2 0.237344 0.138921 0.335767 4.726402 0.000002
type3 0.725361 0.576989 0.873733 9.581877 0.000000
Standard errors: MLE 
summary(mE)

Call:
MASS::glm.nb(formula = los ~ type2 + type3, data = medpar, init.theta = 2.233022138, 
    link = log)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  2.17823    0.02224  97.927  < 2e-16 ***
type2        0.23734    0.05022   4.726 2.29e-06 ***
type3        0.72536    0.07570   9.582  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for Negative Binomial(2.233) family taken to be 1)

    Null deviance: 1685.1  on 1494  degrees of freedom
Residual deviance: 1568.1  on 1492  degrees of freedom
AIC: 9608.4

Number of Fisher Scoring iterations: 1

              Theta:  2.2330 
          Std. Err.:  0.0991 

 2 x log-likelihood:  -9600.4380

从输出的分析报告结果来看，模型估计的不同住院形态(elective, urgent, and emergency) 的平均住院时长分别是\(\exp(2.1782) = 8.83\) 天，\(\exp(2.1782 + 0.2373) = 11.20\) 天，\(\exp(2.1782 + 0.7253) = 18.24\)天。要估算模型的预测方差，可以通过手工计算。已知如果期待值是 \(\mu\)，那么其方差是 \(\mu(1 + \alpha \mu)\)。那么依据这个公式，就可以估算不同住院形态(elective, urgent, and emergency) 的住院时长的方差分别是\(8.83 \times (1 + 0.4478 \times 8.83) = 43.74\) days\(^2\)，\(11.2 \times (1 + 0.4478 \times 11.2) = 67.3\) days\(^2\)，\(18.24 \times (1 + 0.4478 \times 18.24) = 167.2\) days\(^2\)。

和观测值相比较：

var(medpar$los[medpar$type1 == 1])
[1] 41.68005
var(medpar$los[medpar$type2 == 1])
[1] 77.87802
var(medpar$los[medpar$type3 == 1])
[1] 424.8788

可见，对于前两种住院形态来说，模型推测的方差和观测方差还是比较接近的。但是对于第三种住院形态 (emergency)，观测方差远远大于模型推测方差。当然，这可能由于因为紧急住院的患者人数较低 (n = 96)。也就是说，样本量太低的组使用该模型推测的方差准确度就比较低。另一种原因值得考虑的就是，负二项回归模型中使用的Gamma 随机效应分布可能并不适用与本数据(对于不同住院形态的住院时长可能需要不同的alpha，而不是强迫他们都是相同的)。另外一个被忽视的是，住院时长这一数据其实是只能取正值的数据。然而泊松分布和负二项分布本身其实是允许有 0 这样的数字的。如此一来，或许我们应该把住院时长定义为允许 0 存在的数字。也就是把住院时长的数据减去1。 (number of days from the day of admission)