New ways of estimating excess mortality of chronic diseases: Insights from the illness-death model
Ralph Brinks
q-bio.PE stat.ME





New ways of estimating excess mortality of chronic diseases: Insights from the illness-death model
    Ralph Brinks
Hiller Research Center for Rheumatology
University Duesseldorf
    
==================================================================================================




Recently, we have shown that the age-specific prevalence of a 
disease can be related to the transition rates in the illness-death model 
via a partial differential equation (PDE). In case of a chronic disease, 
we show that the PDE can be used to estimate excess mortality from prevalence 
and incidence. Applicability of the new method is demonstrated in a simulation and claims data
about diabetes in German men.




§ INTRODUCTION
 
Recently, we have shown that the age-specific prevalence of a 
health state or disease can be related to the transition rates in the illness-death model (IDM) 
via a partial differential equation (PDE) <cit.>. In case of a chronic disease, 
this relation can be used to estimate the incidence from a sequence of cross-sectional 
studies if information about mortality is available <cit.>. 

In this paper, we demonstrate that it is also possible to estimate excess mortality from prevalence 
and incidence of a chronic disease, which can be useful for the analysis of data from disease 
registers or health insurance claims. For this, we examine the relations of the illness-death model and
associated PDEs. In this context, we derive a new PDE which generalises 
the PDE of Brunet and Struchiner <cit.>. In a simulation study, the new PDE is used 
to demonstrate how the excess mortality can be estimated directly. Furthermore, we present
an estimation method in the framework of Bayesian statistics. In an application of the Bayesian approach, 
we estimate the excess mortality of diabetes from claims data comprising 70 million Germans.



§ ILLNESS-DEATH MODEL AND ASSOCIATED PARTIAL DIFFERENTIAL EQUATIONS
  
We consider the illness-death model for chronic (i.e., irreversible) diseases 
shown in Figure <ref>. The considered population is split 
into the relevant disease states Healthy (H) and Ill (I). From either states
people can transit into the state Dead (D). 
The transition rates between the three states are the incidence
rate (i), the mortality rate of the healthy (m_0) and the mortality rate 
of the diseased (m_1). These rates depend on the calendar time t and 
on the age a. Additionally, the mortality rate m_1 depends 
on the duration d of the disease. 



Let the numbers H and I denote the numbers of people in the 
respective states. To be more specific, H(t, a) is the number of healthy
people aged a,  a ≥ 0, at time t; and I(t, a, d) is the number of
diseased people aged a,  a ≥ 0, at t who are diseased for the duration 
d,  d ≥ 0. We assume that the considered population is sufficiently large
that H and I can be considered as smooth functions. The total number of subjects 
aged a at t who have the chronic 
disease is I^⋆(t, a) = ∫_0^a I(t, a, δ)  dδ.  



Let us furthermore assume that the considered population is closed, i.e., 
there is no migration and that the disease is contracted after birth. 
The later assumption implies I^⋆(t, 0) = 0 for all t. Then, we can formulate 
following equations for the change rates of H and I: 


    (∂_t + ∂_a)   H(t, a)    =    -  ( m_0(t, a) + i(t, a)  )   H(t, a) 
    
    (∂_t + ∂_a + ∂_d)   I(t,a,d)    =     - m_1(t, a, d)   I(t, a, d),


where ∂_x means the partial derivative with respect to 
x, i.e., ∂_x = ∂∂ x for x ∈{t, a, d}.

In addition to the PDEs (<ref>) – (<ref>)), which describe the outflows from the 
states Healthy and Ill, we need the following initial conditions:

    H(t - a, 0)     = H_0(t - a). 
    
    I(t,a, 0)     = i(t, a)   H(t, a).

The first initial condition represents the number of (disease-free) newborns H_0, 
and the second initial condition describes the number of newly diseased persons, the incident cases.



In <cit.>, we have shown that the age-specific prevalence 
p(t, a) = I^⋆(t, a)I^⋆(t, a) + H(t, a) is the solution of a scalar PDE
that can be derived from the two-dimensional system (<ref>) – (<ref>). Here,
we choose a different approach. Instead of considering the age-specific prevalence p(t, a), 
we follow the idea of Brunet and Struchiner and 
examine the prevalence-odds π(t, a) = I^⋆(t, a)H(t, a) <cit.>. 
Using the terminology ∂ = ∂_t + ∂_a, we obtain 

    ∂π   =   ∂ I^⋆/H + π  ( i + m_0  ) 
       =   -m_1^⋆  I^⋆ + i   H/H + π  ( i + m_0  ) 
       =    i + π  ( i + m_0 - m_1^⋆ ).


For the second equality we used ∂ I^⋆ = - m_1^⋆  I^⋆ + i   H, which has been
proven in the Appendix of <cit.>. The rate m_1^⋆ is defined as

    m^⋆_1(t, a) := 
    ∫_0^a  m_1(t, a, δ)   I(t, a, δ) dδ/∫_0^a  I(t, a, δ) dδ    for  I^⋆(t, a) > 0
    
    0     for  I^⋆(t, a) = 0.





  
Remark: The rate m_1^⋆ may be accessible in epidemiological surveys
by choosing a sample population with representative distribution of disease duration.
However, in most practical cases, it is unknown because the distribution 
I(t, a, d)∫_0^a  I(t, a, δ) dδ in Eq. (<ref>)
is not known.



Thus, we obtain following linear scalar PDE

    ∂π = i - π  ( i + m_0 - m_1^⋆ ),


which shows how the temporal change of the prevalence-odds
∂π is governed by the rates in the illness-death model in Figure <ref>
and the value of the prevalence-odds itself.



  
Remark: If m_1 does not depend on the duration d, we have m_1^⋆ = m_1
and Eq. (<ref>) is equivalent to the PDE (3) in Brunet and Struchiner <cit.>. Hence, Eq. (<ref>)
is a generalisation of the PDE of Brunet and Struchiner.



  
Remark: Eq. (<ref>) is equivalent to 

    ∂ p = (1-p)   ( i - p  (m_1^⋆ - m_0  )  ),

which has been proven in <cit.>. With the definition R = m_1^⋆m_0
and letting m be the overall mortality, m = p   m^⋆_1 + (1 - p)   m_0, 
then Eq. (<ref>) becomes

    ∂ p = (1-p)   ( i - m  p   (R - 1)/1 + p   (R - 1) ).




For our purpose of estimating the excess mortality Δ m = m_1^⋆ - m_0, 
Eq. (<ref>) is very useful, because it holds

    Δ m = i   (1 + π) - ∂π/π.




An advantage of the approach of Brunet and Struchiner lies in an explicit representation 
of the prevalence-odds in case the rates i, m_0 and m_1 are given. Then starting
from Eq. (<ref>) – (<ref>) combined with the initial conditions of above, 
we obtain following equation by using calculus:


    π(t, a)    = ∫_0^a i(t-δ, a-δ)  ×
           exp ( - ∫_0^δ{ m_1(t - δ + τ, a - δ + τ, τ) 
                                                 - (i + m_0)  (t - δ + τ, a - δ + τ ) }dτ ) dδ.


With π = p1-p, we see that Eq. (<ref>) is a generalisation of 
Eq. (1) in <cit.>. One advantage of the explicit representation of 
π in (<ref>) is the possibility to (numerically) calculate π with a 
prescribed accuracy, e.g. by Romberg integration <cit.>, which we will use 
in the examples below. Numerical solutions
of differential equations usually do not allow prescribed levels of accuracy. 






§ EXAMPLES AND DEMONSTRATION



 §.§ Direct estimation of excess mortality

The first example is about a hypothetical chronic disease with all time-scales 
t, a and d playing a role. The incidence of the chronic disease is assumed to 
be i(t, a) = (a-30)_+3000, which implies that the
disease affects only people aged 30 and older. 
The age-specific mortality rate of the non-diseased is 
chosen to be m_0(t,a) = exp(-10.7 + 0.1 a + t ln(0.98)). 
In addition, we assume that the mortality m_1 of the diseased can be written as
a product of m_0 and a factor that depends only on the duration d:

    m_1(t, a, d) = m_0(t, a) × ( 0.04 (d - 5)^2 + 1).

Except for the time trend in m_0, this example is the same as Simulation 2 in 
<cit.>.

For the example, we mimic the situation that we have three cross-sectional studies in the years
t_0 = 95, t_1 = 100, and t_2 = 105. We calculate the prevalence odds π for these years via
Eq. (<ref>). Figure <ref> shows the prevalence-odds for the three years. 
Until age of about 70 years, the three prevalence-odds are virtually the same. For our example, we
additonally assume that we have the age-specific incidence rate available for the year t_1 = 100.
The aim is to estimate the excess mortality in t_1.




The proposed method to estimate the excess mortality Δ m in the year t_1=100 
is direct application of Eq. (<ref>). As assumed the incidence i for t_1 = 100
is assumed to be given. The partial derivative ∂π is approximated by 
following finite difference:

    ∂π (t_1, a) ≐1t_2-t_0  [ π(t_2, a + t_2-t_02) - 
                π(t_0, a - t_2-t_02)  ].


Then, we the excess mortality Δ m can be estimated by plugging these numbers into
Eq. (<ref>). In case the mortality rate
m_0 of the non-diseased is known, Δ m is often expressed in terms of the hazard ratio 
R = m_1^⋆m_0 which can be obtained from

    R = 1 + Δ m/m_0.


The age-specific HR expresses the mortality rate of the diseased people relative to the non-diseased
at the same age. For the hypothetical chronic disease we find the age-specific HR as in Figure <ref>.
The age-specific HR is peaking between age a = 70 and a = 80 and falling with increasing age.



In case, the mortality rate m_0 of the non-diseased population is not known, Δ m can also be
compared to the mortality rate m of the general population. It holds 
m = (1-p)   m_0 +  p   m_1^⋆ = 11+π  m_0 +  π1+π  m_1^⋆. 
Usually, the mortality of the general population is accessible from vital statistics of the federal 
statistical offices.



 §.§ Bayes estimation of excess mortality

The second example is about claims data from Germany during the years 2009 to 2015. Goffrier 
and colleagues reported the age-specific
prevalence of diabetes of German men in the years t_0 = 2009 and t_1 = 2015 as shown in Figure 
<ref>, <cit.>.



Based on the incidence rate (i) in 2012 (reported in Table 5 of <cit.>), our aim is to estimate
the age-specific hazard ratio R(a) for the same year. For this, we use a Bayes approach. Motivated 
by empirical findings from the Danish diabetes register, we assume that the logarithm of the age-specific
HR R approximately is a straight line in the age range 50 to 90 years of age (see Figure 5 in <cit.>). 
Thus, we make the approach

    ln R(a) = ln R(50) + ln R(90) - ln R(50)/90 - 50× a.


For R(50) and R(90) we use weakly informative prior distributions R(50) ∼ U(1.5, 4) 
and R(90) ∼ U(1, 2.5), where U means the uniform distribution.
In Bayes terminology, our aim is to estimate a-posteriori distributions for 
R(50) and R(90). 

This done by randomly drawing R(50) and R(90) from the prior distributions, solving
the PDE (<ref>) with initial condition p(2009, a). For solving the PDE,
we use the Method of Characteristics <cit.> to convert Eq. (<ref>)
into an ordinary differential equation (ODE) and then solve the ODE by the Runge-Kutta Method
of fourth order <cit.>.
The calculated prevalence in 2015, p(2015, a), is then compared with the 
observed prevalence in 2015 as shown as blue line in Figure <ref>.
Instead of the joint a-posteriori distribution, the log-likelihood of the deviation between observed
and calculated prevalence is computed. The results are shown in Figure <ref>. 
The black cross indicates the maximum a-posteriori 
(MAP) estimator for these data, which is given by R_MAP(50) = 2.19 and 
R_MAP(90) = 1.46.







§ DISCUSSION
 
In this work, we have described how the illness-death model can be
used to obtain information about mortality in case prevalence and incidence are given. 
This allows insights into the excess mortality of people with chronic diseases compared
to the people without the disease or the general population. 

We presented two methods of estimating excess mortality, one direct estimation method and one
method in the context of Bayesian statistics. The methods can be used
if prevalence and incidence of chronic disease are given, e.g. in claims data <cit.> 
or in the setting of disease registers <cit.>.















Contact: 

Ralph Brinks 

Hiller Research Center for Rheumatology

University Clinics Duesseldorf

D-40225 Duesseldorf


