life_tables_mc.R

# this code calculates a life table using mortality rates by single-year-of-age
# confidence intervals are calculated assuming that each year-of-age is an independent sample

#-----------------
# make sample data
#-----------------

set.seed(56)
dat <- data.frame(age = 0:99)
dat$person_years <- 50000 / (1 + exp(-(seq(10, -5, length.out = nrow(dat)))))
dat$person_years <- round(dat$person_years, 0)
dat$theoretical_mortality_rate <- (1.09^(0:99) / (1.09^99 * 4))
dat$deaths_target <- round(dat$person_years * dat$theoretical_mortality_rate, 0)
dat$deaths <- rpois(nrow(dat), dat$deaths_target)
dat <- dat[, c('person_years', 'deaths')]

#----------------
# life expectancy
#----------------

# life table function

life.table <- function(mx, cohort = 100000, EX = T) { # if EX is true, just return life expectancy
  n <- length(mx) + 1
  qx <- 2 * mx / (2 + mx)
  qx <- c(qx, 1) # forced method - mortality rate max age + 1 is 100%
  lx <- c(1, cumprod(1 - qx)) * cohort
  dx <- -c(diff(lx), lx[n] * qx[n])
  t <- (lx + c(lx[-1], 0)) / 2
  Tx <- rev(cumsum(rev(t)))
  ex <- Tx / lx
  if (EX) {
  	return(ex[1])
  } else {
  	return(data.frame(lx = lx, dx = dx, t = t, Tx = Tx, ex = ex)[1:n,])
  	}
}

dat$mortality_rate <- dat$deaths / dat$person_years

lt <- life.table(dat$mortality_rate, EX = F)
# plot(lt$lx) # survival curve
lt$ex[1] # life expectancy = 80.4

#---------------------------------
# monte-carlo confidence intervals
#---------------------------------

# if you have mortality rates rather than person-years and counts of deaths, you'll have to assume a distribution around the mortality rates

B <- 10000 # number of simulations
l <- nrow(dat)
mmc <- rpois(B * l, dat$deaths)
mmc <- matrix(mmc, ncol = B)
mmc <- mmc / dat$person_years
mmc <- apply(mmc, 2, life.table)
quantile(mmc, probs = c(0.025, 0.5, 0.975))

# life expectancy = 80.4 (95% CI 80.2-80.6)