Estimation of the transition rates in the illness-death model for chronic diseases from aggregated current status data: a feasibility and simulation study

Recently, it has been shown that the transition rates of the illness-death model (IDM) for chronic conditions are related to the age-specific prevalence by a partial differential equation (PDE). Given mortality, the PDE could be used to estimate incidence rates from cross-sectional data. The aim of this article is to extend the IDM and introduce a novel method to estimate the age-specific incidence rate together with the two mortality rates from aggregated current status (ACS) data. By ACS data we mean counts of people in the four states of the extended IDM at different points in time. ACS data stem from epidemiological studies where only current disease status and vital status data need to be collected without following-up people (as, for example, in cohort studies). To demonstrate feasibility of the method, we use a simulation study from the context of diabetes in Germany. Two estimation methods are introduced, a least squares estimator and a maximum likelihood estimator. We find a good agreement between the estimates and the input parameters used to set up the simulation.

Introduction

It could recently be shown that the transition rates i, m₀, and m₁ in the illness-death model (IDM, see left part of Figure 1) for chronic diseases are linked to the age-specific prevalence of the disease via a differential equation (DE) (1). Given information about mortality rates, the DE could be used to estimate the incidence from the prevalence (2) or to project future case-numbers of chronic diseases based on assumptions about the future transition rates (“scenarios”) (3). In the case that prevalence and incidence of the chronic condition are given, measures of excess mortality could be derived from the DE (4), although diagnostic accuracy turned out to be very important (5).

Figure 1

Figure 1. Two illness-death models (IDMs) for chronic diseases with transition rates i, m₀, and m₁ between the states. The conventional IDM (left) consists of the three states Non-diseased, Diseased and Dead. The extended IDM (right) distinguishes between Dead without disease and Dead with disease.

The aim of this article is to extend the IDM by an additional state and derive a system of DEs to estimate all three rates, i, m₀, and m₁ from aggregated current status (ACS) data. ACS data is closely connected to a state model (such as the IDM) and the information that we have available is the number of people in each of the states of the associated state model at specific points in time. The points in time are called monitoring times (6), and at each of the monitoring times, we only have the total counts of people in the health states. This is advantageous in settings with high data protection standards like, for example, in health insurance claims data [see e.g., (2)].

Methods

Extended IDM and system of DEs

We extend the conventional IDM as depicted in the left part of Figure 1 by separating the state Dead in Figure 1 into two different (disjunct) states, Dead without disease and Dead with disease. This means, we differentiate whether a deceased subject has contracted the disease before death or not. Hence, each subject of the population under consideration is assigned to exactly one the four states shown in the right part of Figure 1: Non-diseased (state 1 in Figure 1), Diseased (2), Dead without disease (3) and Dead with disease (4). The numbers, i.e., counts, of people in the respective states 1–4 are denoted by S, C, D₀, and D₁. Subjects may change their state along the arrows in Figure 1 as time t evolves.

Simulation setup and estimation

To apply the theory in a realistic setting, we choose a test example from (7) motivated by the situation of type 2 diabetes in Germany. Diabetes is assumed to be irreversible. In the example, we have only one time scale and the time variable t, which is interpreted as age (7). For setting up the simulation, the age-specific incidence rate i is given by i(t) = max(0, t−30)/2,000 and the mortality rates m₀ and m₁ are assumed to be of Gompertz-type with m₀(t) = exp(−10.7 + 0.1 t) and m₁(t) = exp(−10 + 0.1 t). The simulation starts with a population of N = 10,000 study participants, for whom the changes of the states in the extended IDM as in the right part of Figure 1 are simulated by a micro-simulation (8). In short, for each subject possible transition times between the states in the IDM are randomly drawn by inverse transform sampling. Then, we mimic a sequence of eleven cross-sectional studies at time points t₁ = 0, t₂ = 10, …, t₁₁ = 100. Each subject is randomly allocated to exactly one of the eleven cross-sections t_k, where the vector X(t_k) = (S(t_k), C(t_k), D₀(t_k), D₁(t_k)) of count data, is determined for k = 1,…, 11.

We assume that the functional form of the rates is known (incidence rate and Gompertz mortality rates) and the aim of our simulation study is to estimate the numeric parameters underlying the rates. This means, in statistical terms we have a six-dimensional parametric model [Definition 2.1 (9)] with parameters ϑ = (ϑ₁, ϑ₂, ϑ₃, ϑ₄, ϑ₅, ϑ₆), which we want to estimate from the simulated ACS data. The associated rates read as i(t) = ϑ₂ × max(0, t−ϑ₁), m₀(t) = exp(ϑ₃ + t × ϑ₄) and m₁(t) = exp(ϑ₅ + t × ϑ₆). The parameter ϑ₁ is the age of onset of diabetes, which for type 2 diabetes in Germany is at the age of about 30 years. The parameter ϑ₂ is the slope of the age course of the incidence rate. As a coarse approximation of the age-specific incidence rate in Germany, we choose an incidence rate that increases linearly with age starting at age 30. The remaining parameters ϑ₃ to ϑ₆ describe the two Gompertz mortality rates. The true parameter vector ϑ^(true) underlying the simulation is ϑ^(true) = (ϑ₁, ϑ₂, ϑ₃, ϑ₄, ϑ₅, ϑ₆) = (30, 1/2000, −10.7, 1/10, −10, 1/10). We will use the findings of the seminal work of (10) to derive a system of DEs that describe the extended IDM. The system of DEs will be the ground for two estimation methods for ϑ from the ACS data: the first is based on a least squares fit, the second is a maximum likelihood approach. The estimated ϑ for both methods will be compared to the true ϑ^(true) in terms of the absolute relative error.

The source code of the simulation and both estimation approaches in the free statistical programming language R (The R Foundation for Statistical Computing) is available from the open access public repository Zenodo under DOI 10.5281/zenodo.16789935 (11).

Results

Extended IDM and differential equations

We will use the extended IDM as in the right part of Figure 1. Let S(t), C(t), D₀(t), and D₁(t) denote the number of subjects in the respective state at time t. The non-negative transition rates may depend on time t and the fractions of people in the four states (numbered 1–4 in Figure 1) are denoted by p_j = p_j(t), j = 1, 2, 3, 4 according to the numbers of the states. For example, at some time t* the fraction p₃(t*) denotes the fraction of people, who are deceased without disease at t*. The fractions p₁, p₂, p₃, and p₄ always add up to 100%: p₁(t) + p₂(t) + p₃(t) + p₄(t) = 1 for all t. Furthermore, let N(t) = S(t) + C(t) + D₀(t) + D₁(t) and for brevity define the integer valued, four dimensional vector X(t): = (S(t), C(t), D₀(t), D₁(t)).

The work of Jahnke & Huisinga (10) from the context of state models proves that X(t) follows a multinomial distribution where the parameter p = p(t) of the multinomial distribution is the solution of the following system of ordinary differential equations (ODEs):

${\displaystyle \frac{d}{dt}p(t)=p^{\mathrm{\prime }}(t)=A(t)p(t)}(1)$

with initial condition p(0) = p₀. The matrix A(t) in Equation (1) is given by

$A(t)=\left[\begin{array}{cccc}-i(t)-m_0(t)& 0& 0& 0\\ i(t)& -m_1(t)& 0& 0\\ m_0(t)& 0& 0& 0\\ 0& m_1(t)& 0& 0\end{array}\right]$

Using p(t) = (p₁(t), p₂(t), p₃(t), p₄(t)) and omitting the time dependency, we obtain the following linear system of ODEs:

$p_1^{\mathrm{\prime }}=-(i+m_0)p_1(\mathrm{2a})$

$p_2^{\mathrm{\prime }}=i{p}_1-m_1{p}_2(\mathrm{2b})$

$p_3^{\mathrm{\prime }}=m_0{p}_1(\mathrm{2c})$

$p_4^{\mathrm{\prime }}=m_1{p}_2.(\mathrm{2d})$

Note that p(t) ∈ (0,1)⁴ is four-dimensional, but p₁(t), p₂(t), p₃(t), p₄(t) ∈ (0,1) are scalars.

We will use the four-dimensional system as in Equation (1) or the four equivalent scalar Equations (2a–d) to estimate the parameter ϑ underlying the rates i, m₀, and m₁ from the four-dimensional count vector X(t) = (S(t), C(t), D₀(t), D₁(t)).

Simulation

After running the micro-simulation detailed in (8) with the rates i, m₀ and m₁ parameterized with ϑ^(true) as described above, we randomly assign the N = 10,000 simulated individuals to exactly one of the cross-sections at t_k = 0, …, 100. Table 1 shows the resulting number of subjects in each state at each of the eleven cross-sections. For example, at cross-section t_k = 60 a total of N₇ = 899 were counted, with 677, 143, 55, and 24 in the states 1–4 (see the numbers in Figure 1), respectively. The total number N_k of participants at each of the eleven cross-sections ranges from 877 to 959. As the sum over all S(t_k) equals 6,784, we see that the huge majority of sampled individuals are in health state 1, Non-diseased.

Table 1

Table 1. Aggregated current status data at each of the eleven cross-sections.

The resulting fractions p^(obs) are shown in Figure 2 as filled dots. For comparison, the associated components of the solution p of the ODE system (1) with the true parameters ϑ^(true) for i, m₀ and m₁ are shown as lines. Deviations of the dots from the solution show the random character of the ACS data [which come from a multinomial distribution (10)].

Figure 2

Figure 2. The four components of p^(obs) from the ACS data in Table 1 are shown as filled dots. For comparison, the associated components of the solution p of Equation (1) are shown as lines. The colors refer to the components, for example, green presents the third of the four components.

We use the system given in Eqs. (2a–d) to determine the parameters ϑ in the rates i, m₀ and m₁ from ACS data as shown in Table 1 and graphically depicted as dots in Figure 2.

Least squares and maximum likelihood estimation

For a given parameter vector ϑ, we can solve Equations (2a–d) and obtain the solution p(t; ϑ). Then for the first estimation method, we may calculate the squared residual || p(; ϑ)−p^(obs) ||² of the difference between p(; ϑ) and the observed fractions p^(obs) in the ACS data. The norm || · ||² sums over all observed points in time t_k, k = 1, …, K; for example in Table 1 from t₁ = 0 to t₁₁ = 100 in units of 10. The optimal ϑ* in the sense of least squares is given by

${\vartheta }^{\ast }=\arg \min ||p(;\vartheta )-p^{(\mathrm{obs})}|{|}_{.}^2(3)$

In other words, ϑ* minimizes the squared difference between the modelled solution p(t; ϑ) of system Equations (2a–d) and the observed fractions p^(obs) ∈ (0,1)⁴; hence, ϑ* is a least squares estimator [Definition 3.1 in (9)].

In case that study participants at the K cross-sections are all independent, the log-likelihood function ℓ(ϑ) of the multinomial distribution can be used to estimate the unknown parameter ϑ = (ϑ₁, ϑ₂, ϑ₃, ϑ₄, ϑ₅, ϑ₆). Under independence we obtain

$\ell (\vartheta )=\sum _{k=1}^K\left[log(n_k!)+\sum _{j=1}^4\{x_{k,j}log(p_j(t_k;\vartheta ))-log(x_{k,j}!)\}\right](4)$

which can be maximized with respect to ϑ. This provides a second estimation method, a maximum likelihood estimator [Definition 4.3 in (9)]. In Equation (4) x_k,j denotes the number of subjects in state j at time t_k, for j = 1, 2, 3, 4 and k = 1, …, K.

To sum up, for estimating ϑ we can use two different estimation methods: the least-squares minimization as in Equation (3) and the maximum likelihood estimation as in Equation (4). In both estimation approaches, we need to calculate the solution p(t; ϑ) for a given ϑ, which will be done by the classical Runge-Kutta method of fourth order [rk4 in the R package deSolve (12)]. The optimization method optim of the statistical software R provides the least squares estimates and the maximum likelihood estimates according to Equations (3, 4), respectively. Table 2 shows the results for both estimation methods based on the data in Table 1 together with the absolute relative error.

Table 2

Table 2. Least squares and maximum likelihood estimates for the parameter ϑ with absolute relative errors (as percentages in brackets).

It is visible from Table 2 that both estimation methods provide estimates with an absolute relative error below 10%. It seems that the maximum likelihood estimates are slightly superior to the least squares estimates, because the absolute relative errors of the maximum likelihood method are consistently below 4%.

Discussion

In this article, we have sketched how aggregated current status (ACS) data from the extended illness-death model (IDM) can be utilized to obtain information about the transition rates, the incidence (i) and the mortality rates (m₀ and m₁). IDM. We showed the feasibility of two novel estimation methods based on a system of ordinary differential equations (ODE): a least squares estimator and a maximum likelihood estimator. The least squares estimator uses the prevalence proportions p = (p₁, p₂, p₃, p₄), while the maximum likelihood estimator requires count data X = (S, C, D₀, D₁) (as in Table 1). The slight superiority of the maximum likelihood estimates compared to the least squares estimates in terms of absolute relative errors can be seen in the light of statistical efficiency. In short, the maximum likelihood estimation is the standard for asymptotic efficiency because it fully utilizes the distributional assumption of the data [see Theorem 6.2.2 in (13)]. The least squares estimator does not require a full distribution but therefore sacrifices efficiency. In case of the maximum likelihood estimation, 95% confidence bounds can be obtained with standard arguments from asymptotic statistics via inversion of the Fisher information matrix [Section 4.3 from (14)]. Unfortunately, the least squares estimator does not allow such an easy derivation of secondary order statistics. The least squares estimator requires further considerations for inferential measures, for example multidimensional resampling (15). However, secondary order statistics are beyond the scope of this feasibility study.

We think that the novel methods will be helpful in settings where follow-up data are not available. To estimate transition rates, usually cohort studies are run. ACS data, however, allude to a study design, where following-up of study participants is not necessary. Compared to simple prevalence data π: = p₂/(p₁ + p₂), which is frequently reported in epidemiological studies, for the novel method described here, the number or percentage of deceased persons is necessary. These data about fatalities may be obtained from official residents’ registries (in case they exist). In case only the simple prevalence π is reported, we cannot estimate all rates of the IDM. Inserting Eqs. (2a,b) yields the ODE πʹ = (1−π)[i−π (m₁−m₀)], which compared to the four dimensional system (1) has the advantage of being scalar (making computations easier), but can only be used to estimate one rate from π (and πʹ) given the other two rates.

The possibility of using one study for estimating all three rates (incidence i, moralities m₀ and m₁) has an advantage for considering stochastic dependencies between the prevalence and the three rates. When two (or more) data sources are used to estimate measures of the IDM like e.g., in (2), usually these dependencies cannot be accounted for, because figures of covariance between prevalence and rates are not provided among the different data sources. If only one data source is used, it is possible to consider dependencies e.g., in multidimensional resampling (15) or bootstrapping (16).

The question arises, where ACS data play a role in practical applications. ACS data has the advantage that in many cases individual study subjects are non-identifiable, because data are aggregated and need not be reported on the individual level. This can be important for data protection reasons—especially if the health state Ill refers to particularly sensible diseases like sexually transmitted diseases. As an example for an application of ACS data, one might think of people with permanent need for long-term care in elderly population. In Germany, the mortality rate ratio (MRR) for these people is unknown. Follow-up studies are difficult because elderly people can be difficult to contact, withdrawal of consent is frequent and transport to examination centers imposes more problems than in a younger study population. Studies, which we think are suitable for application of the proposed method, are for example the Survey of Health, Ageing and Retirement in Europe (SHARE) and the Health and Retirement Study (HRS). Both studies provide large samples in an older population and collect comprehensive information about chronic conditions and mortality. Aggregated count data similar to Table 1 can be extracted after taking the survey weighting scheme into account. However, the practical aspects with weighting and secondary order statistics necessary for an appropriate analysis of these studies is beyond the scope of this feasibility study.

Usually when a study participant dies, the subject leaves the survey (“drop out”) and data is used only up to that point in time. In the proposed method, this is different: data can enter the analysis many years later after the fatality occurred. To illustrate this, we provide details of the most extreme case of the simulation described above. Study subject 3,512 (of the N = 10,000) dies young without contracting the disease at age 19.2 years. The data of subject 3,512 enter the ACS data of Table 1 more than 80 years later at the cross-section at t_k = 100. The subject is one of the subjects at the cell D₀(100) = 488 in Table 1 above, which means that the statistical information of this young fatality is used no earlier than 80 years after death. The fact that information about the state (death) can be used for a data point more than 80 years later, requires that the health state does not change in the meantime (which is obvious in case of death). In case of incorrect assignments to health states, the estimates may become biased (misclassification). An example may be if a subject is counted towards D₀ although the subject contracted the disease (and should hence be counted toward D₁). This situation has been called misclassification of disease state at death (17).

The question arises how the method described here can be generalized to the case, when more than one time scale is involved. Our micro-simulation can be considered as a birth cohort, which implies that one time variable (here t) is sufficient to capture all relevant temporal changes. In epidemiological practice, frequently two times scales, calendar time and age, are involved, see for example (18). In this case, the ODE system in Equation (1) will become a partial differential equation (PDE) using the multidimensional chain rule [as described in (19)]. Since the PDE can be reduced to an ODE similar to Equation (1) by the method of characteristics (20), the techniques for treating two time scales are very similar to the methods described in this manuscript.

To sum up, we could show that the extended IDM can be used to obtain insights into all three transition rates. Two methods for estimating the parameters have been developed, one method is based on a least squares minimization and the other is based on maximum likelihood optimization. Applicability has been demonstrated in a simulation study motivated by diabetes in Germany. We find a good agreement between both estimation methods and the input used to set up the simulations.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants legal guardians/next of kin in accordance with the national legislation and the institutional requirements.

Author contributions

RB: Writing – original draft, Resources, Investigation, Formal analysis, Software, Visualization, Data curation, Project administration, Conceptualization, Validation, Writing – review & editing, Methodology, Supervision. MM: Writing – review & editing. SV: Writing – review & editing.

Funding

The author(s) declared that financial support was not received for this work and/or its publication.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The authors RB, SV declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Abbreviations

ACS, aggregated current status; DE, differential equation; IDM, illness-death model; ODE, ordinary differential equation; PDE, partial differential equation.

References

1. Brinks R, Landwehr S. A new relation between prevalence and incidence of a chronic disease. Math Med Biol. (2015) 32(4):425–35. doi: 10.1093/imammb/dqu024

PubMed Abstract | Crossref Full Text | Google Scholar

2. Tamayo T, Brinks R, Hoyer A, Kuß OS, Rathmann W. The prevalence and incidence of diabetes in Germany. Dtsch Arztebl Int. (2016) 113(11):177–82. doi: 10.3238/arztebl.2016.0177

PubMed Abstract | Crossref Full Text | Google Scholar

3. Tönnies T, Röckl S, Hoyer A, Heidemann C, Baumert J, Du Y, et al. Projected number of people with diagnosed type 2 diabetes in Germany in 2040. Diabet Med. (2019) 36(10):1217–25. doi: 10.1111/dme.13902

PubMed Abstract | Crossref Full Text | Google Scholar

4. Tönnies T, Hoyer A, Brinks R. Excess mortality for people diagnosed with type 2 diabetes in 2012—estimates based on claims data from 70 million Germans. Nutr Metab Cardiovasc Dis. (2018) 28(9):887–91. doi: 10.1016/j.numecd.2018.05.008

Crossref Full Text | Google Scholar

5. Brinks R, Tönnies T, Hoyer A. Impact of diagnostic accuracy on the estimation of excess mortality from incidence and prevalence: simulation study and application to diabetes in German men. F1000Res. (2021) 10:49. doi: 10.12688/f1000research.28023.1

PubMed Abstract | Crossref Full Text | Google Scholar

6. Jewell NP, van der Laan M. Current Status data: review, recent developments and open problems. Handbook of Statistics 23. Amsterdam: Elsevier (2003). p. 625–42. doi: 10.1016/S0169-7161(03)23035-2

Crossref Full Text | Google Scholar

7. Brinks R. Illness-death model in chronic disease epidemiology: characteristics of a related, differential equation and an inverse problem. Comput Math Methods Med. (2018):5091096. doi: 10.1155/2018/5091096

PubMed Abstract | Crossref Full Text | Google Scholar

8. Brinks R, Landwehr S, Fischer-Betz R, Schneider M, Giani G. Lexis diagram and illness-death model: simulating populations in chronic disease epidemiology. PLoS One. (2014) 9(9):e106043. doi: 10.1371/journal.pone.0106043

PubMed Abstract | Crossref Full Text | Google Scholar

9. Shao J. Mathematical Statistics. New York: Springer (2003).

Google Scholar

10. Jahnke T, Huisinga W. Solving the chemical master equation for monomolecular reaction systems analytically. J Math Biol. (2007) 54:1–26. doi: 10.1007/s00285-006-0034-x

PubMed Abstract | Crossref Full Text | Google Scholar

11. Brinks R. Example for Estimation of Rates in Illness-death model for Chronic Diseases Based on Aggregated Current status Data. Zenodo (2025). doi: 10.5281/zenodo.16789935

Crossref Full Text | Google Scholar

12. Soetaert K, Petzoldt T, Setzer RW. Solving differential equations in R: package deSolve. J Stat Softw. (2010) 33(9):1–25. doi: 10.18637/jss.v033.i09

PubMed Abstract | Crossref Full Text | Google Scholar

13. Hogg RV, McKean JW, Craig AT. Introduction to Mathematical Statistics. 8th ed. Harlow: Pearson (2020).

Google Scholar

14. Wood S. Core Statistics. Cambridge: Oxford University Press (2015).

Google Scholar

15. Oakley JE, O'Hagan A. Probabilistic sensitivity analysis of complex models: a Bayesian approach. J R Stat Soc Ser B Stat Methodol. (2004) 66(3):751–69. doi: 10.1111/j.1467-9868.2004.05304.x

Crossref Full Text | Google Scholar

16. Efron B, Tibshirani RJ. An introduction to the Bootstrap. New York: Chapman and Hall/CRC (1994).

Google Scholar

17. Voß S, Hoyer A, Landwehr S, Pavkov ME, Gregg E, Brinks R. Estimation of mortality rate ratios for chronic conditions with misclassification of disease status at death. BMC Med Res Methodol. (2024) 24(1):2. doi: 10.1186/s12874-023-02111-3

PubMed Abstract | Crossref Full Text | Google Scholar

18. Carstensen B, Rønn PF, Jørgensen ME. Prevalence, incidence and mortality of type 1 and type 2 diabetes in Denmark 1996–2016. BMJ Open Diabetes Res Care. (2020) 8:e001071. doi: 10.1136/bmjdrc-2019-001071

PubMed Abstract | Crossref Full Text | Google Scholar

19. Brinks R, Hoyer A. Illness-death model: statistical perspective and differential equations. Lifetime Data Anal. (2018) 24(4):743–54. doi: 10.1007/s10985-018-9419-6

PubMed Abstract | Crossref Full Text | Google Scholar

20. Olver PJ. Introduction to Partial Differential Equations. New York: Springer (2014).

Google Scholar