第 96 章 Crude and stratified rate ratios

本章內容不討論任何理論的東西,着重強調用 R 進行實際數據的分析,並加強對輸出結果的理解。

此次實戰演練的目的是學會怎樣計算死亡率比 (Rate Ratios, RR)。學會用 Mantel-Haenszel 法總結 RR,並討論其意義。

whitehal <- read_dta("backupfiles/whitehal.dta")
whitehal$followupyrs <- (whitehal$timeout - whitehal$timein)/365.25
max(whitehal$followupyrs*365.25) # time difference in days
## Time difference of 7078.9927 days
summary(whitehal$followupyrs <- as.numeric(whitehal$followupyrs)) # time difference in years
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.15063 17.16357 17.96020 16.46116 18.56262 19.38123
# categorize agein into groups (40-44, 45-49, 50-54, ... , 65-69)
whitehal$agecat <- cut(whitehal$agein, breaks = seq(40, 70, 5), right = FALSE)
with(whitehal, table(agecat))
## agecat
## [40,45) [45,50) [50,55) [55,60) [60,65) [65,70) 
##     277     445     362     340     215      38
# examine how mortality rates change with age at entry
# 
# with(whitehal %>% group_by(agecat) %>%
#   summarise(D = sum(all),
#             Y = sum(followupyrs)),
#   cbind(whitehal$agecat, pois.exact(x = D, pt = Y/1000)))


## rate ratios and 95% CIs for each age category compare with [40,44) age group
Model0 <- glm(all ~ agecat + offset(log(followupyrs)), family = poisson(link = "log"), data = whitehal); ci.exp(Model0)
##                   exp(Est.)          2.5%         97.5%
## (Intercept)    0.0048794197  0.0032705278  0.0072797841
## agecat[45,50)  1.1740971320  0.7154062934  1.9268827913
## agecat[50,55)  2.7732622719  1.7597191900  4.3705743919
## agecat[55,60)  4.5784055712  2.9519678778  7.1009572062
## agecat[60,65)  6.6882727538  4.2856745766 10.4377949444
## agecat[65,70) 17.1110624279 10.1140257533 28.9487553770
## The rate ratios are increasing with age although there is no statistical evidence
## at 5% level that the rate among 45-49 year olds is different to the rate among men
## who are <40 years


# with(whitehal %>% group_by(grade) %>%
#   summarise(D = sum(all),
#             Y = sum(followupyrs)),
#   cbind(whitehal$grade, pois.exact(x = D, pt = Y/1000)))

Model1 <- glm(all ~ factor(grade) + offset(log(followupyrs)), family = poisson(link = "log"), data = whitehal); ci.exp(Model1)
##                  exp(Est.)         2.5%       97.5%
## (Intercept)    0.010865404 0.0095233128 0.012396632
## factor(grade)2 2.305446287 1.8947528380 2.805158793
## There is strong evidence that the all cause mortality rate differs between high
## and low grade workers.

## To examine whether the estimated RR for grade is confounded by age at entry
## we compare the crude RR =2.31 (1.90, 2.81) with the Mantel-Haenszel summary
## estimate.

whitehal_table <- aggregate(cbind(all, followupyrs) ~ grade + agecat, data=whitehal, sum)
stmh_array <- array(c(4, 20,   693.1284,4225.4893,
                      10,35,   1363.821,6491.072,
                      30,52,  1399.63, 4660.12,                                                        51,67,   1832.169,3449.846,
                      59,42,   1660.597,1434.251,
                      28,5,  316.23840, 79.00879),
                      dim=c(2,2,6),
                      dimnames = list(
                      Grade=c("2","1"),
                      c("death", "Person_years"),
                      Agecat=names(table(whitehal$agecat))
                    ))
stmh_array
## , , Agecat = [40,45)
## 
##      
## Grade death Person_years
##     2     4     693.1284
##     1    20    4225.4893
## 
## , , Agecat = [45,50)
## 
##      
## Grade death Person_years
##     2    10     1363.821
##     1    35     6491.072
## 
## , , Agecat = [50,55)
## 
##      
## Grade death Person_years
##     2    30      1399.63
##     1    52      4660.12
## 
## , , Agecat = [55,60)
## 
##      
## Grade death Person_years
##     2    51     1832.169
##     1    67     3449.846
## 
## , , Agecat = [60,65)
## 
##      
## Grade death Person_years
##     2    59     1660.597
##     1    42     1434.251
## 
## , , Agecat = [65,70)
## 
##      
## Grade death Person_years
##     2    28    316.23840
##     1     5     79.00879
mhgrade_age <- epi.2by2(stmh_array, method = "cohort.time", units = 1000)
mhgrade_age
##              Outcome +    Time at risk        Inc rate *
## Exposed +          182            7266              25.0
## Exposed -          221           20340              10.9
## Total              403           27605              14.6
## 
## Point estimates and 95 % CIs:
## -------------------------------------------------------------------
## Inc rate ratio (crude)                       2.31 (1.88, 2.82)
## Inc rate ratio (M-H)                         1.43 (1.16, 1.76)
## Inc rate ratio (crude:M-H)                   1.61
## Attrib rate (crude) *                        14.18 (10.27, 18.10)
## Attrib rate (M-H) *                          6.28 (2.42, 10.15)
## Attrib rate (crude:M-H)                      2.26
## -------------------------------------------------------------------
##  Wald confidence limits
##  M-H: Mantel-Haenszel
##  * Outcomes per 1000 units of population time at risk
## Overall estimate and Wald 95% confidence intervals,
## controlling for agecate
mhgrade_age$massoc$IRR.mh.wald
##        est     lower     upper
## 1 1.429211 1.1598631 1.7611078
mhgrade_age$massoc$chisq.mh ## p-value for age-adjusted MH rate ratio
##   test.statistic df      p.value
## 1      10.760807  1 0.0010367217
## The Mantel-Haenszel summary estimate RR = 1.43 (1.16, 1.76).
## The result shows that the crude estimate of the effect of grade was
## partly confounded by age at entry.

## To assess whether there is effect modification betwee grade and
## agecat we examine the stratum specific estimates and assess
## whether there is evidence of important variation between them.
mhgrade_age$massoc$IRR.strata.wald
##         est      lower     upper
## 1 1.2192515 0.30302945 3.6397113
## 2 1.3598500 0.60057144 2.8059055
## 3 1.9208868 1.18309546 3.0676307
## 4 1.4332751 0.97576189 2.0941999
## 5 1.2132872 0.80308583 1.8474726
## 6 1.3991002 0.53338106 4.6404658
## The result indicates that the data are compatible with the assumption
## of no interaction/effect modification (p=0.79)

## test for unequal RRs (effect modification):
mhgrade_age$res$RR.homog
##   test.statistic df    p.value
## 1       2.407279  5 0.79038965
## Hence, we do not need to present the stratum-specific estimates.