# 第 47 章 計數型因變量 Poisson regression

1. 某個呼吸科診所的患者中，每個人在過去一個月中哮喘發作的次數；
2. 癲癇患者在過去一年中癲癇發作次數；
3. 接受腦部 CT 掃描的患者中，每個人被診斷出顱內腫瘤個數。

## 47.1 泊松 GLM

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

\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}

\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)。

## 47.2 泊松迴歸實例

p <- read_csv("backupfiles/poisson_sim.csv")
## Parsed with column specification:
## cols(
##   id = col_double(),
##   num_awards = col_double(),
##   prog = col_double(),
##   math = col_double()
## )
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.000
##  2      :  1   1st Qu.:0.00   Academic  :105   1st Qu.:45.000
##  3      :  1   Median :0.00   Vocational: 50   Median :52.000
##  4      :  1   Mean   :0.63                    Mean   :52.645
##  5      :  1   3rd Qu.:1.00                    3rd Qu.:59.000
##  6      :  1   Max.   :6.00                    Max.   :75.000
##  (Other):194

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)
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -1.41421  -0.69282  -0.63246   0.00000   3.39133
##
## Coefficients:
##                Estimate Std. Error z value  Pr(>|z|)
## (Intercept)    -1.60944    0.33333 -4.8283 1.377e-06 ***
## progAcademic    1.60944    0.34733  4.6338 3.590e-06 ***
## progVocational  0.18232    0.44096  0.4135    0.6793
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 287.672  on 199  degrees of freedom
## Residual deviance: 234.460  on 197  degrees of freedom
## AIC: 416.515
##
## Number of Fisher Scoring iterations: 6
##
## Call:
## glm(formula = num_awards ~ prog, family = poisson(link = log),
##     data = p)
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -1.41421  -0.69282  -0.63246   0.00000   3.39133
##
## Coefficients:
##                Estimate Std. Error z value  Pr(>|z|)
## (Intercept)    -1.60944    0.33333 -4.8283 1.377e-06 ***
## progAcademic    1.60944    0.34733  4.6338 3.590e-06 ***
## progVocational  0.18232    0.44096  0.4135    0.6793
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 287.672  on 199  degrees of freedom
## Residual deviance: 234.460  on 197  degrees of freedom
## AIC: 416.515
##
## 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
## progVocational 1.2  0.51   2.94

## 47.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

$E(Y_i) = E(E(Y_i | \mu_i)) = E(\mu_i) = \mu$

\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}

### 47.3.1 過度離散怎麼查？

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.45997 197 0.034961171

Goodness of fit 檢驗結果 提示本模型可能存在過度離散，數據擬合度不理想。值得注意的是如果樣本很大時，模型偏差的檢驗統計量將不再服從卡方分佈，應用的時候一定要慎重。

### 47.3.2 負二項式分佈模型 negative binomial model

R 裏擬合負二項式分佈模型的函數 glm.nb 在基本包 MASS 裏。

m1 <- glm.nb(num_awards ~ prog, data = p)
m2 <- glm(num_awards ~ prog, family=poisson(link = log), data=p)
summary(m1)
##
## Call:
## glm.nb(formula = num_awards ~ prog, data = p, init.theta = 1.72267107,
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -1.25581  -0.67036  -0.61517   0.00000   2.32349
##
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -1.60944    0.35215 -4.5703 4.87e-06 ***
## progAcademic    1.60944    0.37291  4.3159 1.59e-05 ***
## progVocational  0.18232    0.46793  0.3896   0.6968
## ---
## 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.264  on 199  degrees of freedom
## Residual deviance: 171.066  on 197  degrees of freedom
## AIC: 406.532
##
## 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)
##
## Deviance Residuals:
##      Min        1Q    Median        3Q       Max
## -1.41421  -0.69282  -0.63246   0.00000   3.39133
##
## Coefficients:
##                Estimate Std. Error z value  Pr(>|z|)
## (Intercept)    -1.60944    0.33333 -4.8283 1.377e-06 ***
## progAcademic    1.60944    0.34733  4.6338 3.590e-06 ***
## progVocational  0.18232    0.44096  0.4135    0.6793
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
##     Null deviance: 287.672  on 199  degrees of freedom
## Residual deviance: 234.460  on 197  degrees of freedom
## AIC: 416.515
##
## 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.06608 197 0.90901874
with(m2, cbind(res.deviance = deviance, df = df.residual,
p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance  df           p
## [1,]    234.45997 197 0.034961171

m2 <- glm(num_awards ~ prog, family=poisson(link = log), data=p)
cov.m2 <- 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.60943791 0.29814240 6.7305731e-08 -2.19379701 -1.0250788
## progAcademic    1.60943791 0.32296809 6.2517920e-07  0.97642045  2.2424554
## progVocational  0.18232156 0.42426407 6.6738767e-01 -0.64923602  1.0138791