age | lx | dx | mx | qx |
---|---|---|---|---|
0 | 100000 | 553.22 | 0.00556 | 0.00553 |
1 | 99447 | 38.18 | 0.00038 | 0.00038 |
2 | 99409 | 23.16 | 0.00023 | 0.00023 |
age | lx | dx | mx | qx |
---|---|---|---|---|
0 | 100000 | 553.22 | 0.00556 | 0.00553 |
1 | 99447 | 38.18 | 0.00038 | 0.00038 |
2 | 99409 | 23.16 | 0.00023 | 0.00023 |
lx
is the total number living in a hypothetical cohort (of 100,000).dx
is the number of deaths in the age interval.age | lx | dx | mx | qx |
---|---|---|---|---|
0 | 100000 | 553.22 | 0.00556 | 0.00553 |
1 | 99447 | 38.18 | 0.00038 | 0.00038 |
2 | 99409 | 23.16 | 0.00023 | 0.00023 |
mx
is the death rate within the interval.qx
is the probability of death within an interval.age | lx | dx | mx | qx |
---|---|---|---|---|
0 | 100000 | 553.22 | 0.00556 | 0.00553 |
1 | 99447 | 38.18 | 0.00038 | 0.00038 |
2 | 99409 | 23.16 | 0.00023 | 0.00023 |
age | lx | dx | mx | qx |
---|---|---|---|---|
0 | 100000 | 553.22 | 0.00556 | 0.00553 |
1 | 99447 | 38.18 | 0.00038 | 0.00038 |
2 | 99409 | 23.16 | 0.00023 | 0.00023 |
The R packages MortalityLaws
and demography
facilitate life table data aquisition and modeling.
Common sources:
Highlighted column: Remaining life expectancy at exact age x.
age | lx | dx | ex | qx | mx |
---|---|---|---|---|---|
0 | 100000 | 553.221 | 78.963 | 0.00553 | 0.00556 |
1 | 99447 | 38.180 | 78.402 | 0.00038 | 0.00038 |
2 | 99409 | 23.160 | 77.432 | 0.00023 | 0.00023 |
3 | 99385 | 17.689 | 76.450 | 0.00018 | 0.00018 |
Highlighted column: Remaining life expectancy at exact age x.
age | lx | dx | ex | qx | mx |
---|---|---|---|---|---|
0 | 100000 | 553.221 | 78.963 | 0.00553 | 0.00556 |
1 | 99447 | 38.180 | 78.402 | 0.00038 | 0.00038 |
2 | 99409 | 23.160 | 77.432 | 0.00023 | 0.00023 |
3 | 99385 | 17.689 | 76.450 | 0.00018 | 0.00018 |
(We’ll try to replicate these numbers with a discrete time Markov model in a bit.)
qx
) values.params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)
params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)
params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)
params =
list(
t_names = c("lifetable"), # Strategy names.
n_treatments = 1, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 0, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt # Processed HMD life table data (2019, USA)
)
lifetable
“strategy” will draw on the HMD life table data processed earlier.We’ll next construct a transition probability matrix using the age-specific probability of death field (qx
) from the life table data.
To do so, we need to define a function that, for a given age and cycle length, calculates the probability of death and places this within a transition matrix.
Life Table:
sim_cohort <- function(params) {
params$t_names %>% map(~({
tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
res <- do.call(rbind,lapply(params$mP, function(tp) {
tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
}))
res <- rbind(c(params$n_cohort,0),res)
dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
res
})) %>%
set_names(params$t_names)
}
age | qx | ex |
---|---|---|
0 | 0.00553 | 78.963 |
age | qx | ex |
---|---|---|
0 | 0.00553 | 78.963 |
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp
# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = cycle_adj(life_exp_cycle)
Age | Life Expectancy (Life Table data) | Life Expectancy (Markov-Life Table) |
---|---|---|
0 | 78.963 | 78.924 |
MortalityLaws
package has a number of mortality models we can draw from:Name | Model | Code |
---|---|---|
Infant mortality | ||
Opperman | \(\mu[x] = A/\sqrt{x} - B + C\sqrt{x}\) | opperman |
Weibull | \(\mu[x] = 1/\sigma (x/M)^{(M/\sigma - 1)}\) | weibull |
Accident hump | ||
Inverse-Gompertz | \(\mu[x] = [1-e^{(M-x)/\sigma}]/[e^{(M-x)/\sigma}-1]\) | invgompertz |
Inverse-Weibull | \(\mu[x] = 1/\sigma (x/M)^{[-M/\sigma - 1]} / [e^{(x/M)^{(-M/\sigma)}} - 1]\) | invweibull |
Name | Model | Code |
---|---|---|
Adult mortality | ||
Gompertz | \(\mu[x] = A e^{Bx}\) | gompertz |
Gompertz | \(\mu[x] = 1/\sigma e^{(x-M)/\sigma}\) | gompertz0 |
Makeham | \(\mu[x] = A e^{Bx} + C\) | makeham |
Makeham | \(\mu[x] = 1/\sigma e^{(x-M)/\sigma} + C\) | makeham0 |
Perks | \(\mu[x] = [A + BC^x] / [BC^{-x} + 1 + DC^x]\) | perks |
Strehler-Mildvan | \(\mu[x] = K e^{-V_0 (1-Bx)/D}\) | strehler_mildvan |
Adult and/or old-age mortality | ||
Van der Maen | \(\mu[x] = A+Bx+Cx^2 + I/[N-x]\) | vandermaen |
Beard | \(\mu[x] = e^{Bx}/[1+KA e^{Bx}]\) | beard |
Beard-Makeham | \(\mu[x] = A e^{Bx}/[1 + KA e^{Bx}]\) | beard_makeham |
Gamma-Gompertz | \(\mu[x] = A e^{Bx}/(1+AG/B[e^{Bx}-1])\) | ggompertz |
Name | Model | Code |
---|---|---|
Old-age mortality | ||
Van der Maen | \(\mu[x] = A + Bx + I/[N-x]\) | vandermaen2 |
Quadratic | \(\mu[x] = A + Bx + Cx^2\) | quadratic |
Kannisto | \(\mu[x] = A e^{Bx}/[1+A e^{Bx}]\) | kannisto |
Kannisto-Makeham | \(\mu[x] = A e^{Bx}/[1+A e^{Bx}]+C\) | kannisto_makeham |
Name | Model | Code |
---|---|---|
Full age range | ||
Thiele | \(\mu[a] = A e^{-Bx} + C e^{-\frac{1}{2} D(x-E)^2}+F e^{Gx}\) | thiele |
Wittstein | \(q[x] = (1/B) A^{-[(Bx)^N]} + A^{-[(M-x)^N]}\) | wittstein |
Heligman-Pollard | \(q[x]/p[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + G H^x\) | HP |
Heligman-Pollard | \(q[x] = A^{[(x + B)^C]} + D e^{-E \log(x/F)^2} + GH^x / [1 + GH^x]\) | HP2 |
Rogers-Planck | \(q[x] = A_0 + A_1 e^{-Ax} + A_2 e^{B(x-u)-e^{-C(x-u)}}+A_3 e^{Dx}\) | rogersplanck |
Martinelle | \(\mu[x] = [A e^{Bx} + C]/[1+D e^{Bx}]+K e^{Bx}\) | martinelle |
Carriere | \(l[x] = P1 l[x](weibull) + P2 l[x](invweibull) + P3 l[x](gompertz)\) | carriere1 |
Kostaki | \(q[x]/p[x] = A^{[(x+B)^C]} + D \exp[-(E_i \log(x/F_))^2] + G H^x\) | kostaki |
Key takaways:
MortalityLaws
has a number of mortality models it can fit to life table data, depending on the context and population you are modeling.
You can easily experiment around with different models to see how well they fit your mortality data.
Another nice feature of the package is that each mortality model type has its own defined function to calculate the mortality rate (or probability) given an age and model coefficients.
So rather than assume a constant mortality rate within age bins, you can model the mortality rate directly for a 40.3 year old, or a 90.0001 year old.
Generally speaking, we need three inputs:
age
: Ages for lifetable (e.g., 0, 1, 5, …, 100; or 0,1,2,3,…,100).dx
: The number of deaths between exact ages x and x+1.lx
: Number of survivors to exact age x.\[ mx = A \exp(Bx) \]
A B
0.00010772 0.07534542
# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )
A
0.0021938
# Can simply use the supplied function to calcualte
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)
$hx
[1] 0.0021938
$par
A B
0.00010772 0.07534542
$Sx
[1] 0.97269
ages <- lt$age[lt$age>=40 & lt$age<100]
deaths <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]
gom_fit40 <- MortalityLaw(
x = ages, # vector with ages
Dx = deaths, # vector with death counts
Ex = exposure, # vector containing exposures
law = "gompertz",
opt.method = "LF2")
While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP
in lieu of gompertz
!
While you could draw from any of the aforementioned mortality models, the Heligman-Pollard model tends to fit the entire age distribution reasonably well.
Easy to use: just use HP
in lieu of gompertz
!
Mortality rate for a 40 year old:
$hx
[1] 0.0020409
$par
A B C D E F_
0.000386408 0.021271772 0.107422357 0.001025285 2.871760263 33.201213579
G H
0.000022277 1.102611045
[1] 70.467
age | mx | Sx |
---|---|---|
70 | 0.01827 | 0.76590 |
71 | 0.01995 | 0.75077 |
72 | 0.02162 | 0.73472 |
y
for a given value of x
.approxfun()
:
x
is quantile in CDF for survival, i.e., 1 - Survival1y
is the age.qHP <-
approxfun(y = 0:115,
x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)
[1] 71.069
Survival in life table:
age | mx | Sx |
---|---|---|
70 | 0.01827 | 0.76590 |
71 | 0.01995 | 0.75077 |
72 | 0.02162 | 0.73472 |
params =
list(
t_names = c("lifetable", # Strategy names
"heligman-pollard",
"gompertz"),
n_treatments = 3, # Number of treatments
s_names = c("Alive", "Dead"), # State names
n_states = 2, # Number of states
n_cohort = 100000, # Cohort size
initial_age = 40, # Cohort starting age
n_cycles = 110, # Number of cycles in model.
cycle = 1, # Cycle length
life_table = lt, # Processed HMD life table data (2019, USA)
gom_fit = gom_fit, # Fitted Gompertz model object
hp_fit = hp_fit # Fitted Heligman-Pollard model object
)
fn_mPt <- function(t, params) {
with(params, {
h = cycle # Time step
lapply(t, function(tt){
# Life table method
current_age_lt = pmin(initial_age + (tt)*h - 1,max(lt$age))
p_death <- lt[lt$age==current_age_lt,"qx"]
# parametric mortality methods
current_age = initial_age + (tt)*h - 1
r_death_gompertz <- gompertz(current_age,params$gom_fit$coefficients)$hx
p_death_gompertz <- 1 - exp(-r_death_gompertz)
r_death_hp <- HP(current_age, params$hp_fit$coefficients)$hx
p_death_hp <- 1 - exp(-r_death_hp)
array(data = c(1 - p_death, 0,
p_death,1,
1 - p_death_hp, 0,
p_death_hp,1,
1 - p_death_gompertz, 0,
p_death_gompertz,1),
dim = c(n_states, n_states, n_treatments),
dimnames = list(from = s_names,
to = s_names,
t_names))
})
})
}
# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)
# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr*(1 / params$n_cohort)) %*%
payoff_life_exp)
# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
49, 43, 59, 17) / 48
cycle_adj <- function(x) sum(alt_simp_coef(length(x)) * x)
total_life_exp = lapply(life_exp_cycle,cycle_adj)
Constructing the cause-deleted life table requires two necessary inputs:
We extract the percentage of overall deaths that are from CVD by age bin:
ihme_cvd <-
tibble::tribble(
~age_name, ~val,
1, 0.038771524,
5, 0.038546046,
10, 0.044403585,
15, 0.033781126,
20, 0.035856165,
25, 0.053077797,
30, 0.086001439,
35, 0.130326551,
40, 0.184310334,
45, 0.21839762,
50, 0.243705394,
55, 0.256334637,
60, 0.26828001,
65, 0.272698709,
70, 0.28529754,
75, 0.310642009,
0, 0.016750489,
80, 0.353518012,
85, 0.399856716,
90, 0.447817792,
95, 0.495305502
) %>%
mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE)) %>%
select(age_ihme, pct_cvd = val)
Merge these percentages into the underlying life table data and use them to calculate the total number of CVD deaths by age:
age | age_ihme | pct_cvd | dx | dx_i |
---|---|---|---|---|
0 | [0,1) | 0.01675 | 553.221 | 9 |
10 | [10,15) | 0.04440 | 11.518 | 1 |
25 | [25,30) | 0.05308 | 103.706 | 6 |
50 | [50,55) | 0.24371 | 365.721 | 89 |
75 | [75,80) | 0.31064 | 1997.766 | 621 |
98 | [95,105) | 0.49531 | 1262.593 | 625 |
We next calculate the cause-deleted probability of death (qd
) and death rate (md
), as well as the cause-specific probability of death (qi
) and death rate (mi
):
qi
): \(q \cdot dx_i / dx\)qi
): \(q \cdot dx_i / dx\)age | dx | dxd | dx_i | q | qi | md | mi |
---|---|---|---|---|---|---|---|
0 | 553.221 | 544.221 | 9 | 0.00553 | 0.00009 | 0.00546 | 0.00009 |
1 | 38.180 | 37.180 | 1 | 0.00038 | 0.00001 | 0.00037 | 0.00001 |
2 | 23.160 | 22.160 | 1 | 0.00023 | 0.00001 | 0.00022 | 0.00001 |
3 | 17.689 | 16.689 | 1 | 0.00018 | 0.00001 | 0.00017 | 0.00001 |
4 | 14.209 | 13.209 | 1 | 0.00014 | 0.00001 | 0.00013 | 0.00001 |
5 | 13.213 | 12.213 | 1 | 0.00013 | 0.00001 | 0.00012 | 0.00001 |
ages_ <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_ <- lt_$dxd[lt_$age<100 & lt_$age>=25]
exposure_ <- lt_$lx[lt_$age<100 & lt_$age>=25]
mort_fit_CVDdeleted <- MortalityLaw(
x = ages_,
Dx = deaths_, # vector with death counts
Ex = exposure_, # vector containing exposures
law = "HP2",
opt.method = "LF2")
Healthy CVD CVDDeath Dead
0 100000 0.0 0.0000 0.00
1 90391 9506.2 0.2945 103.00
2 81702 18088.0 1.1404 208.58
3 73847 25833.6 2.4853 316.88
4 66745 32822.6 4.2817 428.11
5 60324 39126.6 6.8544 542.45
6 54519 44809.4 11.5751 660.14
7 49270 49930.7 17.3944 781.42
8 44526 54544.1 23.8195 906.57
9 40236 58697.0 31.3741 1035.87
TK
TK