# Makeham Mortality Models as Mixtures

Silvio C. Patricio<EMAIL_ADDRESS>Interdisciplinary Center on Population
Dynamics, University of Southern Denmark Trifon I. Missov  Interdisciplinary
Center on Population Dynamics, University of Southern Denmark

###### Abstract

Background: The Makeham term serves as a fundamental component in mortality
modeling, offering a constant additive hazard that accounts for background
mortality factors usually unrelated to the aging process. This term, widely
employed in mortality analysis, provides a crucial mechanism for capturing
mortality risks unrelated to age-related deterioration. Objective: The
objective of this paper is to explore the relationship between Makeham
mortality models and competing risk frameworks. It aims to investigate how
Makeham models, which are widely used for studying mortality, can be
understood and analyzed within the context of competing risks. The paper seeks
to provide insights into the mathematical properties, interpretation, and
applicability of Makeham models in modeling mortality risks associated with
various causes of death. Additionally, the paper aims to demonstrate formally
that competing risk models can be represented as mixture models, thereby
facilitating a deeper understanding of mortality dynamics.
Contribution:Makeham mortality models, when represented as mixtures, have a
straightforward specification that can easily be extended to account for
unobserved heterogeneity. They provide a semantically and computationally
convenient platform to disentangle senescent from extrinsic (premature)
mortality. By expressing Makeham models as a convex combination of probability
distributions, we are able to estimate the age profile of premature mortality,
especially at the oldest ones, where we intuitively assume that most deaths
are senescent. We are also able to estimate the senescent mortality component,
which is the one to focus on when studying the aging process.
Keywords: Makeham mortality models; Mixture models; Mortality modeling;
Competing risks; Non-senescent mortality.

##  1 Introduction

Death can occur due to various causes, each characterized by its own failure
time. Individuals are exposed to a finite number of such causes, and the
interplay of the latter determines the observed length of life.

Let there be a finite number of causes of death labeled
$c_{1},c_{2},\dots,c_{n}.$ Suppose a non-negative random variable $T_{k}$ with
a hazard function $h_{k}(x)$ captures the age at death from cause $c_{k}$,
$k=1,2,\dots,n$, when $c_{k}$ is the only cause of death. Failure time $T_{k}$
measures the ‘’absolute potency” of cause $c_{k}$. The length of life for an
individual is captured then by the random variable
$Y=\min(T_{1},\dots,T_{n})$, whose hazard function, $\mu(x)$, reflects the
competing $n$ cause-specific risks of dying (see, for example, Chiang, , 1961;
Berman, , 1963; Gail, , 1975).

Makeham models are widely used for studying mortality as they account for a
simple and straightforward competing-risk framework. They are characterized by
a hazard function of the type

$\mu(x)=h(x)+c\,,$ (1)

where $h(x)$ is another hazard function. In the context of mortality, $h(x)$
most often pertains to the aging process, while $c$, the Makeham term,
accounts for extrinsic (or background) hazards that usually capture premature
deaths. When

$h(x)=ae^{bx}\,,$ (2)

(1) reduces to the Gompertz-Makeham (GM) model (Gompertz, , 1825; Makeham, ,
1860), while for

$h(x)=\frac{ae^{bx}}{1+\frac{a\gamma}{b}(e^{bx}-1)}\,,$ (3)

(1) yields the gamma-Gompertz-Makeham ($\Gamma$GM) model (Vaupel et al., ,
1979). Beard,  (1959) and Kannisto,  (1994) provide alternative logistic
curves with different asymptotics. The hazard $h(x)$ can be extended by an
additive component that reflects infant and childhood mortality, as in Siler’s
model (Siler, , 1979)

$h(x)=a_{1}e^{-b_{1}x}+ae^{bx}\,.$ (4)

The GM, $\Gamma$GM, and Siler models reflect a competing risk framework: an
individual dies either as a result of biological processes at early or late
ages or due to some extrinsic risk $c$, whatever strikes first.

The Makeham term in all models described by the general framework (1) accounts
for a competing exponentially distributed risk. Elandt-Johnson,  (1976); Gail,
(1975); Hakulinen and Rahiala,  (1977) discuss mathematical properties and
practical implications of assuming independence among failure times, as well
as the consequences when this assumption is violated. They do not study,
though, the representation of competing risk models as mixtures, in which
every distribution characterizes the distribution of deaths due to a specific
out of of $n$ potential competing risks.

##  2 The Relationship

Suppose a population is subject to $n$ causes of death
$c_{1},c_{2},\dots,c_{n}$, where each individual is characterized by a
corresponding vector random variable
$\bm{T}=(T_{1},T_{2},\dots,T_{n})^{\prime}$ representing the hypothetical
failure time from each cause of death in the absence of other causes. Then, a
non-negative continuous random variable $Y$, capturing the length of life for
an individual, can be represented as $Y=\min(T_{1},T_{2},\dots,T_{n})$, its
hazard function can be expressed as

$\mu(x)=h(x,1)+h(x,2)+\cdots+h(x,n)\,,$ (5)

and its probability density function (PDF) can be represented as

$f(x)=\sum_{j=1}^{n}h(x,j)\exp\left\\{-\sum_{k=1}^{n}H(x,k)\right\\}\,,$ (6)

where $H(x,k)=\int\limits_{0}^{x}h(t,k)dt$ is the $k$-th cumulative sub-hazard
function. (6) implies that the PDF of $Y$ can be expressed as a mixture

$f(x)=\sum_{j=1}^{n}\pi_{j}g_{j}(x)\,,$ (7)

where $g_{j}(x)=f(x\,|\,j)$ denotes the PDF of the lifespan for individuals
who die from $c_{j}$, $j=1,2,\dots,n$ while being exposed to all $n$ risks,
and $\pi_{j}$ is the probability that cause $c_{j}$ strikes first among all
competing causes.

##  3 Proof

We want to demonstrate that a competing risk model can be represented as a
mixture model. In a competing risk setting, individuals are exposed to
multiple causes of death, and death occurs when the first cause strikes. This
is analogous to a mixture model, where each cause of death corresponds to a
component distribution, and lifetimes are modeled as a mixture of these
component distributions.

We will prove that the density of the lifetime distribution, representing the
time until the first from multiple competing causes strikes, can be expressed
as a convex combination of probability distributions corresponding to dying
from each cause of failure in the presence of the risks from all others. Each
component distribution corresponds to the lifetime of individuals who die from
a specific cause. The weights $\pi_{j}$ in the convex combination capture the
probabilities of each cause $c_{j}$ striking first.

###### Proof.

It is given that a population is subjected to $n$ causes of death, labeled
$c_{1},c_{2},\dots,c_{n}$, and each individual is characterized by a vector of
random variables $\mathbf{T}=(T_{1},\dots,T_{n})^{\prime}$, representing the
failure times $T_{j}$ according to each cause of death $c_{j}$ in the absence
of other competing causes, $j=1,\ldots,n$. We denote the hazard function for
each $T_{j}$ by

$h_{j}(t)=\lim_{\Delta\downarrow 0}\frac{1}{\Delta}\mathbb{P}(t<T_{j}\leq
t+\Delta|T_{j}>t)\,,$

and the corresponding survival function by $S_{j}(t)$. The survival function
of $\mathbf{T}$ is given by

$S_{\mathbf{T}}(t_{1},\cdots,t_{n})=\mathbb{P}\left(\bigcap_{j=1}^{n}[T_{j}>t_{j}]\right),$

which satisfies $S_{\mathbf{T}}(0,\cdots,0)=1$ and
$S_{\mathbf{T}}(\infty,\cdots,\infty)=0$, with each $t_{j}\in(0,\infty)$. The
function $S(\cdot)$ is continuous from the right and monotonically non-
increasing in each argument.

Note that we cannot observe the failure times $T_{1},\dots,T_{n}$
simultaneously. As a result, $S_{\mathbf{T}}(t_{1},\cdots,t_{n})$ cannot be
observed, nor can its form be tested (Elandt-Johnson, , 1976). Instead, as
highlighted by Moeschberger and David,  (1971), what we observe is the
individual lifetime denoted by $Y=\min(T_{1},\dots,T_{n})$. The corresponding
survival function is given by

$S_{Y}(t)=\mathbb{P}\left(Y>t\right)=S_{\mathbf{T}}(t,\cdots,t)\,,$ (8)

where $t$ is the observed failure time.

With the assumption that $\mathbb{P}(T_{i}=T_{j})=0$ for all $i\neq j$, we can
define the random variable $J$ as the index of the smallest $T_{i}$.
Therefore, $J=j$ implies that cause $c_{j}$ is responsible for the death, and
$T_{j}<T_{i}$ for all $i\neq j$, $j=1,\dots,n$. By applying the law of total
probability and Bayes’ theorem to (8), we obtain

$S_{Y}(t)=\sum_{j=1}^{n}\mathbb{P}\left(Y>t,J=j\right)=\sum_{j=1}^{n}\mathbb{P}\left(Y>t\,|\,J=j\right)\mathbb{P}\left(J=j\right).$

The function

$S(t\,|\,j)=\mathbb{P}\left(Y>t\,|\,J=j\right)=\frac{\mathbb{P}(t<T_{j}<\min_{i\neq
j}T_{i})}{\mathbb{P}(T_{j}<\min_{i\neq j}T_{i})},$

represents the cause-specific survival function for cause $c_{j}$, while the
quantity

$\pi_{j}=\mathbb{P}\left(J=j\right)=\mathbb{P}(T_{j}<\min_{i\neq j}T_{i})$

is the probability that cause $c_{j}$ strikes first among all competing
causes. Thus,

$S_{Y}(t)=\sum_{j=1}^{n}\pi_{j}S(t\,|\,j)\,.$ (9)

Assuming $S_{\mathbf{T}}$ is differentiable, we can express the density of $Y$
as

$f_{Y}(t)=-\frac{d}{dt}S_{Y}(t)=\sum_{j=1}^{n}f(t,j)=\sum_{j=1}^{n}\pi_{j}f(t\,|\,j)\,,$
(10)

where $f(t\,|\,j)=-\frac{d}{dt}S(t\,|\,j)$ represents the density of lifetimes
due to competing cause $c_{j}$. Integrating $f(t,j)$ with respect to $t$, we
can derive $\pi_{j}$. The hazard function of $Y$ is given by

$h_{Y}(t)=\frac{f_{Y}(t)}{S_{Y}(t)}=\sum_{j=1}^{n}h(t,j)\,,$ (11)

where $h(t,j)=\frac{f(t,j)}{S_{Y}(t)}$. This function is sometimes called the
$j$th sub-hazard rate and should not be confused with the function
$h(t\,|\,j)=-\frac{d}{dt}\log S(t\,|\,j)=f(t\,|\,j)/S(t\,|\,j)$.

(9)-(11) illustrate that, no matter if the failure times $T_{1},\dots,T_{n}$
are independent or not, the distribution of lifetimes can be expressed as a
convex combination of $n$ probability distributions, each representing the
time to death for individuals affected by a given cause $c_{j}$,
$j=1,\ldots,n$.

∎

### Note 1: The assumption on the identity of the forces of mortality

As previously mentioned, the function $S_{\mathbf{T}}(t_{1},\cdots,t_{n})$ is
not directly observable, nor can its form be tested. In practical
applications, an additional assumption on the identity of the forces of
mortality is typically made:

$h(t,j)=h_{j}(t),\quad j=1,\dots,n.$

This assumption implies that $S_{Y}(t)$ can be expressed as a product of the
individual survival functions $S_{j}(t)$:

$S_{Y}(t)=S_{\mathbf{T}}(t,\cdots,t)=\prod_{j=1}^{n}S_{j}(t)\,.$

This assumption suggests that the events $[T_{j}>t],\,j=1,\dots,n$ are
independent. Note that this does not necessarily imply independence of the
random variables $T_{1},\dots,T_{n}$ (see, for example, Gail, , 1975).

As noted by Elandt-Johnson,  (1976), directly observing the ’time-to-death
from a given cause of death’ random variable is not feasible, and estimating
its associated survival function requires additional assumptions. In addition,
Crowder,  (1991) demonstrated that even when the marginal survival functions
$S_{j}(t_{j})$ of the failure time $T_{j}$ are known, observations of $Y$ and
$J$ do not determine their joint survival function $S_{\bm{T}}$. Therefore,
assuming independence between events like $[T_{j}>t],\,j=1,\dots,n$, or even
between the random variables $T_{i}$ and $T_{j}$, for $i\neq j,\,j=1,\dots,n$,
though not necessarily reflecting the true interaction mechanism among
different causes of death, provides a reasonable framework for distinguishing
between extrinsic and senescent mortality.

### Note 2: On independence assumption for cause-specific failure times

When we assume that the random variables $T_{j},j=1,2,\dots,n$, are mutually
independent, the survival function $S_{\bm{T}}(t_{1},\cdots,t_{n})$ simplifies
to $\prod_{j=1}^{n}S_{j}(t_{j})$, making the proof of the relationship rather
straightforward. However, assuming independence of cause-specific failure
times has at least two implications: first, the associated cause-specific
risks of dying act independently of one another, and second, eliminating one
cause of death does not affect the force of mortality of the other cause.
However, the adequacy of this assumption varies depending on the specific
cause of death and its complex relationship with other causes. A comprehensive
model should acknowledge that (a) a simple constant term, such as the Makeham
term, may not accurately capture the complex mechanisms underlying premature
deaths; (b) lifetimes may not be independent, highlighting the need for a more
sophisticated modeling approach; and (c) the dependence of lifetimes,
reflecting behavioral, social, and environmental risk factors, may vary across
different age groups.

According to Cox,  (1959), when considering two risks, data collected on
variables $Y$ and $J$ will not provide evidence contradicting the assumption
that the failure times associated with these risks are independent of each
other. In other words, observations of $Y$ and $J$ will not suggest dependence
between the failure times linked to these risks.

Even when circumventing the assumption of independence of lifetimes by using
the method proposed by Gumbel and Mustafi,  (1967) for random variables
$T_{1}$ and $T_{2}$ following extreme value distributions, the estimation of
correlation parameters is complex (Tawn, , 1988). Thus, applying such a model
to aggregate (life-table) mortality data presents significant challenges.

##  4 Related results

Makeham models are a special case of (5) when $n=2$, $h_{1}(x)=c$ and
$h_{2}(x)=h(x)$. For the Makeham models, we consider a non-negative continuous
random variable $Y$ characterized by a Makeham hazard function as described in
equation (1). This hazard function represents the combined risk of two causes
of death: extrinsic, denoted by $c_{1}$, and biological (senescent), denoted
by $c_{2}$. Each cause is described by its own failure time, captured by the
not necessarily independent random variables $T_{1}$ and $T_{2}$,
respectively. Subsequently, the probability density function (PDF) of
$Y=\min(T_{1},T_{2})$ can be represented as:

$f(x)=c\exp\left\\{-cx-H(x)\right\\}+h(x)\exp\left\\{-cx-H(x)\right\\},$ (12)

where $H(x)=\int\limits_{0}^{x}h(t)dt$ is the cumulative hazard function of
$T_{2}$. From (12), it is easy to represent a Makeham model of type (1) as a
mixture model:

$f(x)=\pi g_{1}(x)+(1-\pi)g_{2}(x),$ (13)

where $g_{j}(x)=f(x\,|\,j)$ denotes the PDF of the lifespan for individuals
who die from $c_{j}$, $j=1,2$, being exposed to both the senescent and the
premature risk of dying.

Each component of the mixture (13) has a straightforward interpretation. For
example,

$\pi=\int_{0}^{\infty}c\exp\left\\{-cx-H(x)\right\\}dx\,,$ (14)

known as the mixing proportion (McLachlan et al., , 2019), represents the
probability of dying from an extrinsic cause in the presence of a competing
senescent risk of dying, i.e., $\pi$ is the premature mortality prevalence. On
the other hand, the functions

$g_{1}(x)=\frac{c}{\pi}\exp\left\\{-cx-H(x)\right\\}\hskip
8.53581pt\text{and}\hskip
8.53581ptg_{2}(x)=\frac{h(x)}{1-\pi}\exp\left\\{-cx-H(x)\right\\}$ (15)

capture the PDFs of the distribution of deaths according to each general
cause: in the presence of competing premature and senescent risks of dying,
$g_{1}$ represents the PDF of non-aging-related (premature) deaths, while
$g_{2}$ denotes the PDF of aging-related (senescent) deaths.

By representing the class of models (1) as a mixture of two probability
distribution functions, we can determine the threshold age, $x^{*}$, that
separates the age interval with prevailing premature deaths from the one with
predominant senescent mortality:

$x^{*}=\max\\{\max\\{x:h(x)=c\\},0\\}\,.$ (16)

We can also express the proportion of non-aging-related deaths at age $x$ as

$p(x)=c/(h(x)+c)\,.$ (17)

$p(x)$ reflects the age-specific prevalence of premature mortality. Table 1
presents closed-form expressions for $x^{*}$ and $p(x)$ for four popular
Makeham mortality models.

Table 1: Expressions for the threshold age, $x^{*}$, separating the age intervals with prevailing premature and senescent mortality, as well as the proportion of non-aging-related deaths at age $x$, $p(x)$, in Gompertz-Makeham, gamma-Gompertz-Makeham, Beard-Makeham, Kannisto-Makeham, and Siler model settings. Model | $h(x)$ | $x^{*}$ | $p(x)$
---|---|---|---
Gompertz-Makeham | $ae^{bx}$ | $\frac{1}{b}\ln\frac{c}{a}$ | $\frac{c}{ae^{bx}+c}$
$\Gamma$-Gompertz-Makeham | $\frac{ae^{bx}}{1+\gamma\frac{a}{b}\left(e^{bx}-1\right)}$ | $\frac{1}{b}\ln\frac{c(b-a\gamma)}{a(b-c\gamma)}$ | $\frac{c}{\frac{ae^{bx}}{1+\gamma\frac{a}{b}\left(e^{bx}-1\right)}+c}$
Beard-Makeham | $\frac{ae^{bx}}{1+kae^{bx}}$ | $\frac{1}{b}\ln\frac{c}{a(1-kc)}$ | $\frac{c}{\frac{ae^{bx}}{1+kae^{bx}}+c}$
Kannisto-Makeham | $\frac{ae^{bx}}{1+ae^{bx}}$ | $\frac{1}{b}\ln\frac{c}{a(1-c)}$ | $\frac{c}{\frac{ae^{bx}}{1+ae^{bx}}+c}$
Siler | $a_{1}e^{-b_{1}x}+a_{2}e^{b_{2}x}$ | no closed-form | $\frac{c}{a_{1}e^{-b_{1}x}+a_{2}e^{b_{2}x}+c}$

Note that assuming a constant extrinsic risk of dying, does not imply a
constant force of mortality for the distribution of premature deaths (in the
presence of a competing senescent risk). The force of mortality for this
subpopulation is given by

$h(x\,|\,1)=\frac{\exp\left\\{-cx-H(x)\right\\}}{\int_{x}^{\infty}\exp\left\\{-cy-H(x)\right\\}dy}\,.$
(18)

In many Makeham mortality models, we can derive closed-form expressions for
different characteristics of the distributions of deaths in the mixture. For
example, the remaining life expectancy at age $x$ is given by

$e_{x}=\frac{1}{S(x)}\int_{0}^{\infty}S(x+t)\,dt\,,$ (19)

where $S(\cdot)$ is the survival function. Substituting
$S(x)=S(x\,|\,1)=\int_{x}^{\infty}g_{1}(y)\,dy$ and
$S(x)=S(x\,|\,2)=\int_{x}^{\infty}g_{2}(y)\,dy$ in (19) yields $e_{x}$ for the
subpopulation struck by the premature and senescent competing risk,
respectively. The modal age at death for these subpopulations corresponds to
the maximum of $g_{1}(x)$ and $g_{2}(x)$, respectively. Note that for any non-
decreasing $h(x)$, the modal age at death for the non-senescent subpopulation
is zero. For the senescent subpopulation (with PDF $g_{2}$), the modal age at
death is the age where the overall force of mortality equals the relative
derivative of $h(x)$ with respect to age:

$h(x)+c=\frac{dh(x)/dx}{h(x)}\,.$ (20)

For the Gompertz-Makeham model, if $b>c$, the senescent modal age at death is
given by

$M_{GM}=\frac{1}{b}\,\ln\frac{b-c}{a}\,,$ (21)

while for the gamma-Gompertz-Makeham model, it is given by

$M_{\Gamma
GM}=\frac{1}{b}\,\ln\left(\frac{b}{a}\cdot\frac{b-a\gamma-c\left(1-\frac{a\gamma}{b}\right)}{b+a\gamma
c}\right)\,.$ (22)

When $c=0$, the expressions for $M_{GM}$ and $M_{\Gamma GM}$ reduce to the
modal age at death for the Gompertz and the gamma-Gompertz model, respectively
(see Missov et al.,  (2015) for details). The modal age at death can also be
expressed for the Beard-Makeham model as

$M_{BM}=\frac{1}{b}\,\ln\frac{b-c}{a(1+ck)}\,,$ (23)

while in the Siler model setting, there is no closed-form expression. Note
that, given the improvements in human mortality and the rise of longevity, we
may observe a decreasing trend for $c$ over time. This leads to a convergence
of the overall modal age at death to the senescent one.

We can also derive closed-form expressions for $\pi$ in a Gompertz-Makeham and
a gamma-Gompertz-Makeham setting. From (14), we can express $\pi=ce_{0}$,
where $e_{0}$ is the life expectancy at birth of the Makeham model. Using the
life expectancy formulae in Castellares et al.,  (2020), we can derive
expressions for $\pi$ in terms of special functions. For the Gompertz-Makeham
model, we have that

$\pi=\frac{c}{a}\left(\frac{a}{b}\right)^{\frac{c}{b}}e^{\frac{a}{b}}\Gamma\left(-\frac{c}{b},\frac{a}{b}\right)$
(24)

where $\Gamma(u,x)=\int_{x}^{\infty}t^{u-1}e^{-t}dx$, $x>0$, $u\in\mathbb{R}$
is a complementary incomplete gamma function. For the gamma-Gompertz-Makeham
model, the closed-form expression for $\pi$ is given by

$\pi=\frac{c\gamma}{b+c\gamma}\,{}_{2}F_{1}\left(\frac{1}{\gamma},1;\frac{1}{\gamma}+\frac{c}{b}+1;1-\frac{a\gamma}{b}\right)$
(25)

where
${}_{2}F_{1}(m,p;q;z)=\frac{\Gamma(q)}{\Gamma(p)\Gamma(q-p)}\int^{1}_{0}u^{p-1}(1-u)^{q-p-1}(1-zu)^{-m}du$
is the Gaussian hypergeometric function.

###  4.1 The Gompertz-Makeham model

The Makeham term has first been used to adjust Gompertz mortality estimates
(Makeham, , 1860), resulting in the Gompertz-Makeham model. A non-negative
continuous random variable $Y$ has a Gompertz-Makeham distribution if its
hazard function is given by

$\mu(x)=ae^{bx}+c,\quad a,b>0,\hskip 5.69046ptc\geq 0;\hskip 11.38092ptx\geq
0,$ (26)

and survival function

$S(x)=\exp\left\\{-\frac{a}{b}\left(e^{bx}-1\right)-cx\right\\}\,.$ (27)

From (26) and (27), it is straightforward to derive the Gompertz-Makeham PDF:

$f(x)=c\exp\left\\{-\frac{a}{b}\left(e^{bx}-1\right)-cx\right\\}+ae^{bx}\exp\left\\{-\frac{a}{b}\left(e^{bx}-1\right)-cx\right\\}.$
(28)

According to Castellares et al.,  (2022), there are four different shapes for
$f(x)$: (a) it may have a local maximum that is not located at the boundary;
(b) it may have a local maximum both at the boundary and inside the parameter
space; (c) it may have a global maximum at the boundary; and (d) it may have a
global maximum at the boundary and an inflection point inside the parameter
space (Norström, , 1997). The four shapes are displayed in Figure 1. While
only panel (a) in Figure 1 corresponds to the pattern observed in human
mortality, panels (b)-(d) were included to illustrate the diverse dynamics of
mortality when decomposing senescence and premature mortality.

Figure 1: PDF of the Gompertz-Makeham distribution as a mixture. The grey line
presents $f$; the red line indicates $g_{1}$, and the blue line represents
$g_{2}$. For (a) we have $a=0.0005$, $b=0.1$ and $c=0.005$, for (b) we have
$a=0.05$, $b=0.85$ and $c=0.2$, for (c) we have $a=0.4$, $b=0.8$ and $c=0.2$,
and for (d) we have $a=0.05$, $b=0.8$ and $c=0.2$.

The distribution of a Gompertz-Makeham random variable $Y$ can be perceived as
a mixture, i.e., its PDF (28) can be represented as a convex combination of
two probability distributions

$f(x)=\pi g_{1}(x)+(1-\pi)g_{2}(x),$ (29)

where

$g_{1}(x)=\frac{c}{\pi}\exp\left\\{-\frac{a}{b}\left(e^{bx}-1\right)-cx\right\\},\quad
g_{2}(x)=\frac{ae^{-(c-b)x}}{1-\pi}\exp\left\\{-\frac{a}{b}\left(e^{bx}-1\right)\right\\}$
(30)

and $\pi$ is given by (24).

Under the assumption of identity of the forces of mortality, the survival
function takes the form

$S(x)=S_{1}(x)\,S_{2}(x)=\pi G_{1}(x)+(1-\pi)G_{2}(x)$ (31)

with $G_{j}(x)=S(x\,|\,j)=\int\limits_{x}^{\infty}g_{j}(t)dt$ and
$S_{j}(x)=\exp\left\\{-\int_{0}^{x}h_{j}(t)dt\right\\}$, for $j=1,2$.

Figure 1 shows the decomposition of the Gompertz-Makeham PDF into senescent
($g_{2}$ density) and non-senescent ($g_{1}$ density) components, from which
we can see how non-senescent deaths shift the mode of the distribution (also
known as modal age at death) to the left, sometimes even leading to a second
peak of the distribution (panel (b) in Figure 1). Note that although the
competing risks are Gompertz and Exponential, (30) suggests that the Gompertz-
Makeham density is not represented by a convex combination of a Gompertz and
an Exponential distribution.

By representing the Gompertz-Makeham as a mixture of two distributions, we are
able to quantify the overall proportion of non-senescent deaths, given by the
quantity $\pi$, and also the proportion of non-senescent deaths at age $x$,
equal to $p(x)$. For the Makeham-Gompertz model, the age-specific proportion
of premature deaths can be expressed as $p(x)=c/\left(ae^{bx}+c\right)$.
Figure 2 shows the proportion of senescent and non-senescent deaths by age in
a Gompertz-Makeham model for different parameter combinations.

Figure 2: Proportion of senescent and non-senescent deaths by age in the
Gompertz-Makeham setting. For (a) we have $a=0.0005$, $b=0.1$ and $c=0.005$,
for (b) we have $a=0.05$, $b=0.85$ and $c=0.2$, for (c) we have $a=0.4$,
$b=0.8$ and $c=0.2$, and for (d) we have $a=0.05$, $b=0.8$ and $c=0.2$.

The function $p(x)$ aids in estimating the threshold age, $x^{*}$ that
separates the ages with predominant non-senescent deaths from the ages with
predominant senescent deaths. For the Gompertz-Makeham model,

$x^{*}=\left\\{\begin{array}[]{ l l }\frac{1}{b}\ln\frac{c}{a}&c\geq a\\\
0&c<a\\\ \end{array}\right.$ (32)

When $x^{*}=0$, $g_{1}(x)$ and $g_{2}(x)$ do not intersect (see, for example,
panel (c) in Figure 2). In all other cases, $x^{*}$ is their intersection
point, i.e., the point at which the proportion of senescent and non-senescent
deaths is equal.

##  5 Applications

Expressing Makeham models (1) as a mixture can be advantageous for analyzing
the components of mortality, e.g., the ones capturing premature and senescent
deaths. To illustrate these advantages, we estimate the Gompertz-Makeham model
using raw death counts and exposures after age 20 from the Human Mortality
Database (HMD, , 2023) for France, Italy, Japan, and Sweden, years 1947 to
2020, males and females separately.

Mazzuco et al.,  (2021) highlight that the share of premature deaths is also
defined by the shape, scale, and location of the senescent deaths
distribution. This suggests a recommendation to maintain consistent senescent
mortality parameters across countries while allowing for variability in
premature mortality parameters. However, we will not keep the senescent
mortality parameters constant for four main reasons: (i) the population in two
different years is comprised of individuals from different cohorts with
diverse life histories; (ii) distinct genetic backgrounds, shaped by
historical migrations and evolutionary processes, introduce variations in
susceptibility to age-related diseases, potentially influencing senescent
mortality; (iii) the interplay between genetic factors and environmental
determinants, including diet, lifestyle, pollution, and healthcare access,
amplifies differences in senescent mortality patterns among countries; and
(iv) disparities in healthcare systems and socioeconomic factors add
complexity, influencing disease prevention, diagnosis, and treatment
effectiveness among countries.

To estimate the parameters, we apply a Bayesian procedure (Gelman et al., ,
2013) and assume a Poisson distribution for the death counts $D_{ijk}$, where
the multi-index $ijk$ represents age $i$ in year $j$ for country $k$ (see, for
example, Brillinger, , 1986). The Bayesian estimates are obtained by the mode
of the posterior distributions, also known as the maximum a posteriori
probability (MAP) estimate (Patricio and Missov, , 2023).

The prior (and hyper-prior) distributions are defined as

$\displaystyle a_{ijk}|\alpha_{1},\beta_{1}$ $\displaystyle\sim$
$\displaystyle InverseGamma(\alpha_{1},\beta_{1});$ $\displaystyle
b_{ijk}|\alpha_{2},\beta_{2}$ $\displaystyle\sim$ $\displaystyle
InverseGamma(\alpha_{2},\beta_{2});$ $\displaystyle
c_{ijk}|\alpha_{3},\beta_{3}$ $\displaystyle\sim$ $\displaystyle
InverseGamma(\alpha_{3},\beta_{3});$ $\displaystyle\alpha_{1}$
$\displaystyle\sim$ $\displaystyle Gamma(2,2);$ $\displaystyle\alpha_{2}$
$\displaystyle\sim$ $\displaystyle Gamma(2,2);$ $\displaystyle\alpha_{3}$
$\displaystyle\sim$ $\displaystyle Gamma(2,2);$ $\displaystyle\beta_{1}$
$\displaystyle\sim$ $\displaystyle Gamma(1,1);$ $\displaystyle\beta_{2}$
$\displaystyle\sim$ $\displaystyle Gamma(1,1);$ $\displaystyle\beta_{3}$
$\displaystyle\sim$ $\displaystyle Gamma(1,1)\,,$

where $a$, $b$ and $c$ are the Gompertz-Makeham model parameters. Since the
parameters are strictly positive, we choose an inverse-gamma prior
distribution. The latter is characterized by a heavy tail and keeps
probability further from zero than the Gamma distribution.

By specifying the prior distributions, we employed the NUTS algorithm, a
variant of Hamiltonian Monte Carlo known as the No-U-Turn Sampler (Hoffman et
al., , 2014; Betancourt, , 2017), via the Rstan R-package (Team et al., ,
2016) to sample from the posterior distributions. Four chains were executed,
each consisting of 6000 iterations (4000 warm-ups and 2000 sampling). We
assessed convergence using the R-hat diagnostic and found that none of the
R-hat values exceeded 1.05 (see, for instance, Vehtari et al., , 2021). We
opted not to thin the chains, as thinning is typically unnecessary with the
efficient Hamiltonian Monte Carlo method employed by Stan. The point estimates
provided are MAP estimates, and the intervals displayed represent 95% Highest
Density Intervals (HPD intervals) for the posterior distributions.

Figure 3 presents the threshold age of senescent mortality within the
Gompertz-Makeham framework. In France, post-2000, the threshold age remains
relatively stable, fluctuating between approximately 30 and 35 years for males
and around 45 years for females. Almost the same holds for Japanese males
after 2010. For other populations, notably after the 1990s, the threshold age
exhibits a consistent linear increase. This trend suggests that in recent
years, senescent mortality prevalence has been postponed at an approximate
rate of 3 months per year for both Italian males and females, French males,
and Japanese females. Swedish females experience a delay of around 2 months
per year, while Swedish males see a more significant postponement of
approximately 4 months per year. Japanese males and French females show a
slower pace of postponement, at about half a month per year.

Figure 3: Threshold age between premature and senescence mortality estimated
through the Gompertz-Makeham model

In Figure 3, we observe a consistently low threshold age for males in France
(from 1960 to 1998) and Italy (from 1965 to 1985). During these periods, the
data show a complexity that the Gompertz-Makeham model, shown in Figure 4,
fails to capture. While the Gompertz-Makeham model is simplistic and fits well
the log-linear increase in the risk of dying, it struggles to accommodate non-
monotonic deviations from this pattern in the starting ages of analysis
(Missov and Németh, , 2015). This limitation leads to an underestimation of
the extrinsic risk parameter $c$. Consequently, while the Gompertz-Makeham
model adequately represents the data pattern for Sweden, it notably fails to
provide an accurate fit for France and Italy. This discrepancy is primarily
attributed to the complex data pattern observed between ages 20 and 45, which
deviates from both the log-linear increase in the risk of dying and the
assumption of constant premature mortality risk. As a result, the estimated
threshold age $x^{*}$ is simply the starting age of analysis (see Figure 3).

Figure 4: Fitted versus observed mortality rates for males from Sweden in
1965, France in 1966, and Italy in 1985.

The prevalence of premature mortality after age 20 decreases over time (see
Figure 5), faster until about 1970 and more moderately afterward, reaching an
almost constant level in recent years. Sex-specific premature mortality seems
to be balanced, with Sweden having a prevalence for males slightly higher than
the one for females over the entire period. For French and Italian males, we
observe a sudden increase from 1985 to 1990, leading to a persistent gap
between male and female premature mortality prevalence until the late 1990s
for France and about a decade later for Italy. This period coincides with the
emergence of the global HIV/AIDS epidemic, impacting, in particular, France
and Italy (see, for example Hamers et al., , 1998). During this period, both
countries experience a substantial rise in HIV/AIDS-related deaths. Meanwhile,
although Sweden and Japan were also affected by the epidemic, the prevalence
of HIV/AIDS was much lower compared to Italy and France (Gibney et al., ,
2006). Although advancements in HIV/AIDS treatment emerged, access to these
treatments varied. It was not until the later years of the 1990s that
antiretroviral therapy became more widely accessible in France and Italy
(Piot, , 2007; Ippolito et al., , 2001).

Figure 5: Overall prevalence of premature mortality after age 20 estimated
through the Gompertz-Makeham model

Despite the overall prevalence of premature mortality converging for some
countries and diverging for others, when we look at the age-specific
prevalence of premature mortality (Figure 6), we see a different picture. As
expected, the prevalence of premature mortality decreases with age. However,
for France, while the inter-sex difference in overall prevalence seems to
converge to zero, we do not observe this for the age-specific prevalence at
younger ages. The same trend is seen in Italy and Japan. For Sweden, however,
what we see is the opposite. The prevalence seems to be the same for males and
females between ages 30 and 45, and slightly higher for males from ages 60 to
75.

Figure 6: Age-specific prevalence of premature mortality estimated through the
Gompertz-Makeham model

Over time, the age-specific prevalence of premature mortality for the French
population seems to be stable for both sexes after 2000. The same holds for
Japanese males after 2010 and for Italian males after 1990. For Sweden, the
trend in age-specific prevalence of premature mortality seems to be identical
for both sexes. Other mortality measures, such as the premature and senescent
distributions of deaths, their force of mortality at ages 40, 60, and 80, and
their modal age at death, are presented in the appendix.

##  6 Conclusion

In this paper, we present a formal relationship that represents competing risk
models as a mixture of distributions, where each cause of death corresponds to
a component distribution, and an individual’s lifespan is a mixture of these
component distributions. The weights in this mixture correspond to the
probabilities of each cause of death being the first to occur among all
competing causes. This representation allows for a clear understanding of how
competing risks influence overall mortality and provides a framework for
modeling complex mortality patterns.

The mixture-model specification aids in representing the distribution of
deaths as a convex combination of distributions for risk-specific
subpopulations. This facilitates calculating various mortality and longevity
measures for each subpopulation, as well as assessing the overall and age-
specific prevalence of each cause of death.

We focus on the special case of Makeham mortality models that incorporate a
constant extrinsic risk of dying. We represent these models as a mixture that
reflects individual lifetimes in a competing-risk setting. The interpretation
is straightforward: an individual dies either according to a baseline
mortality mechanism or an exponential distribution.

In the case of a two-risk model, the Gompertz-Makeham, when components capture
senescent and premature deaths, we can estimate the threshold age that marks
the change of death prevalence from premature to senescent. We are also able
to reconstruct the age-specific profiles of the premature and senescent
mortality components. To illustrate these findings, we take advantage of a
Bayesian approach to estimate the Gompertz-Makeham model for the French,
Italian, Japanese, and Swedish populations from 1947 to 2020, after age 20.
The results suggest a postponement of senescent mortality prevalence at a pace
between 2 and 4 months per year for most of the populations included in our
analysis. The overall prevalence of premature mortality differs across age
groups and between sexes.

##  7 Acknowledgements

The research leading to this publication is part of a project that has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 research and innovation program (Grant agreement No.
884328 – Unequal Lifespans). Silvio C Patricio gratefully acknowledges the
support from the AXA Research Fund through funding the ”AXA Chair in Longevity
Research.”

## References

  * Beard, (1959) Beard, R. (1959).  Notes on some mathematical mortality models.  In Wolstenholme, G. and O’Connor, M., editors, The lifespan of animals, pages 802–811. Boston: Little, Brown.
  * Bebbington et al., (2007) Bebbington, M., Lai, C.-D., and Zitikis, R. (2007).  Modeling human mortality using mixtures of bathtub shaped failure distributions.  Journal of Theoretical Biology, 245(3):528–538.
  * Berman, (1963) Berman, S. M. (1963).  Note on extreme values, competing risks and semi-markov processes.  The Annals of Mathematical Statistics, 34(3):1104–1106.
  * Betancourt, (2017) Betancourt, M. (2017).  A conceptual introduction to hamiltonian monte carlo.  arXiv preprint arXiv:1701.02434.
  * Brillinger, (1986) Brillinger, D. R. (1986).  A biometrics invited paper with discussion: The natural variability of vital rates and associated statistics.  Biometrics, pages 693–734.
  * Cabral et al., (2022) Cabral, S., Castellares, F., and Queiroz, B. L. (2022).  Mortality modeling at old age: A mixture-model approach.  In PAA 2022 Annual Meeting. PAA.
  * Castellares et al., (2022) Castellares, F., Patrício, S., and Lemonte, A. J. (2022).  On the gompertz–makeham law: A useful mortality model to deal with human mortality.  Brazilian Journal of Probability and Statistics, 36(3):613–639.
  * Castellares et al., (2020) Castellares, F., Patrício, S. C., Lemonte, A. J., and Queiroz, B. L. (2020).  On closed-form expressions to gompertz–makeham life expectancy.  Theoretical Population Biology, 134:53–60.
  * Chiang, (1961) Chiang, C. L. (1961).  On the probability of death from specific causes in the presence of competing risks.  In Proc Fourth Berkely Symp Math Stat Probability, volume 4, pages 169–80.
  * Cox, (1959) Cox, D. R. (1959).  The analysis of exponentially distributed life-times with two types of failure.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 21(2):411–421.
  * Crowder, (1991) Crowder, M. (1991).  On the identifiability crisis in competing risks analysis.  Scandinavian Journal of Statistics, pages 223–233.
  * De Moivre, (1731) De Moivre, A. (1731).  Annuities upon lives, or, the valuation of annuities upon any number of lives, as also, of reversions: To which is added, an appendix concerning the expectations of life, and probabilities of survivorship.  London printed: and, Dublin re-printed, by and for Samuel Fuller.
  * Elandt-Johnson, (1976) Elandt-Johnson, R. C. (1976).  Conditional failure time distributions under competing risk theory with dependent failure times and proportional hazard rates.  Scandinavian Actuarial Journal, 1976(1):37–51.
  * Gail, (1975) Gail, M. (1975).  A review and critique of some models used in competing risk analysis.  Biometrics, pages 209–222.
  * Gelman et al., (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  CRC press.
  * Gibney et al., (2006) Gibney, L., DiClemente, R. J., and Vermund, S. H. (2006).  Preventing HIV in developing countries: biomedical and behavioral approaches.  Springer Science & Business Media.
  * Gompertz, (1825) Gompertz, B. (1825).  Xxiv. on the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. in a letter to francis baily, esq. frs &c.  Philosophical transactions of the Royal Society of London, (115):513–583.
  * Gumbel and Mustafi, (1967) Gumbel, E. J. and Mustafi, C. K. (1967).  Some analytical properties of bivariate extremal distributions.  Journal of the American Statistical Association, 62(318):569–588.
  * Hakulinen and Rahiala, (1977) Hakulinen, T. and Rahiala, M. (1977).  An example on the risk dependence and additivity of intensities in the theory of competing risks.  Biometrics, pages 557–559.
  * Hamers et al., (1998) Hamers, F. F., Downs, A. M., Infuso, A., and Brunet, J.-B. (1998).  Diversity of the hiv/aids epidemic in europe.  Aids, 12:S63–S70.
  * Heligman and Pollard, (1980) Heligman, L. and Pollard, J. (1980).  The age pattern of mortality.  Journal of the Institute of Actuaries, 107(1):49–80.
  * HMD, (2023) HMD (2023).  The human mortality database.  http://www.mortality.org/.
  * Hoffman et al., (2014) Hoffman, M. D., Gelman, A., et al. (2014).  The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.  J. Mach. Learn. Res., 15(1):1593–1623.
  * Ippolito et al., (2001) Ippolito, G., Galati, V., Serraino, D., and Girardi, E. (2001).  The changing picture of the hiv/aids epidemic.  Annals of the New York Academy of Sciences, 946(1):1–12.
  * Jodrá, (2013) Jodrá, P. (2013).  On order statistics from the gompertz–makeham distribution and the lambert w function.  Mathematical Modelling and Analysis, 18(3):432–445.
  * Kannisto, (1994) Kannisto, V. (1994).  Development of oldest-old mortality, 1950–1990: evidence from 28 developed countries.  Odense: Odense University Press.
  * Kozusko and Bajzer, (2003) Kozusko, F. and Bajzer, Ž. (2003).  Combining gompertzian growth and cell population dynamics.  Mathematical biosciences, 185(2):153–167.
  * Makeham, (1860) Makeham, W. M. (1860).  On the law of mortality and the construction of annuity tables.  Journal of the Institute of Actuaries, 8(6):301–310.
  * Makeham, (1867) Makeham, W. M. (1867).  On the law of mortality.  Journal of the Institute of Actuaries, 13(6):325–358.
  * Makeham, (1890) Makeham, W. M. (1890).  On the further development of gompertz’s law.  Journal of the Institute of Actuaries, 28(4):316–332.
  * Marshall and Olkin, (2007) Marshall, A. W. and Olkin, I. (2007).  Life distributions, volume 13.  Springer.
  * Mazzuco et al., (2018) Mazzuco, S., Scarpa, B., and Zanotto, L. (2018).  A mortality model based on a mixture distribution function.  Population Studies, 72(2):191–200.
  * Mazzuco et al., (2021) Mazzuco, S. S., Suhrcke, M. M., and Zanotto, L. L. (2021).  How to measure premature mortality? a proposal combining “relative” and “absolute” approaches.  Population health metrics, 19:1–14.
  * McLachlan et al., (2019) McLachlan, G. J., Lee, S. X., and Rathnayake, S. I. (2019).  Finite mixture models.  Annual review of statistics and its application, 6:355–378.
  * Missov and Németh, (2015) Missov, T. and Németh, L. (2015).  Sensitivity of model-based human mortality measures to exclusion of the makeham or the frailty parameter.  Genus, LXXI(2-3):113–135.
  * Missov and Németh, (2016) Missov, T. and Németh, L. (2016).  Sensitivity of model-based human mortality measures to exclusion of the makeham or the frailty parameter.  Genus, 71(2-3):113–135.
  * Missov et al., (2015) Missov, T. I., Lenart, A., Nemeth, L., Canudas-Romo, V., and Vaupel, J. W. (2015).  The gompertz force of mortality in terms of the modal age at death.  Demographic Research, 32:1031–1048.
  * Moeschberger and David, (1971) Moeschberger, M. and David, H. (1971).  Life tests under competing causes of failure and the theory of competing risks.  Biometrics, pages 909–933.
  * Norström, (1997) Norström, F. (1997).  The gompertz-makeham distribution.
  * Olshansky and Carnes, (1997) Olshansky, S. J. and Carnes, B. A. (1997).  Ever since gompertz.  Demography, 34(1):1–15.
  * Patricio and Missov, (2023) Patricio, S. and Missov, T. I. (2023).  Using a penalized likelihood to detect mortality deceleration.  Plos one, 18(11):e0294428.
  * Pereyra, (2019) Pereyra, M. (2019).  Revisiting maximum-a-posteriori estimation in log-concave models.  SIAM Journal on Imaging Sciences, 12(1):650–670.
  * Perks, (1932) Perks, W. (1932).  On some experiments in the graduation of mortality statistics.  Journal of the Institute of Actuaries, 63(1):12–57.
  * Piot, (2007) Piot, P. (2007).  The HIV pandemic: local and global implications.  OUP Oxford.
  * Plummer et al., (2006) Plummer, M., Best, N., Cowles, K., Vines, K., et al. (2006).  Coda: convergence diagnosis and output analysis for mcmc.  R news, 6(1):7–11.
  * Rainville, (1960) Rainville, E. D. (1960).  Special functions.  Macmillan.
  * Siler, (1979) Siler, W. (1979).  A competing-risk model for animal mortality.  Ecology, 60(4):750–757.
  * Siler, (1983) Siler, W. (1983).  Parameters of mortality in human populations with widely varying life spans.  Statistics in medicine, 2(3):373–380.
  * Tawn, (1988) Tawn, J. A. (1988).  Bivariate extreme value theory: models and estimation.  Biometrika, 75(3):397–415.
  * Team et al., (2016) Team, S. D. et al. (2016).  Rstan: the r interface to stan.  R package version, 2(1):522.
  * Thiele, (1871) Thiele, T. N. (1871).  On a mathematical formula to express the rate of mortality throughout the whole of life, tested by a series of observations made use of by the danish life insurance company of 1871.  Journal of the Institute of Actuaries, 16(5):313–329.
  * Vaupel et al., (1979) Vaupel, J., Manton, K., and Stallard, E. (1979).  The impact of heterogeneity in individual frailty on the dynamics of mortality.  Demography, 16:439–454.
  * Vehtari et al., (2021) Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Bürkner, P.-C. (2021).  Rank-normalization, folding, and localization: An improved $\hat{R}$ for assessing convergence of mcmc (with discussion).  Bayesian analysis, 16(2):667–718.
  * Zanotto et al., (2021) Zanotto, L., Canudas-Romo, V., and Mazzuco, S. (2021).  A mixture-function mortality model: illustration of the evolution of premature mortality.  European Journal of Population, 37(1):1–27.

## Appendix

Figure A-1: Estimated force of mortality at ages 40, 60, and 80 through the
Gompertz-Makeham model for premature and senescent groups Figure A-2:
Estimated senescence and overall modal age at death through the Gompertz-
Makeham model Figure A-3: Estimated non-senescent death distribution after age
20 through the Gompertz-Makeham model Figure A-4: Estimated non-senescent
death distribution through the Gompertz-Makeham model