# Gompertz-Makeham Model Life Expectancy using Numerical Integration and Analytical Solution of Castellares et al. (2020)

I am trying to calculate life expectancies for the Gompertz-Makeham model, but can't replicate the results of the paper which gives the formulas. I use R to calculate the formula using numerical integration. I provide the formulas and my code below and would be happy if someone could check if I am missing something.

Castellares et al. (2020) provide an analytical solution for the life expectancy of the Gompertz-Makeham model:

$$\mathcal{M}_{GM} = \frac{1}{b} \exp\left( \frac{a e^{bx}}{b} \right ) \left( \frac{a e^{bx}}{b} \right)^{c/b} \Gamma\left( - \frac{c}{b}, \frac{a e^{bx}}{b} \right) \qquad \text{with }$$ $$\Gamma(\cdot, \cdot) = \int_x^\infty t^{u-1} e^{-t} dt, \; x>0, \quad u \in \mathbb{R}$$

$$\Gamma(\cdot, \cdot)$$ is called the upper incomplete gamma distribution. Alternatively, the same result can be achieved via numerical integration of: $$\int_0^\infty \exp\left[-ct - \frac{a e^{bx}}{b} (e^{xb} - 1) \right]$$

They estimate the parameters $$a, b$$ and $$c$$ and denote the age by $$x$$:

With my implementation of these formulas in R i can't replicate these findings. The analytical solution gives:

# Analytical Solution for age 30
a = 0.00014
b = 0.11521
c = 0.00033
t = 30

gamma_u = -(c/b)              # First parameter for the upper incomplete gamma distribution
gamma_x = (a * exp(b*t)) / b  # Second parameter for the upper incomplete gamma distribution

# Function to integrate over to get the upper incomplete gamma distribution
intfun = function(t, u){
t^(u-1) * exp(-t)
}

# Result for the upper incomplete gamma distribution
gammainc = integrate(f = intfun,
lower = gamma_x,
upper = Inf,
u = gamma_u)$value # Analytic formula: (1/b) * exp(gamma_x) * (gamma_x^(c/b)) * gammainc # Here the result of the computation: [1] 24.4011  and the numerical integration gives the same: ## Numerical integration intfun2 = function(t, a, b, c, x){ exp( -c*t - (a* exp(b*x))/(b) * (exp(b*t) - 1) ) } integrate(intfun2, lower = 0, upper = Inf, a = 0.00014, b = 0.11521, c = 0.00033, x = 30) # Here the result of the computation: 24.4011 with absolute error < 5.7e-06  ## 1 Answer After some research I came upon a GitHub repository of one of the authors in which he used Bell Regression to model mortality. Albeit the parameters $$a, b$$ and $$c$$ are derived by different methods as in Castellares et al. (2020), the script also contains code that implements life expectancies using the above formula. # Makeham # Gamma-Imcompleta superior Gamma_Inc_sup = function(s, z){ integrate(function(t){ t^(s-1)*exp(-t) }, z, Inf)$value
}

ex_mak = function(t, a, b, c){
(exp(a*exp(b*t)/b)/b)*((a*exp(b*t)/b)^(c/b))*Gamma_Inc_sup(-c/b,a*exp(b*t)/b)
}


This code is equivalent to my code in the above question.

What I found is that the parameter estimation uses mortality data for the ages 30-60, which are the ages in which it is log-linear. Consequently, the code calculates the life expectancy at age 30 using something like ex_mak(0, a, b, c), i.e. t = 0corresponds to age 30.

Using then the values 0, 30 and 60 for the ages 30, 60 and 90 and the parameters as in my above question I find the results:

ex_mak(0, 0.00014, 0.1152, 0.00033)
[1] 52.84945

ex_mak(30, 0.00014, 0.1152, 0.00033)
[1] 24.40485

ex_mak(60, 0.00014, 0.1152, 0.00033)
[1] 4.506404


These are still slightly different to the values in the paper: but I strongly suppose can attribute these to rounding errors.