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$:

Parameter values

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) )


          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 1


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){
  }, z, Inf)$value

ex_mak = function(t, a, b, c){

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: Paper Results but I strongly suppose can attribute these to rounding errors.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.