The incidence_rates() function calculates standardised
incidence rates, including standardising by age-at-risk.
The function uses Epi::Lexis() and
Epi::splitLexis() to expand each participant row into
separate rows for each age-at-risk group (defined by the
ages argument). A new column called agegrp is
created which can be used to calculate separate rates for each
age-at-risk group or to standardise the rates by age group.
By default, the function uses the whole data set to create the
weights. Or the weights argument can be used to provide
weights from a reference/standard population.
The function returns a list of two data frames:
incidence_rates_in_all_strata containing a row for
each stratum formed by the grouping variable(s) (specified by the
group argument) and the standardisation variable(s)
(specified by the adjust_for argument), with results
columns:
standardised_incidence_rates containing a row for each group formed by the grouping variable(s), with results columns:
The incidence incidence rates in this data frame are multiplied by a
scale factor set by the per_pyar argument (default:
1000).
The function gives a warning when any stratum contains fewer than 20 events, and when any stratum contains zero events. Results may be unreliable in these cases. See Confidence intervals for available methods for calculating confidence intervals which may be more suitable when there are few events.
To calculate incidence rates identify the columns that contain
participants’ date of birth (argument: col.dob, default:
“dob_anon”), date of study entry (col.entry_date,
“study_date”), date of study exit (col.exit_date,
“endpoint_date”) and status at exit (col.status,
“endpoint”).
Not specifying the group and adjust_for
arguments will calculate overall incidence rate. The standardised
incidence rate will equal the crude incidence rate.
# $incidence_rates_in_all_strata
# # A tibble: 1 × 6
# all std_pyar pyar nevents incidencerate weight
# <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 1 264327. 264327. 2271 0.00859 1
#
# $standardised_incidence_rates
# # A tibble: 1 × 7
# crude_nevents crude_pyar crude_ir std_ir std_ir_se std_ir_lci std_ir_uci
# <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 2271 264327. 8.59 8.59 0.180 8.24 8.94Not specifying the adjust_for argument will calculate
the crude incidence rate for each group. The standardised incidence rate
for each group will equal the crude incidence rate for that group. For
example, by BMI group:
incidence_rates(ckbtools_participant_data,
col.exit_date = "ihd_date",
col.status = "ihd",
group = "bmi_grp")# $incidence_rates_in_all_strata
# # A tibble: 4 × 7
# all std_pyar bmi_grp pyar nevents incidencerate weight
# <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 1 264327. <18.5 22826. 136 0.00596 1
# 2 1 264327. 18.5-25 168758. 1378 0.00817 1
# 3 1 264327. 25-30 68371. 683 0.00999 1
# 4 1 264327. 30+ 4372. 74 0.0169 1
#
# $standardised_incidence_rates
# # A tibble: 4 × 8
# bmi_grp crude_nevents crude_pyar crude_ir std_ir std_ir_se std_ir_lci
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 <18.5 136 22826. 5.96 5.96 0.511 4.96
# 2 18.5-25 1378 168758. 8.17 8.17 0.220 7.73
# 3 25-30 683 68371. 9.99 9.99 0.382 9.24
# 4 30+ 74 4372. 16.9 16.9 1.97 13.1
# # ℹ 1 more variable: std_ir_uci <dbl>To calculate incidence rates by age-at-risk groups (30-49, 50-59,
60-69, 70-79 years), set group = "agegrp". Any follow-up
time outside these age groups is ignored.
incidence_rates(ckbtools_participant_data,
col.exit_date = "ihd_date",
col.status = "ihd",
ages = c(30, 50, 60, 70, 80),
group = "agegrp")# $incidence_rates_in_all_strata
# # A tibble: 4 × 7
# all std_pyar agegrp pyar nevents incidencerate weight
# <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 1 260331. [30,50) 78644. 184 0.00234 1
# 2 1 260331. [50,60) 85910. 530 0.00617 1
# 3 1 260331. [60,70) 64178. 688 0.0107 1
# 4 1 260331. [70,80) 31599. 708 0.0224 1
#
# $standardised_incidence_rates
# # A tibble: 4 × 8
# agegrp crude_nevents crude_pyar crude_ir std_ir std_ir_se std_ir_lci
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 [30,50) 184 78644. 2.34 2.34 0.172 2.00
# 2 [50,60) 530 85910. 6.17 6.17 0.268 5.64
# 3 [60,70) 688 64178. 10.7 10.7 0.409 9.92
# 4 [70,80) 708 31599. 22.4 22.4 0.842 20.8
# # ℹ 1 more variable: std_ir_uci <dbl>To calculate incidence rates by BMI group, standardised by sex (column “is_female”) and 10 year age groups (and excluding follow-up time where participants are younger than 40 or older than 79):
incidence_rates(ckbtools_participant_data,
col.exit_date = "ihd_date",
col.status = "ihd",
ages = seq(40, 80, 10),
group = "bmi_grp",
adjust_for = c("is_female", "agegrp"),
method = "norm")# Warning: is_female is numeric and has 2 levels, so treating as a group variable.
# ℹ Change is_female to a factor to suppress this warning.# Warning: There are strata with zero events.# $incidence_rates_in_all_strata
# # A tibble: 32 × 8
# is_female agegrp std_pyar bmi_grp pyar nevents incidencerate weight
# <int> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 0 [40,50) 28706. <18.5 2304. 0 0 0.115
# 2 0 [50,60) 36260. <18.5 3272. 14 0.00428 0.145
# 3 0 [60,70) 26611. <18.5 2759. 9 0.00326 0.107
# 4 0 [70,80) 12913. <18.5 1700. 27 0.0159 0.0517
# 5 1 [40,50) 39156. <18.5 2628. 3 0.00114 0.157
# 6 1 [50,60) 49647. <18.5 3681. 19 0.00516 0.199
# 7 1 [60,70) 37568. <18.5 3367. 20 0.00594 0.151
# 8 1 [70,80) 18686. <18.5 1889. 29 0.0154 0.0749
# 9 0 [40,50) 28706. 18.5-25 18612. 42 0.00226 0.115
# 10 0 [50,60) 36260. 18.5-25 23328. 134 0.00574 0.145
# # ℹ 22 more rows
#
# $standardised_incidence_rates
# # A tibble: 4 × 8
# bmi_grp crude_nevents crude_pyar crude_ir std_ir std_ir_se std_ir_lci
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 <18.5 121 21600. 5.60 5.04 0.467 4.13
# 2 18.5-25 1268 159496. 7.95 7.93 0.223 7.50
# 3 25-30 642 64328. 9.98 10.3 0.410 9.54
# 4 30+ 66 4123. 16.0 16.6 2.05 12.6
# # ℹ 1 more variable: std_ir_uci <dbl>incidence_rates(ckbtools_participant_data,
col.exit_date = "ihd_date",
col.status = "ihd",
ages = seq(40, 80, 10),
group = "bmi_grp",
adjust_for = c("is_female", "agegrp"),
method = "gamma")# Warning: is_female is numeric and has 2 levels, so treating as a group variable.
# ℹ Change is_female to a factor to suppress this warning.# Warning: There are strata with zero events.# $incidence_rates_in_all_strata
# # A tibble: 32 × 8
# is_female agegrp std_pyar bmi_grp pyar nevents incidencerate weight
# <int> <fct> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 0 [40,50) 28706. <18.5 2304. 0 0 0.115
# 2 0 [50,60) 36260. <18.5 3272. 14 0.00428 0.145
# 3 0 [60,70) 26611. <18.5 2759. 9 0.00326 0.107
# 4 0 [70,80) 12913. <18.5 1700. 27 0.0159 0.0517
# 5 1 [40,50) 39156. <18.5 2628. 3 0.00114 0.157
# 6 1 [50,60) 49647. <18.5 3681. 19 0.00516 0.199
# 7 1 [60,70) 37568. <18.5 3367. 20 0.00594 0.151
# 8 1 [70,80) 18686. <18.5 1889. 29 0.0154 0.0749
# 9 0 [40,50) 28706. 18.5-25 18612. 42 0.00226 0.115
# 10 0 [50,60) 36260. 18.5-25 23328. 134 0.00574 0.145
# # ℹ 22 more rows
#
# $standardised_incidence_rates
# # A tibble: 4 × 7
# bmi_grp crude_nevents crude_pyar crude_ir std_ir std_ir_lci std_ir_uci
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 <18.5 121 21600. 5.60 5.04 4.17 6.06
# 2 18.5-25 1268 159496. 7.95 7.93 7.50 8.38
# 3 25-30 642 64328. 9.98 10.3 9.56 11.2
# 4 30+ 66 4123. 16.0 16.6 12.8 21.2Note how the standardised incidence rates are exactly the same, but the confidence interval limits are slightly different.
To calculate incidence rates by BMI group, standardised by 10 year age groups with external weights:
std_pop_w <- data.frame(agegrp = c("[40,50)",
"[50,60)",
"[60,70)",
"[70,80)"),
std_pyar = c(50, 30, 20, 10))
incidence_rates(ckbtools_participant_data,
col.exit_date = "ihd_date",
col.status = "ihd",
ages = seq(40, 80, 10),
group = "bmi_grp",
adjust_for = "agegrp",
weights = std_pop_w,
method = "gamma")# Warning: There are strata with fewer than 20 events.# $incidence_rates_in_all_strata
# # A tibble: 16 × 7
# agegrp std_pyar bmi_grp pyar nevents incidencerate weight
# <chr> <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 [40,50) 50 <18.5 4933. 3 0.000608 0.455
# 2 [50,60) 30 <18.5 6953. 33 0.00475 0.273
# 3 [60,70) 20 <18.5 6127. 29 0.00473 0.182
# 4 [70,80) 10 <18.5 3588. 56 0.0156 0.0909
# 5 [40,50) 50 18.5-25 43114. 99 0.00230 0.455
# 6 [50,60) 30 18.5-25 54915. 314 0.00572 0.273
# 7 [60,70) 20 18.5-25 41260. 441 0.0107 0.182
# 8 [70,80) 10 18.5-25 20206. 414 0.0205 0.0909
# 9 [40,50) 50 25-30 18618. 63 0.00338 0.455
# 10 [50,60) 30 25-30 22617. 169 0.00747 0.273
# 11 [60,70) 20 25-30 15774. 194 0.0123 0.182
# 12 [70,80) 10 25-30 7318. 216 0.0295 0.0909
# 13 [40,50) 50 30+ 1197. 6 0.00501 0.455
# 14 [50,60) 30 30+ 1423. 14 0.00984 0.273
# 15 [60,70) 20 30+ 1017. 24 0.0236 0.182
# 16 [70,80) 10 30+ 486. 22 0.0453 0.0909
#
# $standardised_incidence_rates
# # A tibble: 4 × 7
# bmi_grp crude_nevents crude_pyar crude_ir std_ir std_ir_lci std_ir_uci
# <fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 <18.5 121 21600. 5.60 3.85 3.16 4.73
# 2 18.5-25 1268 159496. 7.95 6.41 6.04 6.80
# 3 25-30 642 64328. 9.98 8.50 7.83 9.22
# 4 30+ 66 4123. 16.0 13.4 10.2 17.4Note that there must be columns in the weights data
frame that match the adjust_for variables, and one row for
each combination of levels of those variables, and a column named
std_pyar with the person-years at risk for that stratum in
the reference/standard population.
Multiple columns can be specified in group to calculate
incidence rates for all combinations. For example, to calculate
incidence rates by sex and 10 year age groups (between 30 and 80
years):
incidence_rates(ckbtools_participant_data,
col.dob = "dob_anon",
col.entry_date = "study_date",
col.exit_date = "ihd_date",
col.status = "ihd",
ages = seq(30, 80, 10),
group = c("is_female", "agegrp"))# Warning: is_female is numeric and has 2 levels, so treating as a group variable.
# ℹ Change is_female to a factor to suppress this warning.# Warning: There are strata with fewer than 20 events.# $incidence_rates_in_all_strata
# # A tibble: 10 × 8
# all std_pyar is_female agegrp pyar nevents incidencerate weight
# <dbl> <dbl> <int> <fct> <dbl> <dbl> <dbl> <dbl>
# 1 1 260331. 0 [30,40) 4596. 10 0.00218 1
# 2 1 260331. 0 [40,50) 28706. 62 0.00216 1
# 3 1 260331. 0 [50,60) 36260. 213 0.00587 1
# 4 1 260331. 0 [60,70) 26611. 278 0.0104 1
# 5 1 260331. 0 [70,80) 12913. 290 0.0225 1
# 6 1 260331. 1 [30,40) 6188. 3 0.000485 1
# 7 1 260331. 1 [40,50) 39156. 109 0.00278 1
# 8 1 260331. 1 [50,60) 49647. 317 0.00639 1
# 9 1 260331. 1 [60,70) 37568. 410 0.0109 1
# 10 1 260331. 1 [70,80) 18686. 418 0.0224 1
#
# $standardised_incidence_rates
# # A tibble: 10 × 9
# is_female agegrp crude_nevents crude_pyar crude_ir std_ir std_ir_se
# <int> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
# 1 0 [30,40) 10 4596. 2.18 2.18 0.688
# 2 0 [40,50) 62 28706. 2.16 2.16 0.274
# 3 0 [50,60) 213 36260. 5.87 5.87 0.402
# 4 0 [60,70) 278 26611. 10.4 10.4 0.627
# 5 0 [70,80) 290 12913. 22.5 22.5 1.32
# 6 1 [30,40) 3 6188. 0.485 0.485 0.280
# 7 1 [40,50) 109 39156. 2.78 2.78 0.267
# 8 1 [50,60) 317 49647. 6.39 6.39 0.359
# 9 1 [60,70) 410 37568. 10.9 10.9 0.539
# 10 1 [70,80) 418 18686. 22.4 22.4 1.09
# # ℹ 2 more variables: std_ir_lci <dbl>, std_ir_uci <dbl>If there have been \(n\) events during \(t\) person-years at risk, then the incidence rate is given by \[ r = \frac{n}{t} \text{ .} \] The standardised incidence rate across strata (indexed by \(i\)) is calculated as \[ \textrm{SIR} = \frac{\sum_i{r_i \times s_i}}{\sum_i{s_i}} \] where \(r_i = n_i / t_i\) and \(s_i\) are the person-years at risk for each stratum in a standard population.
By default, this function will calculate the \(s_i\) using all participants in the
produced data, and calculate standard incidence rates for groups defined
by group.
The function returns the weight for each stratum calculated as \(s_i/ \sum_i{s_i}\).
By treating the number of events in each stratum as a Poisson variable, the variance of the SIR can be calculated by \[ \textrm{var}(\textrm{SIR}) = \sum_i \left( \left(\frac{w_i}{\sum w_i}\right)^2 \times \frac{r_i}{t_i}\right) \text{ ,} \] then the standard error is given by \(\sqrt{\textrm{var}}\). Using this standard error to construct confidence intervals by a normal approximation may not be suitable when there are few events or rates are close to zero, or there are a small number of events within strata.
Other methods for calculating confidence intervals including the
gamma method of Fay and Feuer (1997). This is available in the
wspoissonTest() function in the asht package,
and can be used in incidence_rates() by setting the
argument method = "gamma".
Using a normal approximation method is not suitable when there are fewer than 100 events. Methods similar to the gamma method of Fay and Feuer (1997) have been shown to be suitable when there are 10 or more events (Morris et al. 2018).
Dates are converted to a numerical value, giving the calendar year as a fractional number, assuming that years are all 365.25 days long, so inaccuracies may arise see [Epi::cal.yr] for more details