# 第 114 章 粗率比和分層率比 Crude and stratified rate ratios

Epidemiology is like a bikini: what is revealed is interesting; what is concealed is crucial.
Peter Duesberg
The Causal Inference lectures were orgainised and taught by Professor Simon Cousens, Professor Katherine Fielding and Dr Susannah Woodd.

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

## 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(
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; CI: confidence interval ## * 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 11.132745 5 0.00084816981 ## 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
## NULL
## Hence, we do not need to present the stratum-specific estimates.