arriaga sensitivity, an experiment

sen_arriaga_instantaneous.R

library(tidyverse)

mx_to_lx <- function(mx){
  mx[is.na(mx)] <- 0
  lx <- exp(-cumsum(mx))
  lx <- c(1,lx)
  lx[1:length(mx)]
}
lx_to_dx <- function(lx){
  -diff(c(lx,0))
}
mx_to_Lx <- function(mx){
  lx <- exp(-cumsum(mx))
  lx <- c(1,lx)
  (lx[-1] + lx[-length(lx)]) / 2
}

rcumsum <- function(x){
  rev(cumsum(rev(x)))
}

sen_arriaga <- function(mx1,mx2){
  tibble(mx1,mx2) |> 
  mutate(lx1 = mx_to_lx(mx1),
         lx2 = mx_to_lx(mx2),
         Lx1 = mx_to_Lx(mx1),
         Lx2 = mx_to_Lx(mx2),
         Tx1 = rcumsum(Lx1),
         Tx2 = rcumsum(Lx2),
         # dx1 = lx_to_dx(lx1),
         # dx2 = lx_to_dx(lx2),
         direct = lx1 * (Lx2 / lx2 - Lx1 / lx1),
         # direct2 = lx1 * ((lx2 - dx2 * .5)/ lx2 - (lx1 - dx1 * .5)/ lx1),
         indirect = lead(Tx2) * (lx1 / lx2 - lead(lx1) / lead(lx2)),
         indirect = ifelse(is.na(indirect),0,indirect),
         cc = direct + indirect,
         delta = mx2 - mx1,
         # watch out for 0s in delta denominator
         sen = cc / delta) |> 
    pull(sen)
}

sen_arriaga_instantaneous <- function(mx, perturb = 1e-6){
  scale = 1 - perturb
  (sen_arriaga(mx*scale,mx*1/scale) + 
      sen_arriaga(mx*1/scale,mx*scale) ) / 2
}

mx_to_ex <- function(mx){
  lx <- mx_to_lx(mx)
  Lx <- mx_to_Lx(mx)
  Tx <- rcumsum(Lx) 
  ex <- Tx / lx
  ex
}

mx_to_e0 <- function(mx){
  mx_to_ex(mx) %>% '['(1)
}

# Alternative sensitivity version, more in line with Caswell
https://link.springer.com/chapter/10.1007/978-3-030-10534-1_4 , eq 4.2
sen_e0_mx <- function(mx){
  lx <- mx_to_lx(mx)
  dx <- -diff(c(lx,0))
  Lx <- mx_to_Lx(mx)
  ex <- mx_to_ex(mx)
  ex2 <- mean2(ex)
  -Lx * ex2 #+ direct
}

# compare with numDeriv::grad()
mx <- c(0.004121, 0.000331, 0.000202, 0.000139, 0.000169, 0.000128, 
0.000124, 0.000134, 0.000122, 8.6e-05, 0.000159, 5.9e-05, 0.000125, 
0.000168, 0.000145, 0.000232, 0.000203, 0.00027, 0.000241, 0.000282, 
0.000289, 0.000256, 0.000275, 0.000226, 0.000315, 0.000288, 0.000292, 
0.000321, 0.000308, 0.000294, 0.000466, 0.000378, 0.000445, 0.000509, 
0.000506, 0.000546, 0.000707, 0.000705, 0.000708, 0.00078, 0.001039, 
0.000918, 0.001078, 0.001076, 0.001233, 0.001244, 0.001423, 0.001536, 
0.001538, 0.001736, 0.00181, 0.00209, 0.002103, 0.002271, 0.002342, 
0.002805, 0.002821, 0.002897, 0.003246, 0.003815, 0.003974, 0.004457, 
0.004915, 0.005465, 0.005565, 0.006589, 0.007431, 0.008017, 0.009318, 
0.010441, 0.011762, 0.013013, 0.014838, 0.01661, 0.01947, 0.022126, 
0.025082, 0.028847, 0.03375, 0.038257, 0.044347, 0.050758, 0.05875, 
0.068362, 0.076458, 0.086995, 0.098928, 0.111966, 0.129323, 0.144716, 
0.165768, 0.184866, 0.204536, 0.228662, 0.251223, 0.271861, 0.300944, 
0.310618, 0.319114, 0.34236, 0.378892, 0.413136, 0.414993, 0.506128, 
0.532637, 0.485294, 0.5, 0.832117, 0.5, 1.230769, 0.774194)

# tiny residuals, and we don't know which is the less precise
plot(sen_arriaga_instantaneous(mx,1e-6) - numDeriv::grad(mx_to_e0, mx))

# residuals that have a shape qualitatively like a death distribution:
plot(sen_arriaga_instantaneous(mx,1e-6) - sen_e0_mx(mx))

# same
plot(numDeriv::grad(mx_to_e0, mx) - sen_e0_mx(mx))

# when using LTRE, both sen_arriaga_instantaneous() and numDeriv::grad() would 
# give similar precision, whereas sen_e0_mx() would imply some error.

# The difference is not a clean map to a death distribution: 
# 1. the mode is at a higher age, 
# 2. the shape isn't quite the same, 
# 3. the scale is unknown

# observe:
plot(sen_arriaga_instantaneous(mx,1e-6) - sen_e0_mx(mx))
lines((mx_to_lx(mx) |> lx_to_dx())/6) # 

# I also compared with the "second life" age-at-death death distribution,
# also not the same. It's also not an arithmetic or logarythmic mean of the
# first and second death distributions.

# objective: we'd like a version of sen_e0_mx(), i.e. that directly gives the sensitivity
# for a discrete time lifetable and without this systematic error. What further calculation 
# or transformation is necessary to achieve this?