Confidence interval for Population Attributable Fraction with several strata

Question

I have used aggregated data to create a table of person-years (pys) and deaths by social class, age and sex.

If we consider social class to be a modifiable factor, we can calculate the number of 'expected' deaths in a situation where the low class group has the same mortality rate as the high class group. The difference is the number of attributable deaths. In the example below this is 38, and the Population Attributable Fraction for social class is 38 / 182 = 21%.

+-------+-------+--------+------+--------+--------+----------+--------------+
| Class |  Age  |  Sex   | Pys  | Deaths |  Rate  | Expected | Attributable |
+-------+-------+--------+------+--------+--------+----------+--------------+
| High  | Young | Male   |  100 |     10 |    0.1 |       10 |            0 |
| High  | Young | Female |  120 |     12 |    0.1 |       12 |            0 |
| High  | Old   | Male   |   40 |      8 |    0.2 |        8 |            0 |
| High  | Old   | Female |   80 |     12 |   0.15 |       12 |            0 |
+-------+-------+--------+------+--------+--------+----------+--------------+
| Low   | Young | Male   |  200 |     30 |   0.15 |       20 |           10 |
| Low   | Young | Female |  200 |     30 |   0.15 |       20 |           10 |
| Low   | Old   | Male   |  160 |     40 |   0.25 |       32 |            8 |
| Low   | Old   | Female |  200 |     40 |    0.2 |       30 |           10 |
+-------+-------+--------+------+--------+--------+----------+--------------+
| ALL   | ALL   | BOTH   | 1100 |    182 | 0.1655 |      144 |           38 |
+-------+-------+--------+------+--------+--------+----------+--------------+

Do you know how I would calculate a confidence interval for this fraction? It seems straightforward to calculate a PAF and a confidence interval for a 2x2 matrix (e.g. if we were just looking at rows 1 and 5 of the table), but I'm not sure how to approach the problem with a number of different strata.

Answer 1

I've done a bit of reading and couldn't find a formula. I did find an article looking at estimation of PAF in a prospective study - https://academic.oup.com/aje/article/171/7/837/85797 - which I didn't completely understand but I think you need individual-level data for the regression model-based estimation of relative risks.

I therefore used a Monte Carlo method to estimate confidence intervals, and created a general function:

#--------------------------------------------
# Monte Carlo estimation of 95% CIs for a PAF
#--------------------------------------------

# dat is a data frame
# mod_var is the modifiable variable (string)
# ref is the reference category of the modifiable variable (string)
# adj_vars is a vector a non-modifiable variables (e.g. age, sex)
# time is the time variable (e.g. person-years)
# event is the events variable (e.g. deaths)
# N is the number of Monte Carlo simulations (integer)
# point.estimate (MEAN or MEDIAN) is the method for reporting the point estimate

mc_paf <- function(dat, mod_var, ref, adj_vars, time = 'pys', event = 'deaths', N = 10000, point.estimate = 'MEAN') {

  # deaths and mortality rates
  rd <- sapply(dat[,event], function(x) rpois(N, x))
  rd_rate <- t(t(rd) / dat[,time])

  # which are the relevant reference groups?
  refn <- dat[,c(mod_var, adj_vars)]
  refn$nc1 <- seq_len(nrow(refn))
  refn_ref <- refn[refn[,mod_var] == ref,]
  refn_ref[,mod_var] <- NULL
  names(refn_ref) <- c(adj_vars, 'nc2')
  refn <- merge(refn, refn_ref, by = adj_vars)
  refn <- refn$nc2[order(refn$nc1)]

  # expected
  expected <- t(t(rd_rate[,refn]) * dat[,time])

  # attributable
  attributable <- rd - expected
  attributable <- attributable[,dat[,mod_var] != ref] # remove those in reference category
  attributable <- rowSums(attributable)

  # PAF and results
  pafs <- attributable / rowSums(rd)
  point <- if (point.estimate == 'MEAN') mean(pafs) else median(pafs)
  c(point, quantile(pafs, c(0.025, 0.975)))

}

#------------------------
# Example (from question)
#------------------------

dat <- data.frame(
  class = c(rep('high', 4), rep('low', 4)),
  age = rep(rep(c('young', 'old'), each = 2), 2),
  sex = rep(c('male', 'female'), 4),
  pys = c(100, 120, 40, 80, 200, 200, 160, 200),
  deaths = c(10, 12, 8, 12, 30, 30, 40, 40),
  stringsAsFactors = F
)

# Monte-Carlo confidence intervals with 1M simulations

mc_paf(dat, N = 1000000, mod_var = 'class', ref = 'high', adj_vars = c('age', 'sex'), time = 'pys', event = 'deaths')
# PAF = 0.21 (95% CI 0-0.42)

# Manual calculation of PAF

dat$rate <- dat$deaths / dat$pys
dat$expected <- dat$pys * rep(dat$rate[1:4], 2)
dat$att <- dat$deaths - dat$expected
sum(dat$att) / sum(dat$deaths)

#----------------------------------
# Second example (with more levels)
#----------------------------------

dat2 <- structure(list(smoke = c("current", "current", "current", "ex", 
                         "ex", "ex", "never", "never", "never", "current", "current", 
                         "current", "ex", "ex", "ex", "never", "never", "never"), age = c("50-59", 
                                                                                          "60-69", "70-79", "50-59", "60-69", "70-79", "50-59", "60-69", 
                                                                                          "70-79", "50-59", "60-69", "70-79", "50-59", "60-69", "70-79", 
                                                                                          "50-59", "60-69", "70-79"), sex = c("male", "male", "male", "male", 
                                                                                                                              "male", "male", "male", "male", "male", "female", "female", "female", 
                                                                                                                              "female", "female", "female", "female", "female", "female"), 
               pys = c(1578L, 1448L, 1345L, 2529L, 2340L, 2156L, 945L, 921L, 
                       900L, 1785L, 1645L, 1602L, 1932L, 1777L, 1711L, 1445L, 1390L, 
                       1388L), deaths = c(269L, 281L, 366L, 242L, 273L, 331L, 46L, 
                                          77L, 84L, 256L, 289L, 363L, 147L, 178L, 237L, 99L, 80L, 134L
                       )), class = "data.frame", row.names = c(NA, -18L))

mc_paf(dat2, N = 1000000, mod_var = 'smoke', ref = 'never', adj_var = c('age', 'sex'))
# PAF = 0.43 (95% CI 0.38-0.48)

# Manual calculation of PAF

dat2$rate <- dat2$deaths / dat2$pys
dat2_ref <- dat2[dat2$smoke == 'never',]
dat2_ref[,c('pys', 'deaths', 'smoke')] <- NULL
names(dat2_ref)[3] <- 'ref_rate'
dat2 <- merge(dat2, dat2_ref)
dat2$expected <- dat2$pys * dat2$ref_rate
dat2$att <- dat2$deaths - dat2$expected
sum(dat2$att) / sum(dat2$deaths)