第 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 dayssummary(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.00879mhgrade_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.7611078mhgrade_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.