第 48 章 率的廣義線性迴歸 Poisson GLM for rates

48.1 醫學中的率

前章介紹的事件發生次數，使用的是泊松迴歸。本章介紹同樣利用泊松迴歸，對事件發生率類型數據的泊松迴歸模型。常見的率的數據例如：

  • 肺癌發病率
  • 工廠職工的死亡率
  • 術後後遺症的發生率

下列數據來自英國醫生調查 (British doctors study)，研究的是男性醫生中吸菸與否和冠心病死亡之間的關係。最後一列是每組觀測對象被追蹤的人年 (person-year)。

##    agegrp     smokes deaths  pyrs
## 1   35-44     Smoker     32 52407
## 2   45-54     Smoker    104 43248
## 3   55-64     Smoker    206 28612
## 4   65-74     Smoker    186 12663
## 5     75+     Smoker    102  5317
## 6   35-44 Non-smoker      2 18790
## 7   45-54 Non-smoker     12 10673
## 8   55-64 Non-smoker     28  5710
## 9   65-74 Non-smoker     28  2585
## 10    75+ Non-smoker     31  1462

這是一個已經被整理過的數據，我們沒有辦法從這樣的數據還原到每個觀察對象的個人水平數據。冠心病的粗死亡率 (crude death rate) 可以被計算如下表 (忽略年齡分組)，此時默認的前提是死亡事件在追蹤的過程中發生的概率不發生改變。

表 48.1: Death rates due to CHD in smokers and non-smokers, collapsed over age group
Group Person-years of follow-up CHD deaths Death Rate per 1000 person-years Rate Ratios
Non-smokers 39220 101 2.58 1.00
Smokers 142247 630 4.43 1.72

48.2 泊松過程

設 \(Y\) 是代表某段時間 \(t\) 內事件發生次數 (死亡) 的隨機變量。如果可以假設：

  • 每次事件的發生，是互相獨立的，即在沒有重疊的時間線上，每個事件的發生是隨機的。
  • 在一個無限小的時間段 \(\delta t\) 內，事件發生的概率是 \(\lambda\times\delta t\)，其中 \(\delta t \rightarrow 0\)。

那麼根據泊松分佈 (Section 6) 的定義，在這個時間段內，隨機變量 \(Y\) 事件發生次數服從泊松分佈：

\[ \begin{aligned} Y & \sim \text{Po}(\mu) \\ \text{Where } \mu & = \lambda t, \text{ and } \lambda \text{ is the Rate} \end{aligned} \]

所以，從泊松過程可以看到，我們關心的參數是事件發生率 \(\lambda\)。

48.3 率的模型

既然關心的參數只是發生率，且我們已知泊松分佈是指數分佈的家族成員，可以用廣義線性模型的概念來建模。

  1. 因變量分佈，distribution of dependent variable \[Y_i \sim \text{Po}(\mu_i), \text{ where } \mu_i = \lambda_i t_i\]
  2. 線性預測方程，linear predictor \[\eta_i = \alpha + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}\]
  3. 標準鏈接方程，canonical link function \[\text{log}(\lambda_i) = \text{log}(\frac{\mu_i}{t_i})\]

所以，將率的模型整理一下，就變成了

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

你可以看到，時間項的對數部分 \(\text{log}(t_i)\) 其實是被移到線性預測方程的右邊跟參數放在一起的，只是它的迴歸係數被強制爲 \(1\)。這個時間項被叫做 補償項 (offset)。這樣我們就成功地擬合了用於求事件發生率的一個泊松迴歸模型。在 R 裏，你可以用 glm() 命令的 offset = 選項功能，也可以把 offset(log(Person-year)) 作爲線性預測方程的一部分把時間項取對數以後放進模型裏面。

48.4 率的 GLM

所以我們一起來把率的 GLM 正式定義一下，它包含三個部分：

  1. 可被認爲互相獨立的因變量觀測值的分佈服從泊松分佈 \[Y_i \sim \text{Po}(\mu_i)\]
    其中 \(E(Y_i) = \mu_i = \lambda_i t_i\)，\(t_i\) 是第 \(i\) 個觀察對象 (或者觀察組) 的追蹤人年 (person-time)。
  2. 線性預測方程 \[\eta_i = \text{log}(t_i) + \alpha + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}\]
  3. 鏈接方程是均值的對數方程 \[\text{log}(\mu_i) = \eta_i\]

和分組型二項分佈數據相似，如果泊松 GLM 擬合的數據也是分組型數據，如本章開頭的英國醫生隊列數據。那麼模型偏差值 (deviance) 可以用來衡量模型擬合的好壞。在零假設條件下，模型偏差值服從自由度爲 \(n-p\) 的卡方分佈 (這裏的 \(n\) 是分組型數據中的“組的數量”，也就是飽和模型中參數的數量，\(p\) 是擬合的線性預測方程中參數的數量)。

48.5 實戰演練

數據是本章開頭使用的英國醫生隊列

##    agegrp     smokes deaths  pyrs
## 1   35-44     Smoker     32 52407
## 2   45-54     Smoker    104 43248
## 3   55-64     Smoker    206 28612
## 4   65-74     Smoker    186 12663
## 5     75+     Smoker    102  5317
## 6   35-44 Non-smoker      2 18790
## 7   45-54 Non-smoker     12 10673
## 8   55-64 Non-smoker     28  5710
## 9   65-74 Non-smoker     28  2585
## 10    75+ Non-smoker     31  1462
  • 每組的死亡人數用 \(y_i, i=1,\cdots,10\) 標記；
  • 每組追蹤的人年用 \(t_i\) 標記；
  • \(x_{i1} = 0\) 時對象是吸菸者，\(x_{i1} = 1\) 時對象是非吸菸者；
  • \(x_{i2}, x_{i3}, x_{i4}, x_{i5}\) 作爲5個年齡組的啞變量。

分析目的是：

  1. 調查吸菸與冠心病死亡率的關係 (不調整年齡)；
  2. 調查吸菸與冠心病死亡率的年齡調整後關係；
  3. 調查年齡是否對吸菸與冠心病死亡率的關係起到交互作用。

48.5.1 模型 1

第一個模型可以用下面的數學表達式：

\[ \text{log}(\mu_i) = \text{log}(t_i) + \alpha + \beta_1 x_{i1} \]

在 R 裏面用下面的代碼來擬合這個模型，仔細閱讀輸出的結果：

# the following 2 models are equivalent
Model1 <- glm(deaths ~ smokes + offset(log(pyrs)), family = poisson(link = "log"), data = BritishD)
Model1 <- glm(deaths ~ smokes, offset = log(pyrs), family = poisson(link = "log"), data = BritishD)
summary(Model1)
## 
## Call:
## glm(formula = deaths ~ smokes, family = poisson(link = "log"), 
##     data = BritishD, offset = log(pyrs))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -16.5348   -6.0313    4.6116    8.1617   13.6441  
## 
## Coefficients:
##               Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept)  -5.961822   0.099504 -59.9157 < 2.2e-16 ***
## smokesSmoker  0.542221   0.107183   5.0588 4.219e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 935.067  on 9  degrees of freedom
## Residual deviance: 905.976  on 8  degrees of freedom
## AIC: 965.044
## 
## Number of Fisher Scoring iterations: 6

輸出報告中的參數估計部分 Estimate 就是我們擬合模型中參數的估計 \(\hat\alpha, \hat\beta_1\)，他們各自的含義是：

  • \(\hat\alpha = -5.96\)：非吸菸者的冠心病估計死亡率的對數 (the estimated log rate for non-smokers)；
  • \(\hat\beta_1 = 0.547\)：非吸菸者和吸菸者兩組之間冠心病死亡率對數之差 (the estimated difference in log rate between non-smokers and smokers)。

注意看報告中間部分模型偏差部分的數字 Residual deviance: 905.98 on 8 degrees of freedom，如果對 模型 1 進行擬合優度檢驗：

with(Model1, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance df              p
## [1,]    905.97619  8 2.9024145e-190

擬合優度檢驗結果提示，這個模型對數據的擬合非常差 (poor fit)。可能的原因是，模型 1 中忽略了“年齡”這一重要的因素，使得當僅僅使用 吸菸與否 的信息擬合的泊松迴歸模型的擬合值和觀察值之間的差異的波動非常大，大到很可能無法滿足泊松分佈的前提假設。

48.5.2 模型 2

第二個模型的線性預測方程可以寫作：

\[ \text{log}(\mu_i) = \text{ln}(t_i) + \alpha + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + \beta_4 x_{i4} + \beta_5 x_{i5} \]

在 R 裏面用下面的代碼來擬合這個模型，仔細閱讀輸出的結果：

Model2 <- glm(deaths ~ smokes + agegrp + offset(log(pyrs)), family = poisson(link = "log"), data = BritishD)
summary(Model2)
## 
## Call:
## glm(formula = deaths ~ smokes + agegrp + offset(log(pyrs)), family = poisson(link = "log"), 
##     data = BritishD)
## 
## Deviance Residuals: 
##         1          2          3          4          5          6          7          8          9  
##  0.901600   0.510379   0.051347  -0.087318  -0.912369  -2.179780  -1.308000  -0.137907   0.228819  
##        10  
##  1.919020  
## 
## Coefficients:
##              Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept)  -7.91933    0.19176 -41.2978 < 2.2e-16 ***
## smokesSmoker  0.35454    0.10737   3.3019 0.0009604 ***
## agegrp45-54   1.48401    0.19510   7.6063 2.821e-14 ***
## agegrp55-64   2.62751    0.18373  14.3012 < 2.2e-16 ***
## agegrp65-74   3.35049    0.18480  18.1305 < 2.2e-16 ***
## agegrp75+     3.70010    0.19222  19.2494 < 2.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: 935.0673  on 9  degrees of freedom
## Residual deviance:  12.1324  on 4  degrees of freedom
## AIC: 79.2003
## 
## Number of Fisher Scoring iterations: 4

此時可以計算吸菸者與非吸菸者相比時，年齡調整後冠心病死亡率的比爲：

\[ \begin{aligned} e^{0.3545} & = 1.43 \text{ with } 95\% \text{ CI: } \\ (e^{0.3545 - 1.96\times0.1074}, & e^{0.3545 + 1.96\times0.1074}) = (1.16, 1.76) \end{aligned} \]

報告中還包含了對吸菸項迴歸係數的 Wald 檢驗結果 smokesSmoker 0.3545 0.1074 3.302 0.00096 ***，從這一結果來看，數據提供了強有力的證據證明了年齡調整以後，吸菸會引起冠心病死亡率的顯著升高。再利用模型擬合報告中模型偏差部分的數據 Residual deviance: 905.98 on 8 degrees of freedom，模型的擬合優度檢驗結果爲：

with(Model2, cbind(res.deviance = deviance, df = df.residual,
  p = pchisq(deviance, df.residual, lower.tail=FALSE)))
##      res.deviance df           p
## [1,]    12.132366  4 0.016393625

結果依然提示，即使把年齡組放入這個泊松迴歸，模型對數據的擬合程度依然非常的不好。所以，到這裏，在即使調整了年齡之後模型擬合度依然不理想的情況下 (這是需要加交互作用項的證據)，我們需要在模型中加入年齡和吸菸的交互作用項 (結果是加入交互作用項的模型就變成了飽和模型)。

48.5.3 模型 3

Model3 <- glm(deaths ~ smokes*agegrp + offset(log(pyrs)), family = poisson(link = "log"), data = BritishD)
summary(Model3)
## 
## Call:
## glm(formula = deaths ~ smokes * agegrp + offset(log(pyrs)), family = poisson(link = "log"), 
##     data = BritishD)
## 
## Deviance Residuals: 
##  [1]  0  0  0  0  0  0  0  0  0  0
## 
## Coefficients:
##                          Estimate Std. Error  z value  Pr(>|z|)    
## (Intercept)              -9.14793    0.70711 -12.9371 < 2.2e-16 ***
## smokesSmoker              1.74687    0.72887   2.3967  0.016544 *  
## agegrp45-54               2.35737    0.76376   3.0865  0.002025 ** 
## agegrp55-64               3.83016    0.73192   5.2330 1.668e-07 ***
## agegrp65-74               4.62266    0.73192   6.3158 2.688e-10 ***
## agegrp75+                 5.29436    0.72956   7.2569 3.960e-13 ***
## smokesSmoker:agegrp45-54 -0.98662    0.79006  -1.2488  0.211741    
## smokesSmoker:agegrp55-64 -1.36281    0.75619  -1.8022  0.071512 .  
## smokesSmoker:agegrp65-74 -1.44229    0.75653  -1.9065  0.056592 .  
## smokesSmoker:agegrp75+   -1.84699    0.75717  -2.4393  0.014715 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 9.35067e+02  on 9  degrees of freedom
## Residual deviance: 4.44089e-15  on 0  degrees of freedom
## AIC: 75.0679
## 
## Number of Fisher Scoring iterations: 3

此時你會看到模型的偏差已經幾乎接近於零，因爲這已經是一個飽和模型。