I like the indirect standardisation function in epiR. I’ve made a wrapper that lets you put in data rather than fiddly summary tables.

‘dat’ is a data.frame containing the other variables.

‘outcome’ is a character name of the binary outcome variable.

‘strata’ is a character name of the factor variable that you want to stratify by (e.g. region).

‘agegrp’ is a character name of the factor variable that you want to standardise for (e.g. age group).

‘survival’ is a character name of the numeric variable showing survival time. If you include this variabe, the output will be a rate ratio; otherwise it will be a risk ratio.

‘ref_strat’ is a character name of the stratum that you want to use as the reference value. If you include this variable, the specified stratum will have an SMR of 100. Otherwise, the study population will be used as a reference.

SMR.wrapper <- function(dat, outcome, strata, agegrp, survival = NA, ref_strat = NA) { labs <- levels(dat[,strata]) levs <- 1:nlevels(dat[,agegrp]) dat$ind <- 1 strat_deaths <- matrix(aggregate(dat[,outcome], by = list(dat[,strata]), FUN = sum)[,2], dimnames = list(labs, "")) mths_at_risk <- matrix(xtabs(ifelse(rep(is.na(survival), nrow(dat)), dat$ind, dat[,survival]) ~ dat[,strata] + dat[,agegrp]), nlevels(dat[,strata]), nlevels(dat[,agegrp]), dimnames = list(labs, levs)) if(!is.na(ref_strat)) { dat <- dat[dat[,strata] == ref_strat,] } ref_deaths <- aggregate(dat[,outcome], by = list(dat[,agegrp]), FUN = sum) ref_deaths <- c(ref_deaths[,2], sum(ref_deaths[,2])) ref_mths <- aggregate(ifelse(rep(is.na(survival), nrow(dat)), dat$ind, dat[,survival]), by = list(dat[,agegrp]), FUN = sum) ref_mths <- c(ref_mths[,2], sum(ref_mths[,2])) ref_rate <- matrix(ref_deaths / ref_mths, 1, dimnames = list("", c(levs, "Total"))) return(epi.indirectadj(strat_deaths, mths_at_risk, ref_rate, units = 100)) }