Directly Standardized Rates for Recurrent Outcomes

Overview

The dsrrec() function computes directly standardized rates for recurrent and non-recurrent outcomes.

For recurrent outcomes, confidence intervals are computed using the negative binomial approach to variance estimation for directly standardized rates by Stukel et al. (1994).

For non-recurrent outcomes, the gamma distribution approach outlined in Fay et al. (1997) is used.

Pre-Analysis

The dsrrec() function expects person-level event counts and unit-times. On each person-level record, the following variables should be present:

standardization variables (e.g. age, sex, etc).
subgroup variables, if any.

Note: For the standard or reference population, the data must be aggregated by the standardization variables and the unit-time variable name must labeled pop. See example below.

Example Data

The readmission dataset from frailtypack contains rehospitalization times (in days) following surgery in patients diagnosed with colorectal cancer. We will use this data to examine directly (sex) standardized rates for rehospitalizations by Dukes’ tumoral stage, a subgroup variable.

#load necessary packages
require(frailtypack)
require(dplyr)
data(readmission)

Since the data are in a person-period format, we’ll have to calculate the total number of events and total observation time for each person. This will result in a single record per person that summarizes the relevant information.

#Calculate total events and total observation times per person
treadm <- readmission %>%
          group_by(id) %>%
          filter(max(enum)==enum ) %>%
          mutate(events=enum-1,time=t.stop) %>%
          select(id, events, time, sex, dukes)

#View first 6 records
knitr::kable(head(treadm), caption='Person-level Dataset')

Person-level Dataset
id	events	time	sex	dukes
1	2	1037	Female	D
2	1	1182	Male	C
3	1	783	Male	C
4	4	2048	Female	A-B
5	1	1144	Female	C
6	3	1407	Male	A-B

The next step is to form the standardized or reference population, which for this analysis, is the total observation time across our study. Since we are standardizing by sex, we’ll compute total observation time by sex. As noted previously, the unit-time must be labelled pop.

#Make the reference data
tref <- treadm %>% 
        group_by(sex) %>% 
        mutate(pop=sum(time)) %>% 
        select(sex, pop) %>% 
        distinct(sex, pop)


#View
knitr::kable(tref, caption='Reference Dataset')

Reference Dataset
sex	pop
Female	180198
Male	233093