Calculating crude incidence rate

To illustrate, we will now use the colon dataset from the survival package.

library(survival)
library(dplyr)
library(broom)

Running ?survival::colon tells us the following:

Data from one of the first successful trials of adjuvant chemotherapy for colon cancer

OK, so now that we understand the data, let’s start calculating crude incidence rates for death among the different treatment groups:

# Using the colon dataset from the survival package
# Only keep records related to the death outcome
colon_death <- colon %>% dplyr::filter(etype == 2)
# Time is divided by 365.25/100 to get the time in days variable
# first to years, then to 100 person-years
colon_death %>% group_by(rx) %>%
summarise(Events = sum(status),
Time = sum(time/365.25/100),
Rate = Events / Time,
lower = poisson.test(Events, Time)$conf.int[1],
upper = poisson.test(Events, Time)$conf.int[2])
## # A tibble: 3 x 6
## rx Events Time Rate lower upper
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Obs 168 13.8 12.2 10.4 14.2
## 2 Lev 161 13.7 11.7 10.0 13.7
## 3 Lev+5FU 123 15.0 8.22 6.83 9.80

Now we compare to the calculated rates with rates obtained from the survRate() function from the biostat3 package:

library(biostat3)
survRate(Surv(time/365.25/100, status) ~ rx, data = colon_death) %>%
dplyr::select(rx, event, tstop, rate, lower, upper) %>%
as_tibble() %>%
dplyr::rename(Events = event,
Time = tstop,
Rate = rate)
## # A tibble: 3 x 6
## rx Events Time Rate lower upper
## <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Obs 168 13.8 12.2 10.4 14.2
## 2 Lev 161 13.7 11.7 10.0 13.7
## 3 Lev+5FU 123 15.0 8.22 6.83 9.80

Good, the incidence rates are identical. The observational patients had an mortality incidence rate of 12.2 per 100 person-years, compared to the Lev+5-FU treated patients with an incidence rate of 8.22 per 100 person-years. Now, let’s try and repeat these results with poisson regression.

Obtaining estimated incidence rates using poisson regression

Here we use the broom package tidy function to obtain exponentiated estimates:

poisson_fit <- glm(status ~ rx + offset(log(time/365.25/100)),
family = poisson, data = colon_death)
tidy(poisson_fit, exponentiate = T, conf.int = T)
## # A tibble: 3 x 7
## term estimate std.error statistic p.value conf.low conf.high
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 12.2 0.0772 32.4 3.13e-230 10.4 14.1
## 2 rxLev 0.965 0.110 -0.324 7.46e- 1 0.777 1.20
## 3 rxLev+5FU 0.675 0.119 -3.32 9.16e- 4 0.534 0.850

The Intercept estimate here is the estimated IR for the reference level, i.e. the Obs group.

To get estimated IR of Lev+5FU treated:

lev_5fu <- predict(poisson_fit,
newdata = data.frame(rx = "Lev+5FU", time = 36525),
type = "link", se.fit = T)
as_tibble(lev_5fu) %>% summarise(Treatment = "Lev+5FU",
IR = exp(fit),
lower = exp(fit - (1.96 * se.fit)),
upper = exp(fit + (1.96 * se.fit)))
## # A tibble: 1 x 4
## Treatment IR lower upper
## <chr> <dbl> <dbl> <dbl>
## 1 Lev+5FU 8.22 6.88 9.80

Here, the confidence interval needs to be calculated on the \(\log\) scale and then exponentiated back. This will cause the confidence interval to not be centered around the estimate.

A poisson model is actually modeling \(\log\text{incidence rates (ratios)}\) when we use the time variable as an offset. Therefore, we can include covariates in the model to be accounted for, such as age and sex.

Age- and sex-adjusted incidence rates using poisson regression

First, we’ll do it using my age_sex_adjust() function

source("age-sex-adjust.R")
# Usage: age_sex_adjust(data, group, age, sex, event, time)
age_sex_adjust(colon_death, rx, age, sex, status, time)
## # A tibble: 3 x 4
## rx IR lower upper
## <chr> <dbl> <dbl> <dbl>
## 1 Obs 12.2 10.2 14.7
## 2 Lev 11.8 10.1 13.8
## 3 Lev+5FU 8.25 6.99 9.73

Here we see that the adjusted rates are very similar to the crude rates calculated above. Since this data comes from a randomized trial, this is expected and can be taken as a sign that the randomization worked.

Now, let’s do the some thing but without using the ready made function to see how it works under the hood.

poisson_fit <- glm(status ~ rx * I(age^2) + sex,
offset = log(time/365.25/100),
data = colon_death,
family = poisson)
std_fit <- stdGlm(poisson_fit, data = transform(colon_death, time = 36525), X = "rx")
std_sum <- summary(std_fit, CI.type = "log")
std_sum[["est.table"]] %>%
as_tibble(rownames = "rx") %>%
transmute(Treatment = rx,
IR = Estimate,
lower = `lower 0.95`,
upper = `upper 0.95`)
## # A tibble: 3 x 4
## Treatment IR lower upper
## <chr> <dbl> <dbl> <dbl>
## 1 Obs 12.2 10.2 14.7
## 2 Lev 11.8 10.1 13.8
## 3 Lev+5FU 8.25 6.99 9.73

The numbers are identical to the ones obtained from the age_sex_adjust() function, which is logical since we did the same thing as the function does.
A few finishing notes. Here I included age as a quadratic term, and as an interaction with exposure. These are modeling decisions one will have to take, however the model could have been only a main effects model such as:

\[\log\lambda = \beta_0 + \beta_1\text{rxLev} + \beta_2\text{rxLev+5FU} + \beta_3\text{age} + \beta_4\text{sex}\]
Regarding the interaction term, a good explanation was given in the Stata forum in this post.

For anyone who wants to read more, I recommend the course material from the PhD course Biostatistics III at Karolinska Institutet, available here.