Evaluating Semi-Markov Processes and Other Epidemiological Time-to-Event Models by Computing Disease Sojourn Density as Partial Differential Equations

Abstract

Introduction

Epidemiological models benefit from incorporating detailed time-to-event data to understand how disease risk evolves. For example, decompensation risk in liver cirrhosis depends on sojourn time spent with cirrhosis. Semi-Markov and related models capture these details by modeling time-to-event distributions based on published survival data. However, implementations of semi-Markov processes rely on Monte Carlo sampling methods, which increase computational requirements and introduce stochastic variability. Explicitly calculating the evolving transition likelihood can avoid these issues and provide fast, reliable estimates.

Methods

We present the sojourn time density framework for computing semi-Markov and related models by calculating the evolving sojourn time probability density as a system of partial differential equations. The framework is parametrized by commonly used hazard and models the distribution of current disease state and sojourn time. We describe the mathematical background, a numerical method for computation, and an example model of liver disease.

Results

Models developed with the sojourn time density framework can directly incorporate time-to-event data and serial events in a deterministic system. This increases the level of potential model detail over Markov-type models, improves parameter identifiability, and reduces computational burden and stochastic uncertainty compared with Monte Carlo methods. The example model of liver disease was able to accurately reproduce targets without extensive calibration or fitting and required minimal computational burden.

Conclusions

Explicitly modeling sojourn time distribution allows us to represent semi-Markov systems using detailed survival data from epidemiological studies without requiring sampling, avoiding the need for calibration, reducing computational time, and allowing for more robust probabilistic sensitivity analyses.

Highlights

Time-inhomogeneous semi-Markov models and other time-to-event–based modeling approaches can capture risks that evolve over time spent with a disease.

We describe an approach to computing these models that represents them as partial differential equations representing the evolution of the sojourn time probability density.

This sojourn time density framework incorporates complex data sources on competing risks and serial events while minimizing computational complexity.

Graphical Abstract

This is a visual representation of the abstract.

Keywords

cost-effectiveness analysis Markov models probabilistic sensitivity analysis simulation methods survival analysis semi-markov models time-to-event modelling

Health economic modeling evaluates the impact and cost-effectiveness of interventions by synthesizing and extrapolating data on health and cost outcomes. Almost all health economic models are state-transition models,¹ capturing a set of discrete health states an individual can be in and the likelihood of transitioning between these states. These models can simulate health care scenarios to predict outcomes, aiding informed decision making in resource allocation. There are various mathematical and computational approaches to state-transition modeling, each with their own strengths and limitations.²

Many health economic analyses use Markov models, which are easy to comprehend, implement, and compute, and they often have high parameter identifiability (i.e., transition rates can be supported by the available data) due to the clear delineation of health states. However, Markov models have notable disadvantages. The Markov property restricts the likelihood of a state transition to be either constant (time homogeneous) or dependent on only the current model time (time inhomogenous); they cannot depend on the time spent in the current state (the sojourn time). Consequently, shorter sojourn times for any particular state transition at a given time are more likely than any longer times.³ While this may be addressed with the addition of hidden states, this can lead to serious problems such as overfitting⁴ due to lack of identifiability,⁵ leading to poorly fitted models and spurious findings. A poor choice of time discretization/cycle length⁶ can also dramatically bias a Markov model’s outcomes, and altering the cycle length can lead to numerical errors.⁷

Semi-Markov models (or equivalently, Markov renewal processes)^8–10 can capture significantly more detail than Markov models can, allowing the time between state transitions to vary according to any distribution¹¹ rather than being constrained to a geometric distribution. This can dramatically increase model fidelity. This means the risk of developing a disease can depend directly on the length of time spent with a precursor condition,^12,13 the sojourn time. Detailed survival analyses from prior studies can be used directly to inform disease risk, which has been identified as a priority in health technology assessments.^14,15 This approach conserves identifiability by explicitly including states for which survival data are available. By further generalizing to time-inhomogeneous semi-Markov processes,¹¹ transitions can change over time, allowing the model to also capture risks that vary by age and secular trends. Semi-Markov models are also used in fields such as actuarial science^16–18 and engineering.¹⁹ For example, semi-Markov models of cancer control can be used to increase the fidelity of estimates of disease progression and the impact of interventions.^20,21

Unlike Markov models, the dynamics of semi-Markov models cannot usually be analytically computed. Instead, Monte Carlo sampling methods are typically used to estimate patient-level dynamics.^22–24 In this formulation, the semi-Markov model will describe a patient-level simulation approach²⁵ equivalent to a microsimulation or discrete event simulation,^26–29 which allows for a significant level of depth and detail.^2,30 However, Monte Carlo sampling introduces stochastic uncertainty to the model (first-order uncertainty).^25,31 In general, there is no consistent method of estimating the number of samples required to estimate a particular outcome. As a first-order approximation, if outcomes are assumed to be binomially distributed by the model, the length of confidence interval for the outcome is approximately $O (1 / \sqrt{N})$ , meaning an increasing number of simulations is required to accurately estimate an outcome, particularly for rare events. This can dramatically increase the computation time required for modeling, which both slows model development and limits the ability to conduct extensive sensitivity analyses, a priority for health economic estimates.³² The tradeoff between discrete-time Markov-type models and continuous-time semi-Markov simulation models can lead to many potential sources of bias and inaccuracy.³³

We developed the sojourn density model framework for explicit evaluation of semi-Markov and related models. This approach models the dynamics of the probability density of the sojourn time as a system of partial differential equations (PDEs), tracking both the likelihood of being in a given state and the time spent in that state, avoiding the need for stochastic sampling methods. The modeling framework is parametrized by the cause-specific hazard rates^34,35 (analogous to transition intensities³⁶ or the force of increment³⁷), allowing us to exploit survival data and methods³⁸ such as Kaplan–Meier estimators^34,39,40 and risk ratios.⁴¹ This captures the same detail as survival models, with the addition of serial state transitions. This approach is adapted from existing methods used in actuarial science,^11,16–18 which use the Kolmogorov forward equation.⁴² We have modified this for an epidemiology and health economic context, including adapting the notation and simplifying the construction where possible.

In this study, we describe the model structure for sojourn density models, compare their strengths and weaknesses versus other modeling approaches, and describe a numerical scheme for calculating these models. To demonstrate the advantages of this approach, we include an illustrative model of liver disease. In Australia, liver cancer is the seventh most common cause of cancer-related death,⁴³ with increasing incidence and mortality trends.⁴⁴ Routine surveillance of patients with liver cirrhosis (late-stage liver scarring) can increase detection at early stages and improve survival⁴⁵ but requires an understanding of the competing risk of liver cancer versus other-cause mortality. Our modeling approach allows us to incorporate competing risks of decompensation, cancer, and other-cause mortality and complex survival data,^46,47 which have been noted as key considerations in liver disease modeling.⁴⁸

Methods

We now describe the sojourn density model framework and demonstrate its useful mathematical properties. We also develop a numerical scheme for the calculation of these models.

Sojourn Density Model

Consider a state transition model with $N$ possible discrete states labeled $1, 2, . . ., N$ , such that at any time $t \in T \subset R$ (typically $T = [0, \infty)$ ), an individual is in exactly 1 of these states. We not only describe the likelihood of an individual being in a given state but also track the distribution of the length of time spent in that state, the sojourn time $τ \in T$ .

Doing so defines a continuous-time random process ${X_{t}}_{t \in T}$ with discrete state space ${1, 2, . . ., N}$ such that $P (X_{t} = i) = P (individual in state i at time t)$ . We can equivalently define the jump process ${Y_{n}, T_{n}}_{n = 0, 1, 2, \dots}$ such that $Y_{n} = X_{t} \forall t \in [T_{n}, T_{n + 1})$ , $Y_{n} \neq Y_{n + 1}$ .⁴⁹

The cause-specific hazard rates

λ_{i, j} (t, τ),

(1)

depending on both the model time $t$ and sojourn time $τ$ , are the key parameters for the model. These are the instantaneous transition rates for an individual in state $i$ to transition to state $j$ at time $t$ , given that the individual has so far spent time $τ$ in $i$ without any transitions (i.e, entered state $i$ at time $t - τ$ ). Hazard rates are widely used parameters in public health modeling.²² In some cases, $h$ is used for hazard rates.⁵⁰ They are equivalent to the force of increment discussed in equation 4.6 of Hoem¹¹ in the case $s = t$ (i.e., the reference time when the model is in state $i$ is the same as the start of the period over which the transition may occur).

As the transition likelihoods depend on only the current state and the sojourn time, the random process $X_{t}$ is a semi-Markov proccess,³⁶ as the sequence of transitions is Markov, but the time between transitions is not neccesarily exponentially distributed. By allowing the hazard rate to depend on $t$ , this generalizes to time-inhomogeneous semi-Markov models. This $t$ dependency can capture temporal trends as well as evolving risk with age $a = a_{0} + t$ , where $a_{0}$ is the age at $t = 0$ . While $X_{t}$ is a semi-Markov process, the jump process $(Y_{n}, T_{n})$ is a Markov renewal process.⁵¹

For the hazard rate $λ_{i, j} (t, τ)$ to represent the transition rate from state $i$ to state $j$ at time $t$ given sojourn time $τ$ , by equation (3) in Król et al.,²² the following should hold:

\begin{matrix} λ_{i, j} (t, τ) = lim_{Δ t \to 0} \frac{1}{Δ t} P (Y_{n + 1} = j, \\ T_{n + 1} \in [t, t + Δ t) | Y_{n} = i, T_{n} = t - τ) \\ = lim_{Δ t \to 0} \frac{1}{Δ t} P (X_{t + Δ t} = j | X_{s} = i \forall s \in [t - τ, t), X_{t - τ} \neq i) . \end{matrix}

(2)

For convenience, assume $λ_{i, i} (t, τ) = 0$ (i.e., there are no self-loops). If there is no possible transition between $i$ and $j$ , then $λ_{i, j} = 0$ .

To track the model state and the sojourn time, we introduce the dynamic density functions for the sojourn time:

f_{i} (t, τ) : T \times T \to [0, \infty)

(3)

for $i = 1, 2, \dots, N$ . For a fixed time $t \in T$ , these are probability density functions on the space of sojourn times and states $(τ, i) \in T \times {1, 2, . . . N}$ , describing the likelihood of being in a particular state with a particular sojourn at the present time. This notation is based on the notation introduced by Asanjarani et al.,³⁶ who described an equivalent density function with a single time/sojourn time variable.

Then we can define the process $X_{t}$ by the total probability mass in that state:

\begin{matrix} \begin{matrix} P (X_{t} = i) = \int_{0}^{\infty} f_{i} (t, s) ds = : g_{i} (t) . \end{matrix} \end{matrix}

(4)

We should expect the total probability mass to be 1, that is,

\sum_{i} g_{i} (t) = 1 .

(5)

Our goal now is to describe the dynamics of $f_{i} (t, τ)$ and then show this definition satisfies (2) and (5). Let $f_{i} (t, τ)$ satisfy the partial differential equation

\frac{\partial}{\partial t} f_{i} (t, τ) + \frac{\partial}{\partial τ} f_{i} (t, τ) = - (\sum_{j} λ_{i, j} (t, τ)) f_{i} (t, τ)

(6)

with boundary condition

f_{i} (t, 0) = \sum_{j} (\int_{0}^{\infty} λ_{j, i} (t, τ) f_{j} (t, τ) d τ) .

(7)

This equation is equivalent to the Kolmogorov forward equation.^17,42 The left-hand side of (6) represents the transport equation along lines of constant $t - τ$ , representing the likelihood of remaining in state $i$ while both $t$ and $τ$ increase, whereas the right-hand side captures the likelihood of making a transition to a subsequent state. Along the directional derivative $\nabla_{(1, 1)} f_{i}$ , the likelihood decays proportionally to the hazard rate (see the “Method of Characteristics and Survival Analysis” section). More rigorous derivations of this PDE are included elsewhere.^16,17 The findings here have been simplified for a health economic context and assume the hazard rates $λ_{i, j} (t, τ)$ are smooth.

The boundary condition (7) represents the accumulation of probability mass in state $i$ post-transition. By (5), the initial conditions $f_{i} (0, τ)$ should satisfy

\sum_{i} \int_{0}^{\infty} f_{i} (0, τ) d τ = 1 .

(8)

In addition,

lim_{τ \to \infty} f_{i} (t, τ) = 0

(9)

(i.e., supported on finite time).

Compare (6) to equation (2.37.1) in Helwich.¹⁷ Helwich considered the likelihood $p_{i j} (s, t, u, v)$ of transition from state $y$ to state $z$ such that $X_{s} = i$ with sojourn time $u$ and $X_{t} = j$ with sojourn time less than $v$ . In our notation, this is equivalent to fixing the model at state $i$ with sojourn time $u$ at time $s$ (i.e., $f_{k} (s, τ) = δ (τ - u)$ if $k = i$ else $0$ ), then tracking the evolution of $f_{j} (t, τ)$ for $t$ in a small neighbourhood of $s$ .

We now have a complete description of the random process $X_{t}$ through (4), (6), (7) and appropriate initial conditions. This constitutes what we will call a sojourn density model. With this current definition, sojourn density models are a type of semi-Markov model as a specific parametrization using hazard rates and computed through the sojourn density functions. In the sections “Sojourn-Preserving Transitions” and “Inhomogeneous Hazard Rates and Infectious Disease Modeling,” we introduce modifications that expand the sojourn density framework to include some non-semi-Markov models.

We now show that the process defined by (4) is well-defined and driven by the hazards $λ_{i, j} (t, τ)$ as expected.

Theorem 2.1. The random process defined by (4) with dynamics defined by (6) and (7) satisfies the conditions (2) and (5).

Proof. The boundary condition (7) can be justified by calculating the derivative of $g_{i} (t)$ , that is, the rate of change in likelihood of being in state $i$ at time $t$ :

\begin{matrix} \frac{d}{dt} g_{i} (t) = \frac{d}{dt} \int_{0}^{\infty} f_{i} (t, τ) d τ \\ = \int_{0}^{\infty} \frac{\partial}{\partial t} f_{i} (t, τ) d τ \\ = - \int_{0}^{\infty} \frac{\partial}{\partial τ} f_{i} (t, τ) d τ - \int_{0}^{\infty} (\sum_{j} λ_{i, j} (t, τ)) f_{i} (t, τ) d τ \\ = f_{i} (t, 0) - \int_{0}^{\infty} (\sum_{j} λ_{i, j} (t, τ)) f_{i} (t, τ) d τ \\ = \sum_{j} (\underset{entering state i}{\underset{︸}{\int_{0}^{\infty} λ_{j, i} (t, τ) f_{j} (t, τ) d τ}}) - \int_{0}^{\infty} (\underset{leaving state i}{\underset{︸}{\sum_{j} λ_{i, j} (t, τ)}}) f_{i} (t, τ) d τ . \end{matrix}

(10)

This satisfies the definition of the hazard rate: the likelihood of being in state $i$ decreases proportionally to the hazard of every transition leaving state $i \sum_{j} λ_{i, j} (t, τ)$ and increases for every transition entering state $i \sum_{j} λ_{j, i} (t, τ)$ proportional to the weight of those states $f_{j} (t, τ)$ . In addition,

\frac{d}{dt} (\sum_{i} g_{i} (t)) = 0,

(11)

that is, that the total probability mass is conserved. Therefore, with initial conditions satisfying equation (8) and boundary conditions (7), (5) is satisfied.

The propagation of hazards (2) through (6) can be observed through analogy with corollary 2.37 in Helwich¹⁷ as noted above; a proof in the current notation is included in Appendix A. ■

Method of Characteristics and Survival Analysis

We now compute the survival curves using method of characteristics. We refer to survival in the sense of the time until an event of interest, which may or may not be mortality, as is often used in epidemiological analyses.^34,52 The use of characteristics is equivalent to the approach in section 2.5 of Buchardt et al.¹⁶ or the explicit solution in equation 9 from Asanjarani et al.³⁶ As noted, equation (6) is best interpreted along characteristics of constant $t - τ$ corresponding to the state entry time (Figure 1). This allows explicit solutions through method of characteristics as well as analogies to other survival modeling approaches.

Figure 1

Characteristics of constant $t^{*} = t - τ$ for $f_{i} (t, τ)$ (left) and the corresponding elements $F_{i}^{m, n}$ with constant $m - n$ for the numerical scheme (right).

Consider the characteristic $t = t^{*} + τ$ for fixed $t^{*} \in R$ (Figure 1). Then, along this characteristic,

\frac{d}{d τ} f_{i} (t^{*} + τ, τ) = \frac{\partial}{\partial t} f_{i} (t^{*} + τ, τ) \frac{d}{d τ} (t^{*} + τ) + \frac{\partial}{\partial τ} f_{i} (t^{*} + τ, τ)

(12)

= - (\sum_{j} λ_{i, j} (t^{*} + τ, τ)) f_{i} (t^{*} + τ, τ) .

(13)

If $t^{*} > 0$ , then $t > τ$ , which corresponds to a characteristic starting from the boundary $τ = 0$ ; if $t^{*} \leq 0$ , then $t \leq τ$ , which corresponds to a characteristic from the initial condition $t = 0$ .

Defining $t_{0} = \max (t^{*}, 0)$ , $τ_{0} = \max (- t^{*}, 0)$ (so that $t^{*} = t_{0} - τ_{0}$ , $t_{0}, τ_{0} \geq 0$ ), the initial condition (at the observable start of the characteristic; not necessarily at $t = 0$ ) is at $(t_{0}, τ_{0})$ and so (12) has solution

f_{i} (t^{*} + τ, τ) = f_{i} (t_{0}, τ_{0}) e^{- λ_{i} (t^{*} + τ, τ)}

(14)

where

λ_{i} (t, τ) = \sum_{j} (\int_{τ_{0}}^{τ} λ_{i, j} (t - τ + s, s) ds) .

(15)

The function $λ_{i} (t, τ)$ is the overall cumulative hazard⁴¹ along the characteristic $t - τ = t^{*}$ . Then,

S_{i} (t, τ) = e^{- λ_{i} (t, τ)},

(16)

is the overall survival (i.e., the likelihood no transitions occur) along this characteristic, and $f_{i} (t^{*} + τ, τ) = f_{i} (t_{0}, τ_{0}) S_{i} (t^{*} + τ, τ)$ .

We can correspondingly calculate the time-to-event distribution after entering state $i$ at time $t - τ$ as the distribution

λ_{i, j} (t, τ) S_{i} (t, τ)

(17)

across $j = 1, 2, \dots, N$ . This allows translation to time-to-event epidemiological models and access to time-to-event methodology.^35,53,54

It follows that the probability of the next transition being from $i \to j$ before time $t = t^{*} + τ$ given entry to state $i$ at time $t_{0}$ , is

\begin{matrix} \begin{matrix} P (next transition time < t \cap next transition i \to j | \\ entered state i at time t_{0}) \\ = P (Y_{n + 1} = j, T_{n + 1} < t | Y_{n} = i, T_{n} = t_{0}) \\ = \int_{t_{0}}^{t} λ_{i, j} (s, s - t^{*}) S_{i} (s, s - t^{*}) ds . \end{matrix} \end{matrix}

(18)

The limit of this is the likelihood of the next transition being to a particular state:

\begin{matrix} P (next transition i \to j | entered state i at time t_{0}) \\ = \int_{t_{0}}^{\infty} λ_{i, j} (s, s - t^{*}) S_{i} (s, s - t^{*}) ds . \end{matrix}

(19)

We can also calculate the cause-specific survival, assuming other transitions (or events) are censored:

S_{i, j} (t, τ) = e^{- \int_{τ_{0}}^{τ} λ_{i, j} (t - τ + s, s) ds} .

(20)

Where cause-specific survival data (equivalently, cause-specific cumulative incidence) is available with competing events censored, this may provide a calibration target for the hazard rates, which can be directly identified from the survival model (see “Hazard and Cumulative “Hazard in Clark et al.³⁴); parametrization of $λ_{i, j} (t, τ)$ across $t^{*}$ can then be identified.

Compare this with the common approach in analytical epidemiological models in which the likelihood of transition is computed explicitly as the convolution of the $τ$ -dependent hazard rate and the time of entry.^55–57 In our formulation, this is reflected in the boundary conditions

\begin{matrix} \begin{matrix} f_{i} (t, 0) = \sum_{j} (\int_{0}^{\infty} λ_{j, i} (t, τ) f_{j} (t, τ) d τ) \\ = \sum_{j} (\int_{0}^{\infty} f_{j} (t - τ, 0) λ_{j, i} (t, τ) S_{j} (t, τ) d τ) . \end{matrix} \end{matrix}

(21)

If we define entry time distribution $F_{j} (x) = f_{j} (x, 0)$ and survival-adjusted cause-specific hazard rate $H_{(j, i)} (x) = λ_{j, i} (t, x) S_{j} (t, x)$ , then $F_{i} (t) = \sum_{j} (F_{j} * H_{j, i}) (t)$ , the convolution of these. Our approach generalizes this to incorporate serial transitions.

Initial Conditions

The appropriate initial conditions are highly dependent on the system being modeled. The Dirac delta function around $τ = 0$ is a potential initial condition, corresponding to starting in state $i$ with no existing sojourn time; explicitly,

f_{j} (0, τ) = {\begin{matrix} δ (τ) & if i = j \\ 0 & if i \neq j . \end{matrix}

(22)

In practice, modeling may start in situ, that is, after some sojourn time, informing a distribution for the initial condition. For example, this could reflect the time an individual would have spent in a health state by the start of a trial (“Incorporating Lead-Time Biases” section).

Numerical Scheme

The characteristic solutions developed in the “Method of Characteristics and Survival Analysis” section can be used to develop a numerical scheme. For a discretization of $T$ by fixed time steps $Δ t \in R^{+}$ , define

\begin{matrix} \begin{matrix} F_{i}^{m, n} : = f_{i} (m Δ t, n Δ t) \\ = {\begin{matrix} f_{i} ((m - n) Δ t, 0) e^{- λ_{i} (m Δ t, n Δ t)} for m \geq n \\ f_{i} (0, (n - m) Δ t) e^{- λ_{i} (m Δ t, n Δ t)} for m < n \end{matrix} \\ = {\begin{matrix} F_{i}^{m - n, 0} e^{- λ_{i} (m Δ t, n Δ t)} for m \geq n \\ F_{i}^{0, n - m} e^{- λ_{i} (m Δ t, n Δ t)} for m < n \end{matrix} \end{matrix} \end{matrix}

(23)

for $m, n \in ℤ^{+}$ . Each value of $m - n \in Z$ corresponds to a characteristic with $t^{*} = (m - n) Δ t$ . For $m < n$ , the initial condition is $F_{i}^{0, n - m} = f_{i} (0, (n - m) Δ t)$ . For $m \geq n$ , the initial condition $F_{i}^{m - n, 0} = f_{i} ((m - n) Δ t, 0)$ is the boundary condition of the full PDE.

The cumulative hazard functions $λ_{i} (t, τ)$ may be computed by (15) in simple cases or otherwise approximated. Equivalently, each $F_{i}^{m, n}$ element can be calculated iteratively by

F_{i}^{m, n} = F^{m - 1, n - 1} e^{- \sum_{j} \int_{0}^{Δ t} λ_{i, j} ((m - 1) Δ t + s, (n - 1) Δ t + s) ds} .

(24)

(see Figure 1 for an illustration of this; compare Buchardt et al.¹⁶). In many cases, the integral in (24) can be explicitly computed; otherwise, it can be approximated by the midpoint $\int_{0}^{δ t} λ_{i, j} ((m - 1) Δ t + s, (n - 1) Δ t + s) d s \approx Δ t λ_{i, j} ((m - 0.5) Δ t, (n - 0.5) Δ t)$ . The exponential function could similarly be approximated using the appropriate Taylor polynomial for improved performance.

A first-order scheme for the boundary conditions can be computed as

\begin{matrix} \begin{matrix} F_{i}^{m, 0} = f_{i} (m Δ t, 0) \\ = \sum_{j} (\int_{0}^{\infty} λ_{j, i} (m Δ t, τ) f_{j} (m Δ t, τ) d τ) \\ \approx \sum_{j} (\sum_{k > 0} Δ t λ_{j, i} (m Δ t, k Δ t) f_{j} (m Δ t, k Δ t)) \\ = \sum_{k > 0, j} (Δ t λ_{j, i} (m Δ t, k Δ t) F_{j}^{m, k}) \end{matrix} \end{matrix}

(25)

for $m \geq 0$ .

One can subsequently estimate the evolution of the random process $X_{t}$ as

\begin{matrix} \begin{matrix} P (X_{t} = i) = g_{i} (m Δ t) \\ \approx \sum_{k = 0}^{\infty} Δ t f_{i} (m Δ t, k Δ t) \\ \approx \sum_{k = 0}^{\infty} Δ {tF}_{i}^{m, k} . \end{matrix} \end{matrix}

(26)

An improvement can be made to this numerical scheme by adjusting the boundary conditions to exactly preserve the total mass at all timesteps:

F_{i}^{m, 0} = \sum_{j, k} {\tilde{h}}_{j, i}^{m, k} (F_{j}^{m - 1, k - 1} - F_{j}^{m, k})

(27)

where

{\tilde{h}}_{i, j}^{m, k} = \frac{λ_{i, j} ((m - \frac{1}{2}) Δ t, (k - \frac{1}{2}) Δ t)}{\sum_{l} λ_{i, l} ((m - \frac{1}{2}) Δ t, (k - \frac{1}{2}) Δ t)} .

(28)

We can then verify that $\sum_{i, k} F_{i}^{m, k}$ is constant; that is, the total probability mass $\sum_{i} g_{i} (t)$ is conserved. We have here used the midpoint estimates for $λ_{i, j}$ . This scheme can improve accuracy when there is a significant scale difference between hazards for competing events.

In epidemiological contexts, hazards are often modeled as continuous and smooth functions, so this scheme will have high accuracy and be convergent due to the inclusion of attractor death states.

A straightforward way to approach the numerical scheme is to calculate $F_{i}^{m, n}$ across all $i$ and all $m < 0$ , tracking the accumulated boundary conditions $F^{m, 0}$ for $m > 0$ , and then computing $F_{i}^{m, n}$ for each $m \geq 0$ across all $i, n$ . This has been implemented in Supplementary Material: Example Code.

Hazard Functions and Common Parametric Distributions

The sojourn time density modeling scheme is determined by the hazard rates $λ_{i, j} (t, τ)$ . It is common to describe these rates parameterically^25,58 to define a manageable set of model parameters. For example, constant hazard rates $λ_{i, j} (t, τ) = c$ correspond to an exponential time-to-event distribution, while $λ_{i, j} (t, τ) = bk τ^{k - 1}$ and $λ_{i, j} (t, τ) = p e^{r τ}$ correspond to Weibull and Gompertz distributions, respectively.⁵⁹ There are many other potential choices for $λ_{i, j} (t, τ)$ , which may vary with $τ$ , $t$ , or both.

The choice of hazard function is highly dependent on the nature of the system being modeled and the available data sources.⁶⁰ van Wijk et al.⁵⁸ provided a tutorial and general overview of determining appropriate parametrized hazard functions, including periodic functions that can represent seasonal variations or other cyclical patterns. Hazard rates can be estimated from processed survival data, such as smoothed⁶¹ Kaplan–Meier or Fine–Grey estimates or time-to-event distributions.^36,62 Covariates can be directly included in hazard rates; for instance, Cox proportional hazard ratios can be used to define ${\tilde{λ}}_{i, j} (t, τ | x) = λ_{i, j} (t, τ) e^{β \cdot x}$ for some vector of covariates $x \in R^{n}$ and a baseline hazard function $λ_{i, j} (t, τ)$ .^41,63 Constant multipliers for the hazards can be used as calibration targets, to maintain the shape of the distribution.⁶⁴

Asanjarani et al.³⁶ showed how time-to-event/discrete event simulation distributions^38,65 can be translated into hazard rates. Hazard rate can often be defined piecewise,⁶⁴ for example, from life tables.^66,67 There are also more advanced methods for calculating and representing hazard rates, including nonparametric smoothing⁶¹ and rational polynomial fraction parametrization.⁶⁸ Arbitrarily detailed and complex hazard rates can be used in sojourn density models; an appropriate level of detail supported by the data must be chosen to avoid overfitting.

Sojourn-Preserving Transitions

In an epidemiological model, there are often transitions that are independent of the underlying growth of disease. For instance, in a model of cancer development and screening, entry into an “surveillance” cohort would not affect the progress of the disease. This can be reflected in our model framework by introducing transitions that preserve the sojourn time.

To achieve this, we can introduce a new set of hazard rates $μ_{i, j} (t, τ)$ that are the transition rates for an individual in state $i$ to transition to state $j$ , at time $t$ , given that the individual has already spent time $τ$ in the state without any transitions occurring, while conserving the sojourn time. These rates $μ_{i, j}$ will depend on the type of transition being modeled but would typically have form similar to the $λ_{i, j}$ hazard rates (“Hazard Functions and Common Parametric Distributions” section). Then, equation (6) can be modified as follows:

\begin{matrix} \begin{matrix} \frac{\partial}{\partial t} f_{i} (t, τ) + \frac{\partial}{\partial τ} f_{i} (t, τ) = - (\sum_{j} λ_{i, j} (t, τ)) f_{i} (t, τ) \\ - (\sum_{j} μ_{i, j} (t, τ)) f_{i} (t, τ) + (\sum_{k} μ_{k, i} (t, τ) f_{k} (t, τ)) \end{matrix} \end{matrix}

(29)

with boundary condition (7). With this formulation, theorem 2.1 holds with minor alterations.

This approach may be useful for models in which there are multiple categories of states; for example, a health economic model may have a “health outcome” state space ${1, 2, 3, \dots}$ (e.g., no lesions, precancerous lesions, cancer) as well as an “intervention” state space ${a, b, c . . .}$ (e.g., screened, unscreened). The state space corresponding to the Cartesian product of these $(1, a), (1, b), (2, a)$ and so forth can be used for a model that tracks both, with transitions between intervention states governed by the $λ$ parameters and transitions between the intervention states governend by $μ$ .

This can be further extended by continuing to track the sojourn time across multiple independent types of transitions. For instance, one could classify $k$ independent “types” of transitions $λ_{i, j}^{1}$ , $λ_{i, j}^{2}$ , … $λ_{i, j}^{k}$ corresponding to $k$ independent sojourn times $τ^{1}, τ^{2}, . . ., τ^{k}$ which are each tracked by the model. Then (29) can be generalized to keep track of multiple sojourn times as

\begin{matrix} \begin{matrix} \frac{\partial}{\partial t} f_{i} (t, τ) + \sum_{k} \frac{\partial}{\partial τ^{k}} f_{i, j} (t, τ) = - (\sum_{k, j} λ_{i, j}^{k} (t, τ)) f_{i} (t, τ) \end{matrix} \end{matrix}

(30)

where $τ = (τ^{1}, τ^{2}, \dots, τ^{k})$ , with new boundary conditions for each sojourn time $τ^{i}$

f_{i} (t, τ^{1} (1 - δ_{1, n}), τ^{2} (1 - δ_{2, n}), . . ., τ^{k} (1 - δ_{k, n})) = \sum_{k, j} (\int_{0}^{\infty} λ_{i, j}^{k} (t, τ) f_{j} (t, τ) d τ^{k})

(31)

for $n = 1, 2, \dots k$ where $δ_{x, y}$ is the Kronecker delta function. Then, each transition type $λ^{k}$ “sends” the density to the corresponding boundary condition at $τ^{k} = 0$ . This expansion allows for a more general and class of models.

Computing Health Economic Outcomes

In addition to the likelihood of being in each state, there are other outcomes relevant for health economic analyses that can be computed from these models. These can usually be assessed via either the expected time spent in each state or the likelihood of transition into a state.

For example, the quality-adjusted life-years (QALYs⁶⁹) can be calculated based on the expected time spent in each state. If a person in state $i$ has a QALY value of $x_{i}$ , the expected QALY value over the period $t_{0}$ to $t_{1}$ is the weighted sum

\sum_{i} x_{i} \int_{t_{0}}^{t_{1}} g_{i} (t) dt .

(32)

Outcomes such as costs associated with entering a particular state can be calculated via the boundary condition. For instance, the likelihood of entering state $i$ over the period $t_{0}$ to $t_{1}$ is

\int_{t_{0}}^{t_{1}} f_{i} (t, 0) dt .

(33)

If there is a cost $α_{i}$ associated with entering state $i$ , the total expected cost across this period is

\sum_{i} α_{i} \int_{t_{0}}^{t_{1}} f_{i} (t, 0) dt .

(34)

If there are ongoing costs associated with the time spent in a particular state, these can be calculated in the same way as QALYs above.

Inhomogeneous Hazard Rates and Infectious Disease Modeling

By relaxing the assumptions on the form of the hazard functions, the sojourn time density models can be used to simulate a more general class of models than semi-Markov models. For example, we could allow the hazard to depend on the current model state, that is, have the form

λ_{i, j} (t, τ, g (t))

(35)

where $g (t) = {g_{i} (t)}$ is a vector of the values (4). One can then simulate models where the effect size is dependent on the current state of the system.

For example, the standard SIR model⁷⁰ has states $S$ (susceptible), $I$ (infectious), $R$ (recovered). The risk of infection is proportional to the currently infected population, leading to the hazard

\begin{matrix} \begin{matrix} λ_{S, I} (t, τ, g (t)) = β P (X_{t} = I) \\ = β \int_{0}^{\infty} f_{I} (t, τ) d τ \\ = β g_{I} (t) . \end{matrix} \end{matrix}

(36)

The recovery rate for infected individuals is constant, that is,

λ_{I, R} (t, τ, g (t)) = γ .

(37)

These hazards can be used to define a sojourn density model; we can confirm that

\frac{d g_{S}}{dt} = - β g_{S} (t) g_{I} (t), \frac{d g_{I}}{dt} = β g_{S} (t) g_{I} (t) - γ g_{I} (t), \frac{d g_{R}}{dt} = γ g_{I} (t),

(38)

that is, the defining equations for the standard SIR model. Using this starting point, we could modify this to introduce more complex risks of infection and recovery. For instance, a distribution for time to disease recovery could be simulated by modifying $λ_{I, R}$ to be a function of $τ$ , the time in infectious state, to reflect an inflection that can clear after only a few days minimum. In the sojourn density model framework, new dynamics can be included by modifying the hazard functions rather than by adjusting the design of the model.

This addition of heterogeneous generalizes the sojourn density model framework and could allow for the synthesis of different model types into one framework. For example, a model of both a viral infection and a resulting disease (such as hepatitis and liver cancer) could be directly modeled. The addition of proportional terms to the hazard does not affect theorem 2.1 but does mean the method of characteristics described in the Method of Characteristics and survival analysis section may no longer be applied.

Incorporating Lead-Time Biases

As the sojourn time is modeled explicitly in this framework, this allows us to incorporate lead-time bias. Lead-time biases occur in many contexts; for instance, screen detection of cancer artificially increases survival times due to the increased observation window.^71,72 A related bias occurs in modeling that has been informed by survival data measured from the start of a trial rather than onset of disease.

These biases can be addressed by modifying the hazard rate. A highly simplistic approach may be to define the adjusted hazard rate

{\tilde{h}}_{i, j} (t, τ) = {\begin{matrix} 0 & for τ < τ_{0} \\ λ_{i, j} (t, τ - τ_{0}) & τ \geq τ_{0} \end{matrix}

(39)

where $h_{i, j} (t, τ)$ is the observed cause-specific hazard rate from a trial that may include a lead-time bias and $τ_{0}$ is a cutoff time. This assumes that there is no risk of transition until some given time $τ_{0}$ . Alternatively, the hazard could be assumed to increase to the observed rates at time $τ_{0}$ , linearly or otherwise. In practice, a useful estimate for $τ_{0}$ may be difficult to determine, requiring a sensitivity analysis over a range; there are many useful methods for estimating the lead time.^73–76

Example Model: Liver Disease

To demonstrate our new approach to modeling, we now describe a model of liver disease and development of cancer. A key intervention in the prevention of liver cancer death is early detection via routine surveillance. Recommendations like 6-monthly ultrasound for people with cirrhosis (late-stage liver scarring) can increase the likelihood of detection at early stages when curative treatment is possible.⁷⁷ However, people with cirrhosis are also at high risk of both liver decompensation (irreversably worsened liver damage) and noncancer death. People who develop liver decompensation or die before the onset of liver cancer cannot benefit from early detection. It is therefore critical to understand the competing risks of liver cancer, liver decompensation, and other-cause death when assessing whether routine surveillance is advisable or whether the costs and patient discomfort caused by surveillance outweigh the health benefits. This model shows how the risk of liver cancer, and competing risk of decompensation, can be captured in a way that evolves with the length of time spent with liver cirrhosis.

The structure of the simplified model is shown in Figure 2 and captures individuals with liver cirrhosis, decompensation events, onset of liver cancer, and death. As the model is designed to illustrate the potential utility of surveillance for liver cancer in the presence of competing risks of decompensation and death, liver cancer and death were used as terminal states. Hazard rates for relevant transitions (decompensation, onset of liver cancer, death) were calibrated by identifying appropriate parametric forms and calibrating to relevant data sources. With this approach, the design of the model states and flow is kept simple, while the complexity is shifted to computing the hazard rates determining the transitions between states. Each cause-specific hazard was fitted in turn, and the overall fit assessed. For this iteration of the model, a cohort aged 55 y was simulated, based on the mean age of cirrhosis diagnosis.

Figure 2

Schematic of the example liver disease model and the sojourn- and time-dependent hazards for competing risks for people with compensated cirrhosis.

The modeling described in this article is part of a larger model developed to inform the “Clinical Practice Guidelines for Liver Cancer Surveillance for People at High Risk in Australia.”⁷⁸ The development of these guidelines, including the modeling, was led by clinicians in the area of liver disease and liver cancer alongside a multidisciplinary working party, which included health care and clinical representatives, representatives with lived experience, and other community representatives. The full model includes detailed liver cancer stage progression and survival as well as economic details. For further details, see Worthington et al.⁷⁹

For this analysis, published survival curves were used as target data for the modeling. Where patient-level data are available, this could be used as a calibration target instead, using methods in the existing literature on semi-Markov models.²²

Results

Performance of Sojourn Density Models

The described model structure can capture the same level of detail as semi-Markov models and many other patient-level simulations, increasing the level of detail while allowing for deterministic evaluation. The computational complexity of the numerical algorithm described in the “Numerical Scheme” section is $O (T^{2} S^{2})$ , where $T$ is the number of timesteps and $S$ is the number of model states. Vectorization across characteristics can improve performance considerably, for example, via the Python NumPy library or other scientific computing libraries. For models in which not every state has a transition to every other state, the complexity is $O (T^{2} U)$ where $U < S^{2}$ is the number of possible transitions. Compare this with a patient-level simulation, including most implementations of semi-Markov models, in which the computational complexity is driven by the number of patients simulated $O (TN)$ , where $N$ is the number of simulated patients. As noted previously, the precision of such models is approximately $O (1 / \sqrt{N})$ , and so in some applications, $N$ may need to be very large for precise estimates. This is avoided in the sojourn time density model. For comparison, a standard numerical algorithm for a Markov model would have computational complexity $O (T S^{2})$ .

In practice, sojourn density models can run in less than 1 s on a standard computer. This allows for runs with a variety of parameter values to assess sensitivity analyses, particular with parallel computing. In most cases, the optimal timestep $Δ t$ can be identified by assessing a particular parameter set, selecting a small enough value for the results of interest to converge. This choice will be highly dependent on the problem being assessed, roughly scaling with the expected transition time. Future work is planned to assess the effect of the timestep.

With an appropriate implementation of the numerical algorithm described in the “Numerical Scheme” section, computing a sojourn density model requires specifying only the model states and the form of the hazard rates between them. Listing 1 shows the syntax of the code for this model; the full code, including the implementation of the numerical scheme, has been included in Supplementary Material: Example Code.

Listing 1

Excerpt Showing Code Syntax for the Liver Disease Model^a

Liver_Disease_Model=Sojourn_Model.Model()
state_names=[“Compensated cirrhosis”,“Decompensated cirrhosis”,“Liver
↪ cancer”,“Death”]
Liver_Disease_Model.add_state(state_names)
#Liver cancer onset transition
Liver_Cancer_Hazard=Hazard_Functions.Hazard( Hazard_Type =
↪ “gompertz_sojourn”, Parameters={“p”:1.197e-02,“r”:6.654e-02})
Liver_Disease_Model.add_transition(“Compensated cirrhosis”,“Liver
↪ cancer”,Liver_Cancer_Hazard)
///Other transitions omitted for brevity
res_df=Liver_Disease_Model.run_model(dt=0.01, max_t=10,
↪ init_conds=[“Compensated cirrhosis”])

Full code is included in Supplementary Material: Example Code.

Liver Disease Model Calibration

The model was calibrated to reproduce key targets, with model parameters shown in Table 1. For each target, an appropriate parameterized hazard rate was selected, and optimal parameters were calculated using the BFGS algorithm in the SciPy package. We also calibrated distributions for hazard ratios for each of the hazard rate functions. These were calibrated to reproduce the 95% confidence intervals in each data source across the range of hazard ratio values and act as multipliers for each distribution for the probabilistic sensitivity analysis. The model was well-fitted to the target data, reproducing terminal values for the survival curves and remaining within the 95% confidence intervals for the duration. Key targets are shown in Figure 3. For all-cause death rates, a piecewise constant hazard was used based on age-specific rates published by the Australian Bureau of Statistics.⁸²

Table 1

Calibrated Parameters for Example Model of Liver Disease

Parameter	Model Value (95% PI)	Target (95% CI)	Note
Cirrhosis decompensation hazard	$λ = 0.0129 t + 0.00668 for t < 6, 0.0359 for t \geq 6$		Calibrated to Vilar-Gomez et al.⁸⁰
Hazard ratio distribution for PSA	$Lognormal (μ = - 0.0213, σ^{2} = {0.117}^{2})$		Calibrated to Vilar-Gomez et al.⁸⁰
Annual decompensation risk (target)	1.03% (0.9%–1.2%)	1.59% (0.17%–3.01%)	Calibration target⁸⁰
Ten-year decompensation risk	34.5% (29.5%–40.2%)	33.7% (28.4%–39.1%)
Liver cancer onset hazard	$λ = 0.0102 τ^{1.377}$		Calibrated to Vilar-Gomez et al.⁸⁰
Hazard ratio distribution for PSA	$Lognormal (μ = - 0.00171, σ^{2} = {0.158}^{2})$		Calibrated to Vilar-Gomez et al.⁸⁰
Annual liver cancer risk	0.99% (0.8%–1.3%)	1.53% (0.14%–2.93%)	Calibration target⁸⁰
Ten-year liver cancer risk	15.2% (11.8%–19.2%)	15.6% (11.5%–19.7%)
Death hazard (compensated cirrhosis)	$λ = 0.0336 τ^{- 0.238}$		Calibrated to D’Amico et al.⁸¹
Hazard ratio distribution for PSA	$Lognormal (μ = - 0.0254, σ^{2} = {0.147}^{2})$		Calibrated to D’Amico et al.⁸¹
Annual death risk (compensated)	4.72% (3.8%–5.8%)	3.97% (2.62%–5.31%)	Calibration target⁸¹
Ten-year death risk (compensated)	39.4% (35.2%–43.6%)*	39.4% (36.1%–42.8%)
Death hazard (decompensated cirrhosis)	$λ = 0.447 - 0.150 τ + 0.0237 τ^{2}$		Calibrated to D’Amico et al.⁸¹
Hazard ratio distribution for PSA	$Lognormal (μ = 0.0478, σ^{2} = {0.0895}^{2})$		Calibrated to D’Amico et al.⁸¹
Annual death risk (decompensated)	37.8% (33.1%–42.4%)*	38.3% (35.0%–41.6%)	Calibration target⁸¹
Ten-year death risk (decompensated)	90.6% (89.0%–91.9%)*	91.1% (89.2%–93.0%)
Death hazard (other causes)	Age-specific annual rates		Australian life tables⁸²

CI, confidence interval; PSA, probabilistic sensitivity analysis.

95% PI: 95% prediction interval of model outcomes, based on samples from the specified hazard ratio distributions.

Figure 3

Selected calibration outcomes. Top row: Liver cancer incidence rates (left; target Vilar Gomez et al.⁸⁰) and decompensation rates (right). Bottom row: all-cause death rates in patients with compensated (left; includes any decompensation events) and decompensated (right) cirrhosis.

A code snippet showing the syntax to define a model using the included code is shown in Listing 1. Examples of the kinds of outputs that can be generated are shown in Figure 4. This includes the likelihood of $X_{t}$ being in a given state by time (through $g (t)$ ). This was used to calculate the 5-y risk of liver cancer or non–liver cancer death in people who have not yet developed liver cancer or died, a key outcome for assessing the benefits of liver surveillance.

Figure 4

Outcomes from the example model of liver disease. Left: evolving likelihood of liver decompensation, liver cancer, and non–liver cancer death. Right: 5-y risk of liver cancer versus 5-y risk of non–liver cancer death for patients who have not yet developed liver cancer/died, including 95% model range based on probabilistic sensitivity analysis.

A probabilistic sensitivity analysis was run by sampling hazard ratios from the distributions given in Table 1. The sensitivity analysis ran 100,000 parameter sets, which was completed in under 2 h. In the 5-y risk outcomes in Figure 4, this was used to develop a prediction interval; 95% of the simulated parameter sets led to outcomes in the indicated range.

Discussion

We described a novel approach to epidemiological modeling by calculating the evolution of the probability density of the sojourn time. This allows detailed survival data sources to be used to model competing risks. We demonstrated this sojourn density model approach by presenting a model of liver disease. The model was calibrated and reproduced the targets accurately with low computational burden.

The novel sojourn density model structure allows for flexible analyses that rely on detailed survival data and flexible time scales. This gives advantage over typical epidemiological modeling approaches, improving detail without sacrificing interprability while limiting computational burden. Designing models that reduce computational burden to facilitate probabilistic sensitivity analyses has been identified as a priority for health technology assessments,³² often requiring surrogate methods such as meta-modeling/emulators to reduce the reliance on computationally expensive models.^83–86

By using the hazard rates as parameters, the models can be interpreted directly in the context of source epidemiology studies, and parameter values can be directly inferred from survival data. In the context of liver disease modeling, this means the model can reflect the true health benefits of early-stage diagnosis by accurately capturing disease progression. The ability to numerically evaluate our model at any time scale also allows us to accurately model phenomena such as surveillance intervals, which were not been included in real-world trials with a high degree of precision. Avoiding a fixed time discretization reduces the likelihood of numerical errors if chosen poorly⁸⁷ and improves the ability to capture short time scale events.

There is great potential to refine survival analysis estimates designed to guide policy recommendations,¹⁵ and it is hoped that this modeling approach allows for the development of straightforward models that incorporate existing data sources in a clear and flexible way, reducing the need for model calibration and minimizing data postprocessing and model design effort. The benefits of this approach are demonstrated in the closeness of fit to the calibration targets (Figure 3) as well as the flexibility to evaluate surveillance at a wide range of intervals.⁷⁹ By providing a numerical implementation, we have also reduced the implementation burden on health economic practitioners, requiring only the specification of the model states and the hazard rates associated with the transitions between them.

The limitations of this model structure include the semi-Markov property; each state transition depends only on the current state, the sojourn time, and the model time. This approach may therefore not be appropriate in contexts where multiple prior events are required to inform future likelihoods.

The liver disease model described here is a part of a larger model of liver disease,⁷⁹ which was used to evaluate surveillance recommendations for the “Clinical Practice Guidelines for Hepatocellular Carcinoma Surveillance for People at High Risk in Australia.”⁷⁸ The flexible sojourn density model allowed us to analyze combinations of surveillance technologies and intervals to analyze more complex surveillance algorithms and optimize these algorithms through iterative design enabled by the fast computation time and high precision.

There are many potential future developments to build on the work presented here. Higher-order approximations for the integrals in (25) and (24) could be used to derive more accurate numerical schemes, which may help when dealing with pathological hazards such as those that are nonsmooth or have discontinuities. It may also be possible to develop spectral methods to approximate solutions without time discretization.⁸⁸ Further work could also use these models to assess the potential impact of lead-time biases, and the loss in precision when these are not considered in modeling. It is hoped that the method described here will be of general use to health economic practitioners.

Conclusions

To inform health policy and make an impact, epidemiological models must be able to provide reliable, timely, and interpretable results. Sojourn density models can incorporate published survival data and generate competing risk estimates quickly, allowing for simple design and fast calibration. We hope this approach will be of use to modelers in public health.

Supplemental Material

sj-pdf-1-mdm-10.1177_0272989X251333398 – Supplemental material for Evaluating Semi-Markov Processes and Other Epidemiological Time-to-Event Models by Computing Disease Sojourn Density as Partial Differential Equations

Supplemental material, sj-pdf-1-mdm-10.1177_0272989X251333398 for Evaluating Semi-Markov Processes and Other Epidemiological Time-to-Event Models by Computing Disease Sojourn Density as Partial Differential Equations by Joachim Worthington, Eleonora Feletto, Emily He, Stephen Wade, Barbara de Graaff, Anh Le Tuan Nguyen, Jacob George, Karen Canfell and Michael Caruana in Medical Decision Making

Footnotes

Acknowledgements

The authors would like to acknowledge the Expert Advisory Group for the Optimising Liver Cancer Control in Australia project convened by Cancer Council Australia. We would also like to thank Cancer Council Australia for coordination of this project.

Materials reported in this article, including the modeling and results, were developed to support the Optimising Liver Cancer Control in Australia project and were originally developed and reported in “Clinical Practice Guidelines for Hepatocellular Carcinoma Surveillance for People at High Risk in Australia: Appendix D7,” a non–peer-reviewed report.

The authors would like to acknowledge the working group for the “Clinical Practice Guidelines for Hepatocellular Carcinoma Surveillance for People at High Risk in Australia.” A full list of working group members is available in “Clinical Practice Guidelines for Hepatocellular Carcinoma Surveillance for People at High Risk in Australia: Appendix I.”

The authors would like to thank the Daffodil Centre Systematic Review team, who provided assistance in reviewing the literature on existing models of hepatocellular carcinoma and surveillance.

The authors would also like to thank the reviewers and editors at Medical Decision Making for their handling of this article and the substantial contributions and improvements they collectively made to the direction and style of the final article.

This article uses custom data sets from the Australian Institute of Health and Welfare and the Cancer Institute NSW. We thank the Australian Institute of Health and Welfare and the population-based cancer registries of New South Wales, Victoria, Queensland, Western Australia, South Australia, Tasmania, the Australian Capital Territory, and the Northern Territory for the provision of data from the Australian Cancer Database. We would also like to thank the Cancer Institute NSW for their assistance with the data request.

Karen Canfell receives salary funding from the National Health and Medical Research Council Australia (NHMRC Leadership Fellowship APP1194679).

The authors declared the following potential conflicts of interest with respect to the research, authorship, and/or publication of this article: Karen Canfell is co-principal investigator (PI) of an investigator-initiated trial of human papillomavirus (HPV) screening in Australia (Compass), which is conducted and funded by the Australian Centre for the Prevention of Cervical Cancer (ACPCC), a government-funded health promotion charity. The ACPCC has previously received equipment and a funding contribution for the Compass trial from Roche Molecular Systems USA. She is also co-PI on a major implementation program, “Elimination of Cervical Cancer in the Western Pacific,” which receives support from the Minderoo Foundation and equipment donations from Cepheid Inc. Michael Caruana is an investigator on an investigator-initiated trial of cytology and primary HPV screening in Australia (“Compass”; ACTRN12613001207707 and NCT02328872), which is conducted and funded by the Australian Centre for the Prevention of Cervical Cancer, a government-funded health promotion charity. The Australian Centre for the Prevention of Cervical Cancer has received equipment and a funding contribution for the Compass trial from Roche Molecular Systems and operational support from the Australian government. However, neither MC nor his institution on his behalf (the Daffodil Centre, a joint venture between Cancer Council NSW and The University of Sydney) receive direct or indirect funding from industry for Compass Australia. All other authors declare no potential conflicts of interest with respect to the research, authorship, and/or publication of this article. The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Cancer Council Australia financially supported this study as part of the Optimising Liver Cancer Control in Australia project, which was funded through the Australian Government Department of Health and Aged Care. The funding agreement ensured the authors’ independence in designing the study, interpreting the data, writing, and publishing the report.

Author’s Note

Anh Le Tuan Nguyen is also affiliated to Hanoi Medical University, Vietnam

Ethical Considerations

Not applicable. This article does not report on animal or human studies and did not use any individual-level data to inform the analysis.

Consent to Participate

Not applicable. This article does not report on animal or human studies and did not use any individual-level data to inform the analysis.

Consent for Publication

Not applicable. This article does not report on animal or human studies and did not use any individual-level data to inform the analysis.

ORCID iDs

Joachim Worthington

Eleonora Feletto

Stephen Wade

Anh Le Tuan Nguyen

Data Availability

The code used to generate these figures has been provided as a supplementary file at , including all relevant parameters.

References

Siebert

Alagoz

Bayoumi

, et al. State-transition modeling: a report of the ISPOR-SMDM modeling good research practices task force-3. Value Health. 2012;15(6):812–20. DOI: 10.1016/j.jval.2012.06.014

Brennan

Chick

Davies

. A taxonomy of model structures for economic evaluation of health technologies. Health Econ. 2006;15(12):1295–310.

Shappell

Caffo

Pekar

Lindquist

. Improved state change estimation in dynamic functional connectivity using hidden semi-Markov models. NeuroImage. 2019;191:243–57. DOI: 10.1016/j.neuroimage.2019.02.013

Rousseau

Mengersen

. Asymptotic behaviour of the posterior distribution in overfitted mixture models. J R Stat Soc Ser B Stat Methodol. 2011;73(5):689–710.

Shum

. Identifiability and Inference of Hidden Markov Models. Technical Report, Citeseer; 2013.

Chrosny

Stevenson

Munzer

. Comparison of Markov and discrete event simulation modeling techniques with application to cost effectiveness analyses. Value Health. 2013;16(7):A587. DOI: 10.1016/j.jval.2013.08.1623

Chhatwal

Jayasuriya

Elbasha

. Changing cycle lengths in state-transition models: challenges and solutions. Med Decis Making. 2016;36(8):952–64.

Smith

. Regenerative stochastic processes. Proc R Soc A. Math Phys Sci. 1955;232(1188):6–31.

Janssen

Limnios

. Semi-Markov Models and Applications. New York: Springer Science & Business Media; 2013.

10.

Barbu

Limnios

. Semi-Markov Chains and Hidden Semi-Markov Models toward Applications: Their Use in Reliability and DNA Analysis. Vol. 191. New York: Springer Science & Business Media, 2009.

11.

Hoem

. Inhomogeneous semi-Markov processes, select actuarial tables, and duration-dependence in demography. In: Greville TNE, eds. Population Dynamics. Philadelphia: Elsevier; 1972. p :251–96.

12.

Commenges

. Multi-state models in epidemiology. Lifetime Data Anal. 1999;5:315–27.

13.

Cao

Buskens

Feenstra

Jaarsma

Hillege

Postmus

. Continuous-time semi-markov models in health economic decision making: an illustrative example in heart failure disease management. Med Decis Making. 2016;36(1):59–71. DOI: 10.1177/0272989x15593080

14.

Latimer

. Survival analysis for economic evaluations alongside clinical trials—extrapolation with Patient-Level Data: Inconsistencies, Limitations, and a Practical Guide. Medl Decis Making. 2013;33(6):743-754. DOI: 10.1177/0272989X12472398

15.

Bell Gorrod

Kearns

Stevens

, et al. A review of survival analysis methods used in NICE technology appraisals of cancer treatments: consistency, limitations, and areas for improvement. Med Decis Making. 2019;39(8):899–909. DOI: 10.1177/0272989X19881967

16.

Buchardt

Møller

Schmidt

. Cash flows and policyholder behaviour in the semi-Markov life insurance setup. Scand Actuar J. 2015;2015(8):660–88.

17.

Helwich

. Durational Effects and Non-smooth Semi-Markov Models in Life Insurance [doctoral dissertation]. Rostock (Germany): University of Rostock; 2008.

18.

Christiansen

Multistate Models in Health Insurance [PhD thesis]. Ulm (Germany): Institut für Versicherungswissenschaften, Universität Ulm; 2010. Available from: https://www.uni-ulm.de/fileadmin/website_uni_ulm/mawi2/dokumente/preprint-server/2010/1015_multmodinhealthins.pdf

19.

Limnios

Oprisan

. Semi-Markov Processes and Reliability. New York: Springer Science & Business Media; 2012.

20.

Castelli

Combescure

Foucher

Daures

. Cost-effectiveness analysis in colorectal cancer using a semi-Markov model. Stat Med. 2007;26(30):5557–71. DOI: 10.1002/sim.3112

21.

Bullement

Cranmer

Shields

. A review of recent decision-analytic models used to evaluate the economic value of cancer treatments. Appl Health Econ Health Policy. 2019;17(6):771–80. DOI: 10.1007/s40258-019-00513-3

22.

Król

Saint-Pierre

. SemiMarkov: an R package for parametric estimation in multi-state semi-Markov models. J Stat Softw. 2015;66:1–16. DOI: 10.18637/jss.v066.i06

23.

Zucchini

MacDonald

Langrock

. Hidden Markov Models for Time Series: An Introduction Using R. Boca Raton (FL): CRC Press; 2017.

24.

Maya

BIG

Limnios

. Using semi-Markov chains to solve semi-Markov processes. Methodol Comput Appl Probab. 2021;23:1419–31.

25.

Davis

Stevenson

Tappenden

Wailoo

. NICE DSU technical support document 15: cost-effectiveness modelling using patient-level simulation. Report by the NICE Decision Support Unit; 2014.

26.

Caro

. Pharmacoeconomic analyses using discrete event simulation. Pharmacoeconomics. 2005;23:323–32.

27.

Degeling

Franken

May

, et al. Matching the model with the evidence: comparing discrete event simulation and state-transition modeling for time-to-event predictions in a cost-effectiveness analysis of treatment in metastatic colorectal cancer patients. Cancer Epidemiol. 2018;57:60–7. DOI: 10.1016/j.canep.2018.09.008

28.

Degeling

IJzerman

Groothuis-Oudshoorn

CGM

, et al. Comparing modeling approaches for discrete event simulations with competing risks based on censored individual patient data: a simulation study and illustration in colorectal cancer. Value Health. 2022;25(1):104–15. DOI: 10.1016/j.jval.2021.07.016

29.

Çağlayan

Terawaki

Chen

Rai

Ayer

Flowers

. Microsimulation modeling in oncology. JCO Clin Cancer Inform. 2018;2:1–11. DOI: 10.1200/CCI.17.00029

30.

Barton

Bryan

Robinson

. Modelling in the economic evaluation of health care: selecting the appropriate approach. J Health Serv Res Policy. 2004;9(2):110–8.

31.

Groot Koerkamp

Stijnen

Weinstein

Hunink

. The combined analysis of uncertainty and patient heterogeneity in medical decision models. Med Decis Making. 2011;31(4):650–61.

32.

Claxton

Sculpher

McCabe

, et al. Probabilistic sensitivity analysis for NICE technology assessment: not an optional extra. Health Econ. 2005;14(4):339–47. DOI: 10.1002/hec.985

33.

Soares

Canto

Castro

. Continuous time simulation and discretized models for cost-effectiveness analysis. Pharmacoeconomics. 2012;30:1101–17.

34.

Clark

Bradburn

Love

Altman

. Survival analysis part I: basic concepts and first analyses. Br J Cancer. 2003;89(2):232–8. DOI: 10.1038/sj.bjc.6601118

35.

Hosmer

Jr Lemeshow

May

. Applied Survival Analysis: Regression Modeling of Time-To-Event Data. New York: John Wiley & Sons; 2011.

36.

Asanjarani

Liquet

Nazarathy

. Estimation of semi-Markov multi-state models: a comparison of the sojourn times and transition intensities approaches. Int J Biostat. 2021;18(1):243–62. DOI: 10.1515/ijb-2020-0083

37.

Hoem

. Markov chain models in life insurance. Blätter der DGVFM. 1969;9(2):91–107.

38.

Latimer

. Survival analysis for economic evaluations alongside clinical trials—extrapolation with patient-level data: inconsistencies, limitations, and a practical guide. Med Decis Making. 2013;33(6):743–54. DOI: 10.1177/0272989X12472398

39.

Kaplan

Meier

. Nonparametric estimation from incomplete observations. J Am Statist Assoc. 1958;53(282):457–81.

40.

Goel

Khanna

Kishore

. Understanding survival analysis: Kaplan-Meier estimate. Int J Ayurveda Res. 2010;1(4):274. DOI: 10.4103/0974-7788.76794

41.

Cox

Oakes

. Analysis of Survival Data. Vol. 21. Boca Raton (FL): CRC Press; 1984.

42.

Kolmogoroff

. Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung. Math Ann. 1931;104:415–58.

43.

Wallace

Preen

Short

Adams

Jeffrey

. Hepatocellular carcinoma in Australia 1982–2014: increasing incidence and improving survival. Liver Int. 2019;39(3):522–30. DOI: 10.1111/liv.13966

44.

Luo

O’Connell

, et al. Cancer incidence and mortality in Australia from 2020 to 2044 and an exploratory analysis of the potential effect of treatment delays during the COVID-19 pandemic: a statistical modelling study. Lancet Public Health. 2022;7(6):e537–48. DOI: 10.1016/S2468-2667(22)00090-1

45.

Huang

Singal

Kanwal

, et al. Hepatocellular carcinoma surveillance — utilization, barriers and the impact of changing aetiology. Nat Rev Gastroenterol Hepatol. 2023;20(12):797–809. DOI: 10.1038/s41575-023-00818-8

46.

Gupta

Rana

Agarwal

, et al. Assessing the risk of further decompensation and survival in patients with cirrhosis with variceal bleeding as their first decompensation event. Am J Gastroenterol. 2022;118(5):833–9. DOI: 10.14309/ajg.0000000000002018

47.

Fleming

Aithal

Card

West

. The rate of decompensation and clinical progression of disease in people with cirrhosis: a cohort study. Aliment Pharmacol Ther. 2010;32(11–12):1343–50. DOI: 10.1111/j.1365-2036.2010.04473.x

48.

Jepsen

Vilstrup

Andersen

. The clinical course of cirrhosis: the importance of multistate models and competing risks analysis. Hepatology. 2015;62(1):292–302. DOI: 10.1002/hep.27598

49.

Georgoulas

Hillston

Sanguinetti

. Unbiased Bayesian inference for population Markov jump processes via random truncations. Stat Comput. 2017;27:991–1002. DOI: 10.1007/s11222-016-9667-9

50.

Lee

Wang

. Statistical Methods for Survival Data Analysis. Vol. 476. New York: John Wiley & Sons; 2003.

51.

Medhi

. Stochastic Processes. New Delhi: New Age International; 1994.

52.

Nahhas

. Introduction to Regression Methods for Public Health Using R. Boca Raton (FL): CRC Press; 2024.

53.

Seaman

Presanis

Jackson

. Estimating a time-to-event distribution from right-truncated data in an epidemic: a review of methods. Stat Methods Med Res. 2022;31(9):1641–55.

54.

Tierney

Stewart

Ghersi

Burdett

Sydes

. Practical methods for incorporating summary time-to-event data into meta-analysis. Trials. 2007;8:1–16.

55.

Hirschmann

. The Semi-Markov Process. Generalizations and Calculation Rules for Application in the Analysis of Systems. Technical report Eidgenoessisches Inst. fuer Reaktorforschung; 1983.

56.

Pinsky

. An early-and late-stage convolution model for disease natural history. Biometrics. 2004;60(1):191–8. DOI: 10.1111/j.0006-341X.2004.00023.x

57.

Flaxman

Mishra

Gandy

, et al. Estimating the effects of non-pharmaceutical interventions on COVID-19 in Europe. Nature. 2020;584(7820):257–61. DOI: 10.1038/s41586-020-2405-7

58.

van Wijk

Simonsson

USH

. Finding the right hazard function for time to event modelling: a tutorial and Shiny application. CPT Pharmacometrics Syst Pharmacol. 2022;11:991–1001. DOI: 10.1002/psp4.12797

59.

Juckett

Rosenberg

. Comparison of the Gompertz and Weibull functions as descriptors for human mortality distributions and their intersections. Mech Ageing Dev. 1993;69(1–2):1–31. DOI: 10.1016/0047-6374(93)90068-3

60.

Kleinbaum

Klein

. Parametric survival models. In: Survival Analysis: A Self-Learning Text. New York: Springer; 2012. p 289–361.

61.

Wang

J-L

. Smoothing hazard rates. In: Armitage

Colton

, eds. Encyclopedia of Biostatistics. Vol. 7. New York: Wiley; 2005. p :4986–97.

62.

Guo

. Cause-Specific Analysis of Competing Risks Using the PHREG Procedure. Paper No.: SAS2159-2018. Cary (NC): SAS Institute Inc.; 2018.

63.

Fine

Gray

. A proportional hazards model for the subdistribution of a competing risk. J Am Stat Assoc. 1999;94(446):496–509. DOI: 10.2307/2670170

64.

Coldman

Phillips

Brisson

, et al. Using the cancer risk management model to evaluate colorectal cancer screening options for Canada. Curr Oncol. 2015;22(2):41–50. DOI: 10.3747/co.22.2013

65.

Crowther

Lambert

. Parametric multistate survival models: flexible modelling allowing transition-specific distributions with application to estimating clinically useful measures of effect differences. Stat Med. 2017;36(29):4719–42.

66.

World Health Organisation. WHO Methods and Data Sources for Life Tables 1990–2019. Geneva: World Health Organisation; 2020.

67.

Australian Bureau of Statistics. Life Expectancy Methodology. Technical report. Sydney: Australian Bureau of Statistics; 2022. Available from: https://www.abs.gov.au/methodologies/life-expectancy-methodology/2020-2022

68.

Bebbington

Lai

C-D

Murthy

DNP

Zitikis

. Rational polynomial hazard functions. Int J Performability Eng. 2010;6(1):35.

69.

Gold

Stevenson

Fryback

. HALYS and QALYS and DALYS, oh my: similarities and differences in summary measures of population health. Annu Rev Public Health. 2002;23(1):115–34.

70.

Hethcote

. Three basic epidemiological models. In: Levin

Hallam

Gross

, eds. Applied Mathematical Ecology. Berlin: Springer; 1989. p :119–44.

71.

Hutchison

Shapiro

. Lead time gained by diagnostic screening for breast cancer. J Natl Cancer Inst. 1968;41(3):665–81.

72.

Jansen

Alexander

Anderson

Church

. Quantifying lead-time bias in risk factor studies of cancer through simulation. Ann Epidemiol. 2013;23(11):735–41. DOI: 10.1016/j.annepidem.2013.07.021

73.

Rosner

Broemeling

. Bayesian inference for the lead time in periodic cancer screening. Biometrics. 2007;63(3):873–80.

74.

Draisma

van Rosmalen

. A note on the catch-up time method for estimating lead or sojourn time in prostate cancer screening. Stat Med. 2013;32(19):3332–41.

75.

Muehleisen

Bergerson

. Bayesian calibration—what, why and how. In: International High Performance Buildings Conference. West Lafayette, Indiana, USA. 10–14 July, 2022.

76.

Wade

Sarich

Vaneckova

, et al. Using Bayesian evidence synthesis to quantify uncertainty in population trends in smoking behaviour. Stat Methods Med Res. 2025;34(3):545–560.

77.

Hong

Gow

Fink

, et al. Surveillance improves survival of patients with hepatocellular carcinoma: a prospective population-based study. Med J Aust. 2018;209(8):348–54. DOI: 10.5694/mja18.00373

78.

Cancer Council Australia Hepatocellular Carcinoma Surveillance Working Group. Clinical Practice Guidelines for Hepatocellular Carcinoma Surveillance for People at High Risk in Australia. Sydney: Cancer Council Australia; 2023. Available from: https://www.cancer.org.au/clinical-guidelines/liver-cancer/hepatocellular-carcinoma

79.

Worthington

Caruana

, et al. An economic evaluation of routine hepatocellular carcinoma surveillance for high-risk patients. medRxiv. 2024. Available from: https://www.medrxiv.org/content/10.1101/2024.02.20.24303111v1

80.

Vilar-Gomez

Calzadilla-Bertot

Wai-Sun Wong

, et al. Fibrosis severity as a determinant of cause-specific mortality in patients with advanced nonalcoholic fatty liver disease: a multi-national cohort study. Gastroenterology. 2018;155(2):443–57. DOI: 10.1053/j.gastro.2018.04.034

81.

D’Amico

Garcia-Tsao

Pagliaro

. Natural history and prognostic indicators of survival in cirrhosis: a systematic review of 118 studies. J Hepatol. 2006;44(1):217–31. DOI: 10.1016/j.jhep.2005.10.013

82.

Australian Bureau of Statistics. Life Tables. Australian Bureau of Statistics; 2022. Available from: https://www.abs.gov.au/statistics/people/population/life-tables/latest-release#cite-window1

83.

Koffijberg

Degeling

IJzerman

Coupé

VMH

Greuter

MJE

. Using metamodeling to identify the optimal strategy for colorectal cancer screening. Value Health. 2021;24(2):206–15. DOI: 10.1016/j.jval.2020.08.2099

84.

Degeling

IJzerman

Koffijberg

. A scoping review of metamodeling applications and opportunities for advanced health economic analyses. Expert Rev Pharmacoecon Outcomes Res. 2019;19(2):181–7. DOI: 10.1080/14737167.2019.1548279

85.

Vahdat

Alagoz

Chen

Saoud

Borah

Limburg

. Calibration and validation of the Colorectal Cancer and Adenoma Incidence and Mortality (CRC-AIM) microsimulation model using deep neural networks. Med Decis Making. 2023;43(6):719–36. DOI: 10.1177/0272989X231184175

86.

Pineda-Antunez

Seguin

van Duuren

, et al. Emulator-based Bayesian calibration of the CISNET colorectal cancer models. Med Decis Making. 2024;44(5):543–53. DOI: 10.1177/0272989x241255618

87.

Standfield

Comans

Scuffham

. Markov modeling and discrete event simulation in health care: a systematic comparison. Int J Technol Assess Health Care. 2014;30(2):165–72. DOI: 10.1017/S0266462314000117

88.

Delgado-Vences

Flandoli

. A spectral-based numerical method for Kolmogorov equations in Hilbert spaces. Infin Dimension Anal Quantum Probab Relat Top. 2016;19(3):1650020.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.03 MB