Biostatistics III in R

Exercise 23. Calculating SMRs/SIRs

---

You may have to install the required packages the first time you use them. You can install a package by install.packages("package_of_interest") for each package you require.

library(biostat3) # for Surv and survfit
library(dplyr)    # for data manipulation

Load the melanoma data. Restrict the data to the localised status of the melanoma and to the 10 firts years of follow-up. Use the time-on-study as the timescale. define the event as death from cancer.

The standardized mortality ratio (SMR) is the ratio of the observed number of deaths in the study population to the number that would be expected if the study population experienced the same mortality as the standard population. It is an indirectly standardized rate. When studying disease incidence the corresponding quantity is called a standardized incidence ratio (SIR). These measures are typically used when the entire study population is considered ‘exposed’. Rather than following-up both the exposed study population and an unexposed control population and comparing the two estimated rates we instead only estimate the rate (or number of events) in the study population and compare this to the expected rate (expected number of events) for the standard population. For example, we might study disease incidence or mortality among individuals with a certain occupation (farmers, painters, airline cabin crew) or cancer incidence in a cohort exposed to ionising radiation.

In the analysis of cancer patient survival we typically estimate excess mortality (observed - expected deaths). The SMR (observed/expected deaths) is a measure of relative mortality. The estimation of observed and expected numbers of deaths are performed in an identical manner for each measure but with the SMR we assume that the effect of exposure is multiplicative to the baseline rate whereas with excess mortality we assume it is additive. Which measure, relative mortality or excess mortality, do you think is more homogeneous across age?

The following example illustrates the approach to estimating SMRs/SIRs using R. Specifically, we will estimate SMRs for the melanoma data using the general population mortality rates stratified by age and calendar period (derived from popmort) to estimate the expected number of deaths. The expected mortality rates depend on current age and current year so the approach is as follows:

1. Split follow-up into 1-year age bands
2. Split the resulting data into 1-year calendar period bands
3. For each age-period band, merge with popmort.dta to obtain the expected mortality rates
4. Sum the observed and expected numbers of deaths and calculate the SMR (observed/expected) and a 95% CI

(a)

Start by splitting the follow-up into 1 year age bands. We do not have a date of birth available, so we can make a mid-point approximation for the age at diagnosis.

data(melanoma) 
mel <- filter(melanoma, stage == "Localised") %>% 
        mutate( dead = (status %in% c("Dead: cancer","Dead: other") & surv_mm <= 120)+0, 
                surv_mm = pmin(120, surv_mm)
              ) 
head(mel) 

## Define the age at start and end of follow-up 
mel <- mutate( mel, adx = age+0.5,   # age at diagnosis  (mid-point approximation) 
               astart = adx, 
               astop  = adx+surv_mm/12 
              )
## Split by age 
mel.split <- survSplit(mel, cut = 1:105, event = "dead", 
                       start = "astart", end = "astop")
## Quick check: the first two ids 
subset(mel.split, id<=2, select = c(id, astart, astop, dead)) 

(b)

Now split these new records into one year calendar period bands.

# For each age time band from (a), we calculate the start and stop in calendar time 
# We calculate the time since diagnosis as difference between age at start/stop and 
# age at diagnosis, and add that interval to year at diagnosis
mel.split2 <- mutate(mel.split, 
                    ystart = ydx + astart - adx, 
                    ystop  = ydx + astop - adx
                  )
subset(mel.split2, id<=2, select = c(id, adx, astart,astop,dead, ydx, ystart, ystop))

## Now we can split along the calendar time 
## For each of the new age-calendar time bands, we now have to adjust the age at 
## start and end in the same way as above 
mel.split2 <- survSplit( mel.split2, cut = 1970:2000, event = "dead", 
                         start = "ystart", end = "ystop" ) %>%
              mutate( astart = adx + ystart - ydx, 
                      astop  = adx + ystop - ydx
                    )
## Quick check: this seems ok 
subset(mel.split2, id<=2, select = c(id, ystart, ystop, astart, astop, dead)) 

(c)

Descriptives: Each subject’s follow–up is now divided into small pieces corresponding to the agebands and calendar periods the subject passes through. We can make tables of deaths and person-years by age and calendar period with

## We calculate the total person time at risk for each time band
mel.split2 <- mutate( mel.split2, 
                      age  = floor(astart),  # Age at which person time was observed 
                      year = floor(ystart),  # Calendar year during which person time was observed 
                      pt  =  ystop - ystart  # ... or astop - astart, works the same
                    ) 
subset(mel.split2, id<=2, select = c(id, ystart, ystop, astart, astop, dead, age, year, pt))

## Now tabulate: sum of person time across all combinations of age & year 
## (for some years, ages) 
xtabs(pt ~ age + year, data=mel.split2, subset = age>=50 & age<60 & year>=1980 & year<1990)

## Same: sum up 0/1 alive/dead for total count of deaths 
xtabs(dead ~ age + year, data=mel.split2, subset = age>=50 & age<60 & year>=1980 & year<1990) 

As the data have been split in one-year intervals on both time scales the table created above is not so informative. Grouped variables will provide a better overview.

(d)

Descriptives: as the data have been split in 1-year intervals on both time scales in the table created above, this takes up a lot of space. Grouped variables will provide a better overview. To make a table of rates by age and calendar period, we can define group variables and use survRate as usually:

mel.split2 <- mutate(mel.split2, 
                     age10  = cut(age, seq(0, 110 ,by=10), right=FALSE), 
                     year10 = cut(year, seq(1970, 2000, by=5), right=FALSE) 
                    ) 
sr <- survRate(Surv(pt, dead) ~ sex + age10 + year10, data=mel.split2) 
rownames(sr) <-1:nrow(sr) ## Simple rownames for display 
head(sr, n = 20) 

(e)

To calculate the expected cases for a cohort, using reference mortality rates classified by age and calendar period, it is first necessary to merge the population rates with the observed person-time. Then the expected number of cases are calculated by multiplying the follow-up time for each record by the reference rate for that record. The SMR is the ratio of the total observed cases to the total number expected.

First, calculate the total person time at risk and the observed number of deaths for each combination:

pt <- mutate(mel.split2, sex = unclass(sex)) %>%    # make sex integer to be in line with popmort 
      group_by(sex, age, year)               %>%    # aggregate by sex, age, year 
      summarise(pt = sum(pt), observed = sum(dead)) # sum the person time, deaths 

## `summarise()` has grouped output by 'sex', 'age'. You can override using the `.groups`
## argument.

pt <- ungroup(pt)  # For convenience
head(pt) 

This is the data set of reference rates:

head(rstpm2::popmort)
summary(rstpm2::popmort) 

We now merge the observed person time and deaths (in the melanoma cohort) with the corresponding reference rates in the popmort data set:

joint <- left_join(pt, rstpm2::popmort)

## Joining, by = c("sex", "age", "year")

head(joint) 

The expected number of events (assuming the reference mortality rates in the general population) are simply:

joint <- mutate(joint, expected = pt * rate) 
head(joint) 

(f)

Calculate the crude overall SMR, as well as the SMR by sex. Plot the SMRs by calendar year and age.

We can calculate the SMR simply from the data, like so:

SMR_all <- summarise(joint, SMR = sum(observed) / sum(expected) ) 
SMR_all 

By sex:

SMR_bySex <- group_by(joint, sex) %>% summarise( SMR = sum(observed) / sum(expected) ) 
SMR_bySex 

We can do the same by e.g. calendar year, but there it is easier to plot than show:

SMR_byYear <- group_by(joint, year) %>% summarise( SMR = sum(observed) / sum(expected) ) 
SMR_byYear 
plot( SMR_byYear, type = "o")  # quick & dirty plot 

For the ages, we have two many levels and too few events ot produce reliable estimates. We can either use splines or (simpler) define a grouping variable as above.

joint <- mutate(joint, age_group = cut(age, seq(0, 110, by=10), right = FALSE))
SMR_byAge <- group_by(joint, age_group) %>% summarise( SMR = sum(observed) / sum(expected) ) 
SMR_byAge 
plot( SMR_byAge)           # quick & dirty plot 
abline( h = 1:2, lty = 2)  # two reference lines at 1 & 2

Note that we can use survRate to get confidence intervals:

by(joint, joint$sex, function(data) poisson.test(sum(data$observed), sum(data$expected)))

(g)

Model the SMR as a function of age, sex and calendar period.

We use the observed (nominator in SMR) as the outcome count and the log-expected count (denominator in SMR) as offset in the Poisson regression.

Let’s prettify the data for the fit. Here, a good choice of reference level for age & year makes a big difference in how easy we can interpret the model parameters;

joint2 <- transform(joint,
                    sex  = factor(sex, levels = 1:2, labels = c("m", "f")),
                    year =  factor(year) %>% relevel(ref = "1985"),  # mid-study
                    age_group = relevel(age_group, ref = "[70,80)")  # Close to one already
                   )
## Model & parameters
summary(fit <- glm(observed ~ sex + year + age_group + offset(log(expected)), data=joint2, family=poisson)) 
eform(fit)

A simple analysis of deviance, which performs a likelihood ratio test for the hypothesis that removing a variable from the model does not improve model fit shows however that all predictors are statistically significantly associated with the SMR:

drop1(fit, test = "Chisq")