rm(list=ls())
library(magrittr)
library(dplyr)
library(tidyverse)
library(popEpi)
library(biostat3)
library(Epi)Standardized Incidence Ratio (SIR) or Standardized Mortality Ratio (SMR) are common measure of event occurrence in epidemiology. They are used to compare the occurrence of an event in a cohort to that observed in a given population called reference population.
SIR or SMR help the investigators to have a global idea on the occurrence of the event of interest in the followed population (the cohort). They are an indirect method of adjustment that describe in numerical term how the cohort average experience of the event during the follow-up compared with that of the reference population as a whole. E.g. a \(SIR=2\) means that the occurrence of the event event in the cohort is 2 times higher than what is expected in the reference population. But caution, this does not mean any causality in regard of a give exposure in the cohort or a given risk factor for which the cohort is followed for.
Mathematically SIR (SMR) is the ratio of observed events and expected events:
\[SIR = \frac{\sum d_j }{\sum n_j \lambda_j} = \frac{D}{E}\]
where \(D\) is the observed events in the cohort population and \(E\) is the expected number of events.
Observed events are the absolute number of events that occurred in the cohort during the follow up.
The expected number of events are derived by multiplying the cohort person-years with the reference population rate. The population rate should be stratified or adjusted by confounding factors (age or age group, gender, calendar period etc.). The reference population rate in strata \(i\) (\(\lambda_i\)) is defined as: \[\lambda_i = \frac{d_i}{n_i}\] where \(d_i\) is the total observed events and \(n_i\) is the total observed person years (or the size of the population) in the \(i\)th strata of the reference population.
Univariate confidence intervals are based on exact values of Poisson distribution and the formula for p-value is \[ \chi^2 = \frac{ (|O - E| -0.5)^2 }{E} \] Modeled SIR is a Poisson regression model with log-link and cohorts person-years as an offset (we will come to that in the next sections).
Let’s have a quick look to some examples of cohort and reference data:
pop$Rate## [1] 0.0009791478 0.0020579002 0.0041527623 0.0068271197 0.0102507645
## [6] 0.0115581283 0.0172613950 0.0219064535pop$Rate*100000## [1] 97.91478 205.79002 415.27623 682.71197 1025.07645 1155.81283 1726.13950
## [8] 2190.64535| age_group | Number_events | Pers_years |
|---|---|---|
| 35-39 | 0 | 480 |
| 40-44 | 1 | 587 |
| 45-49 | 3 | 680 |
| 50-54 | 5 | 541 |
| 55-59 | 8 | 479 |
| 60-64 | 8 | 356 |
| 65-69 | 4 | 157 |
| 70-74 | 1 | 36 |
In this first example Table 1:
| age_group | Number_events | Pers_years | Rate |
|---|---|---|---|
| 35-39 | 81 | 82725 | 0.0009791 |
| 40-44 | 177 | 86010 | 0.0020579 |
| 45-49 | 336 | 80910 | 0.0041528 |
| 50-54 | 434 | 63570 | 0.0068271 |
| 55-59 | 419 | 40875 | 0.0102508 |
| 60-64 | 599 | 51825 | 0.0115581 |
| 65-69 | 907 | 52545 | 0.0172614 |
| 70-74 | 925 | 42225 | 0.0219065 |
In this second example Table 2 :
Note: make sure the same variables with with the same labels and the same classes in both your cohort data and in your population or reference data.From the Table 2, we can extract the Rate column to merge with the cohort data Table 1 and we can calculate the Expected number of events obtained as population reference Rate in the \(i\)th strata multiplied by the cohort \(j\)th strata Pers_years. This is expressed as :
\(Expected\) = \(n_i\)*\(\lambda_j\)
|
|
The final step is now to calculate the SIR value as sum of Number_events divided by the sum of number of Expected :
\[\begin{equation} \begin{split} SIR &=\frac{\sum d_j }{\sum n_j \lambda_j} \\ & = \frac{D}{E} \\ & =\frac{(0+1+3+5+8+8+4+1)}{(0.4699909+1.207987+2.823878+3.693472+4.910116+4.114694+2.710039+0.7886323)} \\ & = 1.44796 \end{split} \end{equation}\]
Interpretation :
In our example, \(SIR = 1.45\) means that standardizing on the age group, incidence in our studied cohort is 1.45 times higher than in the reference population.
In this post, I will presente you three different way of SIR analyses on R.
First of all, let me remember you that raw data for a given study are never in the form we saw in the Table 1 and Table 2. Most the time, you ought to configure your data in way to obtain Table 1 and Table 2 as in our examples. I named this step of data manipulation or data configuration as a data management step.
Overall I will presente you the following steps:
From this step up to the end, I will be using a simulated cohort data named mySIR and a simulated population reference data my.pop.ref.
Using these data allows us only to compute SMR analysis with the research question:
is the cancer incidence rate among females diagnosed with rectal cancer the same as that in the general population?
The cohort data simulated mySIR is on cancer occurrence in patients undergoing medical diagnostic ionizing radiation exposure from between 2000-2015.
These are variables from the mySIR database :
sex: the gender of the patient (1 = Male, 2 = female)
birth_date: the date of birth (date: dd-mm-yyyy)
entry_date: the date of exposure to ionizing radiation (date: dd-mm-yyyy)
exit_date: the date of exit from follow-up (death or censoring) (date: dd-mm-yyyy)
status: the status of the person at exit; 0 exit (or censored) without cancer; 1 diagnosed with cancer
glimpse(mySIR) ## Rows: 14,060
## Columns: 5
## $ sex <dbl> 2, 1, 1, 1, 2, 2, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1,...
## $ birth_date <date> 2003-03-10, 2007-08-01, 2009-07-30, 1989-05-14, 2006-04...
## $ entry_date <date> 2003-07-17, 2011-05-09, 2010-06-07, 2004-07-22, 2007-10...
## $ exit_date <date> 2014-12-31, 2014-12-31, 2014-12-31, 2007-05-14, 2014-12...
## $ status <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...head(mySIR)## sex birth_date entry_date exit_date status
## 1 2 2003-03-10 2003-07-17 2014-12-31 0
## 2 1 2007-08-01 2011-05-09 2014-12-31 0
## 3 1 2009-07-30 2010-06-07 2014-12-31 0
## 4 1 1989-05-14 2004-07-22 2007-05-14 0
## 5 2 2006-04-10 2007-10-04 2014-12-31 0
## 6 2 2006-08-08 2009-10-19 2014-12-31 0Overview of mySIR dataset :
From these variables in the dataset, let’s create the following variables:
age: the age at exit from the cohort (ex_date - bi_date).
year: the calendar year at exit from the cohort (ex_date - dg_date).
Follow_up : the duration of follow-up in years of each patient of the cohort (exit_date - entry_date).
age_exit : age at exit from the follow-up in years (exit_date - birth_date).
mySIR %<>%
mutate(year=format(exit_date, format="%Y"),# to extract only the year (YYYY) of exit from the follow-up
age_exit= cal.yr(exit_date) - cal.yr(birth_date), # using the function cal.yr from the package "Epi" to compute the age at exit from the follow-up in years.
age = cut(age_exit, breaks = c(0, 1, 5, 10, 15, 20), # To create age groups, coded 1 for <1 year, 2 for 1 to 4 years, 3 for 5 to 9 years, 4 for 10 to 14 years and 5 for >= 15 years old.
labels = c(1, 2, 3, 4, 5), right = FALSE),
Follow_up = (exit_date - entry_date)/365.25) # using the function cal.yr from the package "Epi" to compute the Follow-up in years.The reference population data simulated my.pop.ref is on cancer incidence in the general population from 2000-2015, by age group and by gender. The followings could be found from the dataset :
sex: the gender of the patient (1 = Male, 2 = female)
year: calendar year (format “YYYY”)
age: patients age group at cancer diagnosis, coded 1 for <1 year, 2 for 1 to 4 years, 3 for 5 to 9 years, 4 for 10 to 14 years and 5 for >= 15 years old.
Number_cases: the total number of cases (\(n_i\)) in each age, year and sex strata of the reference population data
Population: the size of the population (\(d_i\)) in each age, year and sex strata of the reference population data
Rate: the average population cancer incidence rate per person-year (\(\lambda_i\)=\(n_i\)/\(d_i\)), where \(n_i\) is the number of deaths and \(d_i\) is the person-years). \(n_i\) and \(d_i\) are not explicitly in the data set, we only the resultant rate \(\lambda_i\).
glimpse(my.pop.ref)## Rows: 128
## Columns: 6
## Groups: sex, year [32]
## $ sex <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ year <dbl> 2000, 2000, 2000, 2000, 2001, 2001, 2001, 2001, 2002, ...
## $ age <dbl> 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, 3, 4, 1, 2, ...
## $ Population <dbl> 380740.0, 1481551.5, 1862478.0, 1979979.0, 386279.0, 1...
## $ Number_cases <dbl> 109, 285, 247, 252, 107, 282, 248, 279, 78, 318, 243, ...
## $ Rate <dbl> 0.0002862846, 0.0001923659, 0.0001326190, 0.0001272741...head(my.pop.ref)## # A tibble: 6 x 6
## # Groups: sex, year [2]
## sex year age Population Number_cases Rate
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2000 1 380740. 109 0.000286
## 2 1 2000 2 1481551. 285 0.000192
## 3 1 2000 3 1862478. 247 0.000133
## 4 1 2000 4 1979979. 252 0.000127
## 5 1 2001 1 386279. 107 0.000277
## 6 1 2001 2 1500713. 282 0.000188Overview of my.pop.ref dataset :
We can also obtain the summary by using the function lexpand of the package popEpi which allows quickly summary of cohort data.
cohort_lex <- dplyr::select(mySIR, -c("age", "age_exit"))
lex.sir <- lexpand(cohort_lex,
birth = birth_date,
entry = entry_date,
exit = exit_date,
status = status==1,
breaks = list(fot=0:20, # cut the follow-up time in periods of one year
age=c(0, 1, 5, 10, 15, 20), # cut the age at exit in subgroups, coded 1 for <1 year, 2 for 1 to 4 years, 3 for 5 to 9 years, 4 for 10 to 14 years and 5 for >= 15 years old.
per=2000:2020), # To cut the calendar year in periods of one year. You can choose also to cut in subgroups of 5 or 10 years
aggre = list(sex,
age,
year=per)# rename the per variable in year as in the population data
)
glimpse(lex.sir)## Rows: 150
## Columns: 7
## $ sex <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ age <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1...
## $ year <int> 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009...
## $ pyrs <dbl> 66.43120, 118.50406, 144.38630, 139.66301, 150.89047, 172....
## $ at.risk <dbl> 0, 112, 138, 135, 139, 149, 191, 184, 165, 139, 152, 171, ...
## $ from0to0 <dbl> 10, 12, 9, 9, 20, 12, 15, 17, 11, 19, 15, 20, 12, 13, 0, 1...
## $ from0to1 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0...The object lex.sir is an aggregated cohort data. The original cohort data mySIR is splited by the function lexpand into to subintervals of time over calendar time (year), age, and follow-up time (fot) with given time breaks using splitMulti. lex.sir
In the lex.sir object, the from0to1 column represents the total number of observed events recorded and the pyrs column represents the total number of person-time recorded in each sex, year and age group strata.
Since we computed our cohort data and population reference data similar to tables 1 and 2, we can now easily calculate the SIR in our cohort according to the first equation.
Remember that both the population reference data my.pop.ref and the cohort data lex.sir should contain the same column names with the same classes and the same codifications for at least the variables of adjustment (sex, age, year).
To calculate the SIRs, we will now merge by “sex”, “year” and “age”, the population reference data my.pop.ref and the cohort data lex.sir and then create a new column “Expected” as the Rate that multiply person-time.
lex.sir.final <- merge(lex.sir, my.pop.ref[ , c("sex", "year", "age", "Rate")],
by=c("sex", "year", "age"), all.x=TRUE)
lex.sir.final$Expected <- lex.sir.final$Rate*lex.sir.final$pyrsNow since we have the “Observed” and “Expected” columns ready in our dataset, we can calculate the overall SIR as the sum of “Observed” divided by the sum of “Expected”.
sum(lex.sir.final$from0to1, na.rm = TRUE)/sum(lex.sir.final$Expected, na.rm = TRUE) # Exclude NA values## [1] 3.376423We can also calculate SIRs in different ways and for different variables as follow :
lex.sir.final %>%
group_by(sex) %>%
summarise(O=sum(from0to1),
E=sum(Expected, na.rm = TRUE),
sir_sex=sum(from0to1)/sum(Expected, na.rm = TRUE))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 2 x 4
## sex O E sir_sex
## <dbl> <dbl> <dbl> <dbl>
## 1 1 26 7.11 3.66
## 2 2 17 5.63 3.02lex.sir.final %>%
group_by(age) %>%
summarise(O=sum(from0to1),
E=sum(Expected, na.rm = TRUE),
sir_age=sum(from0to1)/sum(Expected, na.rm = TRUE))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 5 x 4
## age O E sir_age
## <dbl> <dbl> <dbl> <dbl>
## 1 1 1 0.955 1.05
## 2 2 16 5.16 3.10
## 3 3 9 3.67 2.45
## 4 4 9 2.95 3.05
## 5 5 8 0 Inflex.sir.final %>%
group_by(year) %>%
summarise(O=sum(from0to1),
E=sum(Expected, na.rm = TRUE),
sir_year=sum(from0to1)/sum(Expected, na.rm = TRUE))## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 15 x 4
## year O E sir_year
## <int> <dbl> <dbl> <dbl>
## 1 2000 0 0.0743 0
## 2 2001 0 0.188 0
## 3 2002 3 0.310 9.67
## 4 2003 0 0.426 0
## 5 2004 3 0.560 5.36
## 6 2005 3 0.664 4.52
## 7 2006 2 0.783 2.55
## 8 2007 7 0.893 7.84
## 9 2008 4 0.999 4.00
## 10 2009 3 1.06 2.82
## 11 2010 4 1.18 3.38
## 12 2011 4 1.27 3.16
## 13 2012 4 1.35 2.97
## 14 2013 4 1.49 2.69
## 15 2014 2 1.49 1.34Since observed cases in the cohort are a typical case of Poisson distribution, Poisson model can be used to estimate the SIRs as follow:
fit=glm(from0to1~offset(log(Expected)), data=lex.sir.final, family = "poisson")
eform(fit)## exp(beta) 2.5 % 97.5 %
## 2.748252 1.973267 3.827605fit=glm(from0to1[sex==1]~offset(log(Expected[sex==1])), data=lex.sir.final, family = "poisson") # for males
eform(fit)## exp(beta) 2.5 % 97.5 %
## 2.953930 1.926009 4.530459fit=glm(from0to1[sex==2]~offset(log(Expected[sex==2])), data=lex.sir.final, family = "poisson") # females
eform(fit)## exp(beta) 2.5 % 97.5 %
## 2.488359 1.473737 4.201518Using Poisson regression, the SIRs are little bit reduced (underestimated). But I would not be able to explain why do we have such difference.
We can also use Poisson test to compute our SIRs
by(lex.sir.final, lex.sir.final$sex, function(data) poisson.test(sum(data$from0to1), sum(data$Expected, na.rm = TRUE))) ## lex.sir.final$sex: 1
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 26, time base = 7.1092, p-value = 3.846e-08
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 2.389035 5.358714
## sample estimates:
## event rate
## 3.657247
##
## ------------------------------------------------------------
## lex.sir.final$sex: 2
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 17, time base = 5.6262, p-value = 8.273e-05
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 1.760181 4.837841
## sample estimates:
## event rate
## 3.021579by(lex.sir.final, lex.sir.final$age, function(data) poisson.test(sum(data$from0to1), sum(data$Expected, na.rm = TRUE)))## lex.sir.final$age: 1
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 1, time base = 0.95504, p-value = 0.6152
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 0.02650956 5.83391074
## sample estimates:
## event rate
## 1.047072
##
## ------------------------------------------------------------
## lex.sir.final$age: 2
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 16, time base = 5.1585, p-value = 9.826e-05
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 1.772863 5.036891
## sample estimates:
## event rate
## 3.101653
##
## ------------------------------------------------------------
## lex.sir.final$age: 3
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 9, time base = 3.6721, p-value = 0.01311
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 1.120713 4.652596
## sample estimates:
## event rate
## 2.450913
##
## ------------------------------------------------------------
## lex.sir.final$age: 4
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 9, time base = 2.9497, p-value = 0.003412
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## 1.395190 5.792076
## sample estimates:
## event rate
## 3.051172
##
## ------------------------------------------------------------
## lex.sir.final$age: 5
##
## Exact Poisson test
##
## data: sum(data$from0to1) time base: sum(data$Expected, na.rm = TRUE)
## number of events = 8, time base = 0, p-value < 2.2e-16
## alternative hypothesis: true event rate is not equal to 1
## 95 percent confidence interval:
## Inf Inf
## sample estimates:
## event rate
## Infby(lex.sir.final, lex.sir.final$year, function(data) poisson.test(sum(data$from0to1), sum(data$Expected, na.rm = TRUE)))SIR using the function sir of the popEPi package:
sr <- popEpi::sir(lex.sir.final, coh.obs = 'from0to1',
coh.pyrs = 'pyrs',
ref.data = my.pop.ref[ , c("sex", "year", "age", "Rate")],
ref.rate = Rate,
print = c("sex"),
adjust = c("age", "sex", "year"),
test.type = "homogeneity",
conf.type = "wald",
conf.level = 0.95, EAR = F)
sr## SIR (adjusted by age, sex, year) with 95% confidence intervals (wald)
## Test for homogeneity: p = 0.538
##
## Total sir: 3.38 (2.44-4.55)
## Total observed: 43
## Total expected: 12.74
## Total person-years: 90608
##
##
## sex observed expected pyrs sir sir.lo sir.hi p_value
## 1: 1 26 7.11 46856.07 3.66 2.49 5.37 0
## 2: 2 17 5.63 43751.69 3.02 1.88 4.86 0sr <- popEpi::sir(lex.sir.final, coh.obs = 'from0to1',
coh.pyrs = 'pyrs',
ref.data = my.pop.ref[ , c("sex", "year", "age", "Rate")],
ref.rate = Rate,
print = c("age"),
adjust = c("age", "sex", "year"),
test.type = "homogeneity",
conf.type = "wald",
conf.level = 0.95, EAR = F)## Model fitting failed. Univariate confidence intervals selected. (zero values in expected)sr## SIR (adjusted by age, sex, year) with 95% confidence intervals (univariate)
## Could not test homogeneity
##
## Total sir: 3.38 (2.44-4.55)
## Total observed: 43
## Total expected: 12.74
## Total person-years: 90608
##
##
## age observed expected pyrs sir sir.lo sir.hi p_value
## 1: 1 1 0.96 3879.30 1.05 0.03 5.83 0.6415
## 2: 2 16 5.16 25739.74 3.10 1.77 5.04 0.0000
## 3: 3 9 3.67 30299.17 2.45 1.12 4.65 0.0118
## 4: 4 9 2.95 21901.16 3.05 1.40 5.79 0.0012
## 5: 5 8 0.00 8788.39 Inf Inf Inf 0.0000