I've got survival data with an outcome, an exposure, and a variable I'd like to stratify on but can't find a function to perform Mantel-Haenszel rate ratio. For example with the lung dataset from survival I'd like to look at the outcome of status based on sex but stratified by age. I've set up age brackets with
library(tidyverse)
library(survival)
lung2 <- lung %>%
mutate(agecat = as.factor(case_when(age < 50 ~ 0,
age < 70 ~ 1,
age >= 70 ~ 2)))
epi.2by2 from epiR gets me close with
library(epiR)
epi.2by2(table(as.factor(lung2$status),
as.factor(lung2$sex),
lung2$agecat),
method = "cohort.count")
#> Outcome + Outcome - Total Inc risk * Odds
#> Exposed + 26 37 63 41.3 0.703
#> Exposed - 112 53 165 67.9 2.113
#> Total 138 90 228 60.5 1.533
#>
#>
#> Point estimates and 95% CIs:
#> -------------------------------------------------------------------
#> Inc risk ratio (crude) 0.61 (0.44, 0.83)
#> Inc risk ratio (M-H) 0.61 (0.45, 0.84)
#> Inc risk ratio (crude:M-H) 0.99
#> Odds ratio (crude) 0.33 (0.18, 0.61)
#> Odds ratio (M-H) 0.34 (0.19, 0.63)
#> Odds ratio (crude:M-H) 0.97
#> Attrib risk in the exposed (crude) * -26.61 (-40.70, -12.52)
#> Attrib risk in the exposed (M-H) * -25.88 (-42.58, -9.19)
#> Attrib risk (crude:M-H) 1.03
#> -------------------------------------------------------------------
#> M-H test of homogeneity of PRs: chi2(2) = 0.191 Pr>chi2 = 0.909
#> M-H test of homogeneity of ORs: chi2(2) = 0.394 Pr>chi2 = 0.821
#> Test that M-H adjusted OR = 1: chi2(1) = 12.299 Pr>chi2 = <0.001
#> Wald confidence limits
#> M-H: Mantel-Haenszel; CI: confidence interval
#> * Outcomes per 100 population units
But it doesn't take time to event data (rates) into consideration. It has a method = "cohort.time" option, but I can't seem to get it to work. Ultimately the output I'd like would be similar to STATAs stmh sex, by(agecat) which would give risk ratio estimates for each strata with upper and lower 95% confidence interval as well as overall estimate for risk ratio with chi-square and p-value. mhor from epiDisplay gives output close to what I'm looking for
library(epiDisplay)
mhor(lung2$status,
lung2$sex,
lung2$agecat, design = "cohort",
graph = FALSE)
#>
#> Stratified analysis by Var3
#> OR lower lim. upper lim. P value
#> Var3 0 0.476 0.0532 3.786 0.653417
#> Var3 1 0.363 0.1642 0.788 0.006387
#> Var3 2 0.243 0.0424 1.222 0.060217
#> M-H combined 0.344 0.1880 0.631 0.000453
#>
#> M-H Chi2(1) = 12.3 , P value = 0
#> Homogeneity test, chi-squared 2 d.f. = 0.38 , P value = 0.828
but it only gives odds ratios and not rate ratios.
