Survival Analysis – Calculating Incidence Density Rate After Weighting in R

epidemiologypropensity-scoresrsurvival

How to calculate an incidence density/ rate for weighted data, e.g. using weighting methods such as matching weights, IPW, …?

— UPDATE 06.11.2023:
Regard the following paper as an example of incidence density/ rate calculation:
In Table 2 the crude and weighted incidence density/ rate is given as:

Events / Follow-up, person-years per 1000 person-years.

Thus, the crude incidence density/ rate for DPP-4 inhibitors is:
58 / 265 553 * 1000 which results in: 0.22 (rounded to the second).

Lets move on to example data with follow-up time in days:

library(cobalt)
library(dplyr)

set.seed(123)
lalonde <- cbind(lalonde,
                 event = sample(c(0,1), size=614, replace=TRUE, prob=c(0.84,0.16)),
                 time = runif(614, min=10, max=365))

formula <- treat ~ age + educ + race + married + nodegree + re74 + re75 + re78

# PS
lalonde$pscore <- glm(formula, data = lalonde,
                  family = binomial(link = "logit"))$fitted.values

# Calculate weights
lalonde$weight <- ifelse(lalonde$treat == 1,
                         pmin(lalonde$pscore, 1 - lalonde$pscore) / lalonde$pscore,
                     pmin(lalonde$pscore, 1 - lalonde$pscore) / (1 - lalonde$pscore))

Here, the calculation of a crude incidence density/ rate for the treated per 1000 person-years is given:

crude <- lalonde %>%
  filter(treat == 1) %>%
  summarize(event = sum(event == 1),
            persontime = sum(time) / 365.25)

crude %>% 
  summarize(Incidence_density = round(event / persontime * 1000, 1))

resulting in:

  Incidence_density
1             294.5

When using weighted data to calculate the incidence density/ rate for the treated we have to keep in mind the weighted events and weighted person-time:

weighted.event:
The number of events is the sum of the weights for the individuals who had events and were treated.

weighted.persontime:
Person time is the sum of the weights times the person time for those who were treated.

So its just: weighted.event / weighted.persontime

Thus:

weighted <- lalonde %>%
  filter(treat == 1) %>%
  summarize(weighted.event = sum(weight * event),
            weighted.persontime = sum(weight * time) / 365.25)

weighted %>% 
  summarize(Incidence_density = round(weighted.event / weighted.persontime * 1000, 1))

which gives:

  Incidence_density
1             269.6

Is there a reference for the methodology? I do not find any explanation of that.
Further, calculating the weighted incidence density/ rate may need robust standard errors to calculate the CI.
In addition, I am not aware of any package that can perform these calculations.

Thank you

Best Answer

The survey package can do these weighted calculations (conditional outcome probability weights work the same way as sampling weights)

> d<-svydesign(id=~1, weights=~weight, data=lalonde)
> d<-update(d, time_kyr=time/365.25/1000)
> svyratio(~event,~I(time_kyr),design=d)
Ratio estimator: svyratio.survey.design2(~event, ~I(time_kyr), design = d)
Ratios=
      I(time_kyr)
event    270.0646
SEs=
      I(time_kyr)
event    40.04273

Or, if you prefer, estimate the totals with standard errors and do the division

> totals<-svytotal(~event+time_kyr,design=d)
> totals
            total     SE
event    29.88659 4.3266
time_kyr  0.11066 0.0059
> svycontrast(totals, quote(event/time_kyr))
          nlcon     SE
contrast 270.06 40.043

Ratio estimation is classical in survey sampling; the formulas will be in any survey sampling textbook.

Related Question