A dynamic Bayesian Markov model for health economic evaluations of interventions in infectious disease.

BACKGROUND: Health economic evaluations of interventions in infectious disease are commonly based on the predictions of ordinary differential equation (ODE) systems or Markov models (MMs). Standard MMs are static, whereas ODE systems are usually dynamic and account for herd immunity which is crucial to prevent overestimation of infection prevalence. Complex ODE systems including distributions on model parameters are computationally intensive. Thus, mainly ODE-based models including fixed parameter values are presented in the literature. These do not account for parameter uncertainty. As a consequence, probabilistic sensitivity analysis (PSA), a crucial component of health economic evaluations, cannot be conducted straightforwardly.
METHODS: We present a dynamic MM under a Bayesian framework. We extend a static MM by incorporating the force of infection into the state allocation algorithm. The corresponding output is based on dynamic changes in prevalence and thus accounts for herd immunity. In contrast to deterministic ODE-based models, PSA can be conducted straightforwardly. We introduce a case study of a fictional sexually transmitted infection and compare our dynamic Bayesian MM to a deterministic and a Bayesian ODE system. The models are calibrated to simulated time series data.
RESULTS: By means of the case study, we show that our methodology produces outcome which is comparable to the "gold standard" of the Bayesian ODE system.
CONCLUSIONS: In contrast to ODE systems in the literature, the dynamic MM includes distributions on all model parameters at manageable computational effort (including calibration). The run time of the Bayesian ODE system is 15 times longer.

Background

Vaccines, antibiotics and antivirals against infectious diseases offer health benefits to society [1, 2] and have been instrumental in the prevention and treatment of conditions previously causing egregious burden to public health. Examples include the extremely low prevalence of syphillis, the control of human immunodeficiency virus, the worldwide eradication of smallpox [3] and the extremely low incidence of tetanus, diphteria and congenital rubella syndrome in the Western world [4]. However, despite being frequently successful from a clinical point of view, vaccination programmes and antiretrovirals are often costly to apply. As a pre-requisite to their implementation, health interventions such as vaccines are thus increasingly subject to cost-effectiveness analyses (CEAs) [5, 6].

The National Institute for Health and Care Excellence (NICE) is arguably the leading health technology assessment agency in the world. In the UK, NICE is responsible for providing guidance and advice on whether proposed interventions should be publicly funded. NICE has developed a set of criteria and guidelines that drive the analytic process of CEA [7]. Crucially, these involve the explicit necessity of assessing the impact of parameter uncertainty on the decision making outcome, a process typically known as Probabilistic Sensitivity Analysis (PSA) [6, 8, 9].

In the UK, the appraisal of vaccines falls under the remit of the Joint Committee for Vaccines and Immunisations, an independent expert advisory committee to the ministers and health departments. Since 2009, the Health Protection Regulation obliges the Secretary of State to ensure that recommendations for national vaccination programmes are based on an assessment demonstrating cost-effectiveness [10]. However, there are currently no vaccine-specific guidelines for developing clinical or cost-effectiveness evidence.

One of the reasons for this circumstance is perhaps the intrinsic complexity of infectious disease modelling, which is typically performed through compartmental models. These are highly complicated mathematical tools capable of simulating the natural history of disease infection and progression. More specifically, in pathogens transmissible among humans, these models need to account for population dynamics and herd immunity [11]. Herd immunity implies that due to lower infection prevalence, the introduction of preventive and therapeutic measures such as vaccination, quarantine, antivirals and antibiotics induces a reduced risk of pathogen exposure. Only dynamic models are able to prevent incorrect predictions since they are suitable to incorporate these effects [12-14].

Dynamic compartmental models are commonly fitted by solving systems of Ordinary Differential Equations (ODEs) in continuous time. While these deterministic models usually deal with features such as herd immunity (and thus are considered the “industry standard” in infectious disease modelling), they are characterised by a notable computational effort. One important consequence is that, in most cases, epidemiological and economic modelling for infectious disease performed by means of ODEs is based on the inclusion of fixed, predefined values on the model parameters. These fixed values are usually informed through a point estimate. The joint uncertainty in the parameters is then not considered; these models result in outcome (e.g. on the number of people in the states) which does not include distributions.

Therefore, PSA on the model outcome can only be conducted in retrospect and not in a straightforward way. An additional step using Latin Hypercube Sampling or Monte Carlo sampling is necessary, as shown in [15-17] and [18-20], respectively. Alternative methods that may prove computationally efficient when estimated through polynomial chaos expansions as shown in [21, 22] are provided by the Sobol and Fourier amplitude sensitivity test (FAST) indices. These indices are based on ANOVA techniques and thus estimate the total contribution of each model parameter or a combination of parameters to the variance of the model output [23, 24]. However, in contrast to a full PSA, uncertainty is not propagated through the whole model.

The computational feasibility of PSA in retrospect is limited in models which include a high number of states and model parameters. In addition, in contrast to a Bayesian approach, parameter uncertainty is not propagated through the crucial model parts of pathogen transmission and disease progression. The outcome based on fixed parameter values can differ considerably from the PSA outcome, and the two are in most cases not reported in enough detail to identify possible inconsistencies (with the exception of [19]). The results presented commonly focus on uncertainty in the health economic rather than the prevalence outcome of the models, and PSA results on infection prevalence including the corresponding confidence intervals are often not given [18, 20, 25]. This approach is highly questionable, especially with respect to consistency with validation targets. To ensure that the model outcome on the number of people in the states and infection prevalence is realistic, calibration to high quality data based on large sample size is necessary. This is often not conducted [25], or only conducted on the outcome based on fixed parameter values [18]; the PSA outcome on prevalence is usually not evaluated with respect to fitting high quality data [19, 20, 26].

The more complex an ODE system especially with respect to state space and number of model parameters, the larger the effort on implementation and computation, especially if each model parameter is assigned a suitable distribution. This might be one of the reasons why the International Society for Pharmacoeconomics and Outcomes Research guideline for best modelling practice in infectious disease suggests that PSA is not a fundamental component of health economic assessment [12]. This recommendation is given in contrast with NICE and virtually any other disease area. As a consequence, most economic models for vaccines only consider deterministic sensitivity analysis, which is based on selecting a grid of “plausible” values for a subset of model parameters in order to assess the robustness of the decision-making process. This approach is however not recommended in general, as it fails to account for potential correlation among the parameters [6, 9, 27].

In contrast to ODE-based models, systems of equations can also be defined in discrete-time, which are termed difference equation models [28]. An alternative compartmental specification is given by Markov models (MMs). MMs are used to model progression over time across a finite set of states. Since MMs are described by a stochastic process, these are classified as stochastic models. This is in contrast to ODE systems and difference equation models, which belong to the class of deterministic models. In a deterministic model, the same set of parameter values and initial conditions always results in the same output. In contrast, a stochastic model produces different output each time the model is run, accounting for randomness. Apart from the model class, difference equation models and discrete-time MMs are mathematically comparable.

Although MMs can also be computationally intensive, it is generally feasible to implement even complex models in a Bayesian framework or to use re-sampling methods such as the bootstrap to characterise the uncertainty in the model parameters. Perhaps for this reason, MMs are a very popular tool in health economic evaluation. Nevertheless, a major limitation in infectious disease modelling is that they are intrinsically static, i.e. they do not account for population dynamics [29].

We introduce in this paper an extension to standard MMs, which we term “dynamic Bayesian MM” to indicate that we consider a stochastic model and use a Bayesian framework to estimate its underlying parameters. We directly include the force of infection of the pathogen, which automatically accounts for time-dependent changes in prevalence and thus the effects of herd immunity, into the state allocation algorithm of a MM. In other words, the movement of susceptibles to the state of infection is directly represented by the dynamic force of infection. A direct inclusion of the force of infection into the state allocation algorithm of a difference equation model [30], a MM [31, 32] or its direct consideration in a model based on a stochastic process [33-35] was presented previously by several authors. However, the authors who present stochastic models do not conduct a health economic evaluation. In the health economics literature, to the best of our knowledge, no approach of a dynamic MM including a high number of states and suitable probability distributions on all model parameters is presented.

Our dynamic Bayesian MM combines six advantages in comparison to the few dynamic MMs presented in the literature. Firstly, in contrast to our contribution, the compartmental models in the literature are only suitable to include a low number of states due to computational limitations and the majority consist of no more than four states (apart from [30]). Our methodology is especially suitable to incorporate an extensive number of states as described elsewhere [36] for the application to human papillomavirus (HPV) modelling. In our model on HPV, 36 states in females and 22 states in males are included to account for all known HPV-induced diseases apart from recurrent respiratory papillomatosis. In addition, statisticians or health economic modellers can implement our approach directly in the commonly used software R linking to JAGS or WinBUGS; the corresponding run-time is considerably fast, and therefore, it is not necessary to use a compiled language such as C which is usually not widely used in this field. The four remaining advantages are given through the Bayesian framework, which i) is highly flexible with its probabilistic nature since it considers multiple sources of prior information in terms of evidence synthesis [37], ii) enables propagation of parameter uncertainty through the infection transmission, progression and economic evaluation process, iii) ensures that calibration targets are met through a constant updating process of the outcome on the numbers of people in the states directly in the state allocation algorithm, using available time series data, and iv) simplifies the process of PSA, an essential part of CEAs, avoiding the necessity of applying additional sampling techniques such as Monte Carlo or Latin Hypercube Sampling once the model output is available.

The paper is structured as follows: Firstly, we describe compartmental models in widespread use for CEAs on interventions against infectious diseases. Secondly, we introduce our contribution of the dynamic Bayesian MM. Thirdly, we compare the performance of an ODE system including fixed values on the model parameters and an ODE system in a Bayesian framework to our methodology, using a case study of a chronic sexually transmitted infection. To contrast the three methodologies in practice, we compare the natural history of disease following calibration, and conduct CEAs including PSA, comparing a screening strategy to a hypothetical vaccine. Finally, we discuss advantages and disadvantages of ODE-based methodology in comparison to the dynamic Bayesian MM.

Methods

Compartmental models consist of a set of mutually exclusive and exhaustive states describing disease infection and progression. We indicate the elements of as s=1,…,S. Members of a “virtual” population move across the states over a pre-specified time horizon.

Figure 1 shows an example of a compartmental model incorporating the natural disease history of a chronic sexually transmitted infection (STI) with S=5 states which is deemed similar to HIV [38]. The assumptions encoded by this structure are that the whole population initially is in the state Susceptible (indexed by s=1), from which a proportion can move to the state Infected (s=2). Following this, people move to an Asymptomatic state (s=3). A progression to the state Morbid (s=4) induces the development of disease symptoms. The state Dead (s=5) can be reached from any state; people die due to any cause or as a consequence of being in the state Morbid. Compared to the average population, the latter have a higher risk of death. A transition from one state to another is defined according to transition parameters [39]. They are indicated as ϕ, where represent the original and target state, respectively. People proliferate at a rate χ, resulting in a replenishment of the pool of susceptibles at risk of contracting the infection.

Fig. 1

Model structure of a hypothetical chronic sexually transmitted infection. The arrows represent the possible transitions. These are governed by the parameters ϕ with indices representing origin and target states, respectively. The replenishment of the pool of susceptibles by newborns proceeds at a rate χ

Compartmental models may differ in two characteristics. The first is the specification of time. The most realistic option is to allow transitions among the states at any point in time t>0; this is a so-called “continuous-time approach”. Alternatively, it is possible to assume that transitions occur in discrete time where only one transition is possible within a pre-defined time interval , where κ determines the corresponding interval width, commonly referred to as cycle. Depending on the medical context, κ can be specified in terms of daily, weekly, monthly or yearly cycles. The second difference concerns the way in which population dynamics are considered: models including a force of infection which accounts for changes in prevalence are referred to as dynamic, while those that include a fixed force of infection and thus ignore the effects of herd immunity are termed static.

In addition, different approaches to model parameter specification exist, which may have major impact on the outcome of health economic evaluations. Depending on the methodology used, the induced computational effort might not allow the inclusion of probability distributions on all model parameters. In that case, the model parameters are fixed values. Commonly, these are estimated through a relevant summary, e.g. mean, median or mode, obtained from available data. The point estimate is then used as a plug-in for the corresponding parameter. In contrast, if whole distributions of values are assigned, parameter uncertainty is propagated through the infection progression. While frequentist versions of this strategy exist (e.g. based on bootstrap), this type of modelling is most naturally handled within a Bayesian paradigm.

Ordinary differential equation models

ODE systems model the rate of change in the number of people within a given state in continuous time; thus, the corresponding parameters are transition rates and we denote them as ρ(t), with representing again the origin and target states, respectively. In principle, transition rates can depend on t, but do not necessarily have to. The number of people transitioning in each state at t is multiplied by the corresponding transition rates to obtain the inflow and outflow to and from a state. The difference between the number of people entering and leaving a state corresponds to the derivative of the number of those in the respective state.

Back to our example, we define the vector (t)=(n1(t),…,n(t))′, where n(t) is the number of people in state s at time t. The corresponding ODE system is given by the set of equations

The rate of change in the number of people in each state at each point in time t is subject to population dynamics and exposure to sources of infection. The transition to the state of infection is determined by the dynamic force of infection of the pathogen, indicated by ρ1,2(t) in (1). This is a function of the probability of pathogen transmission, partner acquisition rates and prevalence.

If each model parameter is assigned a fixed value, parameter uncertainty is not accounted for. Scenario analyses are often performed, for example by estimating the parameters through summary statistics at the extremes of the corresponding parameter distribution (e.g. lower or upper quantiles, minima or maxima). As mentioned earlier, this is not equivalent to the application of a full PSA.

Theoretically, it is possible to incorporate a probability distribution for each parameter of an ODE-based model, for example in a Bayesian context. The uncertainty is then propagated through the estimation procedure, which again generates a full distribution of outcomes. This type of model can be analysed using for instance software based on Markov Chain Monte Carlo (MCMC) algorithms such as WinBUGS [40] or Stan [41], a very promising tool, which in general performs extremely well with relatively complex systems. Both include ODE solvers and can be linked to the statistical programming language R.

In realistic problems including a large number of states and complex structures, assigning a suitable distribution to all model parameters in ODE systems may be impractical since the model needs to be run for a large number of simulations to ensure convergence of each parameter and thus the ODEs have to be solved repeatedly for each parameter combination. The increase in the computational time is induced by the length of the observation time horizon, the amount of parameters, the complexity of contacts, and most importantly by the number of states which increases the number of differential equations. Consequently, complex ODE-based models which incorporate distributions on all parameters are rare exceptions in the literature on infectious disease transmission modelling [42].

Discrete-time Markov models

The main characteristic of MMs is the Markov assumption which implies that the transition to a future state is exclusively conditional on the current state or on a limited set of previous states, but not on the full trajectory. However, the Markov assumption can be relaxed by accounting for covariates (e.g. age and sex) or for time-specific prevalence and population dynamics in the transition parameters.

ODE systems and MMs differ in the way they describe the process of transitions. As suggested earlier, in the former, the rates of change are calculated dynamically through differentiation, while in the latter, the transitions are described by a static Markov process (a category of a stochastic process). As for ODE-based models, MMs can be implemented for continuous time.

However, the vast majority of MMs in the health economic literature is based on a discrete-time approach [43]. In this case, members of the population move across the states according to a set of transition probabilities only once per time interval (termed “Markov cycle”). These probabilities can be arranged in a matrix =(π), whose elements represent the transition probabilities for movements from an original state r to a target state s.

For the model structure of Fig. 1, the transition probability matrix is defined as implying that, for example, a susceptible either acquires the infection (with probability π1,2), dies (with probability π1,5), or remains susceptible, which occurs with probability π1,1=1−π1,2−π1,5.

If we define the vector =(n1,…,n)′, where n is the number of people in state s and at each time interval , then transitions across the states from one time interval to the next are calculated as

MMs are relatively straightforward to implement and are commonly used to model the progression of non-communicable conditions such as cardiovascular disease and cancer. Therefore, they are established in the health economic literature and well-known to clinicians and decision makers. However, the process of pathogen transmission is not estimated correctly using standard MMs. A transition of susceptibles to the state of infection is commonly represented by a static transition parameter which does not consider changes in the population prevalence over time. These especially occur after the introduction of a preventive intervention such as vaccination into a fully susceptible population.

The predictions of static MMs on population prevalence are commonly incorrect (although notable exceptions include scenarios with very low vaccine coverage or pathogens that cannot be transmitted between humans, e.g. tetanus). In the worst case, the whole model outcome on infection prevalence and the related CEA can be incorrect, e.g. because of the impact of an unrecognised shift in the age of infection of childhood diseases. Some childhood diseases are relatively harmless in young children but prone to lead to serious health issues in adults. Incorrect predictions of static MMs on population health and induced costs, e.g. through hospitalisation and treatment, can have dire consequences [29].

As for ODE-based models, a dynamic force of infection could be incorporated into the transition probabilities to account for the effects of herd immunity. However, to the best of our knowledge, dynamic Markov models are not commonly used in the health economics literature.

Dynamic Bayesian Markov models

To overcome the limitations discussed above and with a view to extending the modelling framework for health economic evaluation of interventions in infectious disease, our main idea is to add a force of infection which depends on population dynamics and prevalence into a MM setting. As a consequence, the transition probabilities from the state Susceptible to the state Infected are directly defined by the dynamic force of infection. Specifically, we set up our model so that the force of infection is calculated separately within each cycle of the state allocation algorithm corresponding to (2) as a function of

the probability of pathogen transmission per contact, which we indicate as β;

the rate of contacts between susceptible and infectious members of the population ω; and

the time-dependent pathogen prevalence where I represents the number of people in the state of infection and, assuming that state S indicates death, is the number of those alive at time interval .

The force of infection is recalculated at each Markov cycle as

Since ω is a rate, (3) also results in a transition rate. Assuming that λ remains constant within each time interval, the corresponding time-dependent transition probability for the discrete-time MM is estimated as

We acknowledge that the estimation of the force of infection may be only approximate, due to the competing risk of death and the assumption of uniformity within the intervals . This assumption is not likely to hold if the disease is characterised by very fast transmission, or when events associated with the infection are likely to occur in short periods of time. In these cases, it is perhaps advisable to reduce the length κ of the cycles and the duration of the follow-up, to avoid unrealistic estimates for the number of subjects in the states per cycle. For example, for yearly cycle length, immediate death following infection would be delayed by up to one year. These delays would then accumulate through the whole model and introduce a more substantial bias at later follow-up.

However, a likelihood function for interval censored data could be estimated to account for competing risks as shown in [44, 45] for an independent inspection process (IPP) model; this model allows for future movements to be conditional on the history of the data. In Lemma 1, Hudgens et al. [44] present this likelihood under the assumption that only one event can occur per cycle. As an alternative to (4), this likelihood could be derived for the cumulative distribution function of an exponential model.

Moreover, competing risks could be considered through the Kolmogorov forward equations. This alternative approach would allow for the possibility that more than one event could occur per cycle; for instance, one could acquire the infection followed by death within the same cycle. In order to move from the state Susceptible to Infected, one would have to take into account that the infection was acquired and death did not occur within the remainder of that cycle. For the exponential distribution, the corresponding Kolmogorov forward equation to estimate the transition probability of moving between the two states was derived as

accounting for ρ1,5(t), the transition rate of moving from the state Susceptible to Death. Replacing (4) by (5), the accuracy of the approximation of π1,2,t could improve. Notice, however, that in the case study presented in Section “Case study”, the probability of death is very low in susceptibles (who are by definition in good health and of young age) and therefore this additional complication in the estimation of the transition probability is not necessary. The results of the corresponding discrete- and continuous-time models are comparable; the approximation in (4) is therefore deemed sufficient.

The transition probability π1,2,t approximated by (4) is multiplied by the proportion of the population in the state Susceptible to provide an estimation of the contingent of movements to the state Infected, effectively including dynamic, time-dependent changes in prevalence in the corresponding transitions.

The computational effort is reduced by fitting models that do not involve complex ODEs, while still allowing for mixing patterns within the population. Another potential advantage of the dynamic MM framework is that it is fairly simple to incorporate suitable probability distributions on all model parameters, even if the model is complex with an extremely large number of parameters and states. In contrast, the related computational effort in a comparable ODE-based model would be extremely high due to i) numerical integration, e.g. through the Runge-Kutta solver [46]; ii) considerably smaller step sizes of ODE solvers when compared to the cycle length of discrete-time models; and iii) accounting for competing risks in the transition rates as elaborated above. The model can be easily extended to include a high number of age cohorts for infectious diseases with age-specific prevalence such as HPV; we present this application including 24 age cohorts at manageable computational effort elsewhere [36].

Accounting for parameter uncertainty is particularly relevant because, for obvious ethical and practical reasons, it is invariably difficult (if possible at all) to obtain and use experimental evidence to inform the pathogen transmission probability β and the active contact rate ω — arguably the crucial parameters. Often observational studies or expert opinions are the only available information with the consequence that large uncertainty remains over the most likely range, let alone the “true” value of the parameters. A Bayesian approach may provide great benefit in allowing this uncertainty to be fully propagated and perhaps in integrating different sources of evidence (e.g. using evidence synthesis [39]); this indeed has been advocated for MMs in the health economics literature [9, 39, 47]. In [48], β is assigned a Uniform distribution in the interval [0;1], which essentially amounts to allowing any value (between 0 and 100%) as equally possible. We have performed extensive sensitivity analysis to this parameter in the HPV model and found that the prevalence outcome was highly sensitive to this parameter (results not published).

In a Bayesian dynamic MM setting, it is possible to assign prior distributions to the parameters (β,ω) to represent the state of science — if data are available, these are updated into posterior distributions although it is possible to still propagate uncertainty in the priors even when no data on pathogen transmission or active contacts are observed. In addition, the quantity ψ is estimated for each cycle as a function of transition probabilities, which can also be modelled using suitable distributions. This modelling process induces a probability distribution on ψ and a fortiori also on λ, which is defined as a function of the three random parameters (β,ω,ψ). Thus, the corresponding transition probabilities π1,2,t are modelled probabilistically, meaning that uncertainty in the population dynamics is propagated through the economic model.

If β and ω are varied simultaneously and difficult to inform directly using empirical evidence, the force of infection may be subject to issues of non-identifiability, i.e. it is possible that it cannot be distinguished to what extend differences in the corresponding outcome on the number of subjects in the state Infected are induced by which parameter. Individuals become infected as a consequence of meeting an infected subject (determined by β) once or several times (determined by ω). However, we note here that we are mainly interested in the overall number of subjects in the state Infected, irrespective of the number of contacts that are necessary to result in an infection. Therefore, in the economy of our modelling framework, it is less important if changes in the number of subjects in the state Infected are induced by β or ω. Nevertheless, we do acknowledge this issue and suggest careful consideration of the prior assumptions when investigating sensitivity of the results.

Another crucial aspect in infectious disease modelling (and more generally in statistical analysis) is that of calibration of the model output [49, 50]. The Bayesian framework enables the calibration of the numbers of people in the states directly in the state allocation algorithm, using available time series data for a specific time frame of follow-up. The corresponding details are explained in “Model parameters and related distributions” section and Appendix 2.

The BMM is generalizeable to any infectious disease which is transmitted between humans and where interventions (e.g. quarantine, vaccination, antibiotic or antiviral treatment) are available to reduce infection prevalence. It is especially suitable to include an extremely high number of states, model parameters with suitable distributions, and age cohorts as shown in [36].

To estimate the efficacy and cost-effectiveness of an intervention, only dynamic models predict realistic outcome and account for herd immunity in case of vaccination. Yet, any kind of intervention which reduces prevalence has a protective impact on susceptibles who do not directly receive it; this can be considered in the BMM.

Finally and specifically for the purpose of economic evaluation, the dynamic BMM has the advantage that PSA can be performed “for free”, once the model output is produced. In a Bayesian framework, the MCMC simulations for all the model parameters can be combined to obtain a full characterisation of the uncertainty in the decision-making process. This can be post-processed (e.g. using the R package BCEA [51]) to produce relevant summaries such as the cost-effectiveness plane, the cost-effectiveness acceptability curve and the analysis of the value of information (see “Cost-effectiveness analysis” section).

Our approach to a dynamic BMM provides PSA samples directly as part of the probabilistic model output. Therefore, in addition to simpler probabilistic sensitivity analysis (e.g. based on CEACs), the BMM enables the conduct of extensive Value of Information (VoI) analysis, which can be used to prioritise further research on key model inputs, currently driving uncertainty in the decision making process. For example, recent methods have been developed for fast computation of both the Expected Value of Partial Perfect Information (EVPPI) and the Expected Value of Sample Information (EVSI) [52-54] that use simulations from the model parameters obtained directly during the process of PSA. The BMM conforms with this structure and would thus be suitable for VoI analysis based on these methods, without the need to divorce the transmission model from the economic evaluation, as often happens in current practice.

Case study

We consider again the fictional chronic STI described above and compare the dynamic Bayesian MM to both a deterministic and a Bayesian ODE system. We denote the three models as BMM, dODE and BODE, respectively. We evaluate whether our BMM produces results that are in line with the “gold standard” of the BODE. In the three models, we distinguish between sexes as well as high- and low-risk sexual behaviour. The duration of the follow-up is set at 100 years, with a yearly Markov cycle length. We consider a population size of 1,000,000 and initially assume that 600 people are infected, whereas the remainder are susceptible. Males amount to 50% and the high-risk group to 20% of the population; the sex ratio in the two risk groups is constant. The proportion of infected people in both sexes and risk groups is identical. We account for sex-specific differences in sexual behaviour, assuming higher partner acquisition rates in males. The population size changes due to births and deaths. We conduct our analysis for two competing health-care interventions. We assume that in the status-quo, screening takes place at intervals of five years at a pre-defined rate to enable an early detection of the STI. For simplicity, we assume that under the vaccination scenario no screening takes place. We assume that the vaccine is only effective before initial STI infection; thus, susceptibles are vaccinated at a specified vaccine uptake rate at intervals of five years. Following STI diagnosis, treatment is provided in both interventions.

In the dODE and BODE, the force of infection (transition rate from the state Susceptible to the state Infected) shown in (3) has to be adjusted to account for the covariates sex and behaviour with respect to infection prevalence. In the BMM, the sex- and behavioural-specific transition probability is estimated by transforming this adjusted transition rate according to (4). For simplicity, we exclusively account for random mixing. Further details are given in Appendix 1.

Model parameters and related distributions

In addition to the probability of STI transmission β and the partner acquisition rates ω, the model contains a variety of parameters such as those determining the screening and vaccine coverage, the unit costs of STI diagnostics and treatment and the health utilities, which are relevant in context of the cost-effectiveness analysis. We assign informative priors to transition parameters and costs (defined in monetary units of £) and utilities. We specify the distributional assumptions so that the outputs of the health economic evaluation are within reasonable ranges. We assign informative or minimally-informative priors to the remaining parameters and update these through simulated individual-level and aggregate data into the corresponding posteriors. Using simulated data is reasonable since the case study is fictional. For example, we pretend that data on partner acquisition rates ω are available from a large data registry. To infer ω, we update informative Gamma priors into the corresponding posteriors through Poisson-Gamma models. Beta-Binomial models are used to infer probabilities (e.g. the vaccine coverage α and efficacy γ, which are informed by data). We assume to have access to vaccination data of 500 individuals of whom 450 have received the vaccine. In addition, we assume to have the information that in 450 out of 500 people who received the vaccine, it was effective.

Table 1 shows an overview of the model parameters ={ω,χ,β,π,τ,ξ,γ,σ,,u}. The means and 95% credible intervals are displayed for the BMM (the values in the BODE are comparable) and rounded to two decimal places. However, the 95% CI of the parameter χ is in fact defined as [0.009420;0.010596]. Since the BMM is based on a discrete-time approach, the corresponding transition probabilities are modelled using Beta distributions. In contrast, the transition rates of the continuous-time BODE are modelled using Gamma distributions. The transition probabilities for movements from the states Susceptible, Infected and Asymptomatic to Dead are assumed as identical; thus, only π1,5 is shown.

Table 1

Overview of the informative priors and the models used for updating informative and minimally-informative priors

Parameter	Description	Distribution/model BMM	Distribution/model BODE	Mean	2.5/97.5% percentiles
ω MH	Partner acquisition rate (high-risk males)	Poisson-Gamma model	Equivalent to BMM	9.10	[8.77;9.29]
ω ML	Partner acquisition rate (low-risk males)	Poisson-Gamma model	Equivalent to BMM	2.98	[2.82;3.12]
ω FH	Partner acquisition rate (high-risk females)	Poisson-Gamma model	Equivalent to BMM	9.00	[8.71;9.26]
ω FL	Partner acquisition rate (low-risk females)	Poisson-Gamma model	Equivalent to BMM	1.96	[1.86;2.09]
χ	Proliferation parameter	Gamma(1111.1,111111.1)	Gamma(1111.1,111111.1)	0.01	[0.01;0.01]
β	STI transmission probability per partnership	Beta-Binomial model	Equivalent to BMM	0.16	[0.15;0.16]
π 2,3	Transition parameter from state 2 to state 3	Beta(5119.2, 1279.8)	Gamma(25600,32000)	0.80	[0.79;0.81]
π 3,4	Transition parameter from state 3 to state 4	Beta(1842.66, 18631.34)	Gamma(2025,22500)	0.09	[0.09;0.09]
π 4,5	Transition parameter from state 4 to state 5	Beta(1535.96, 36863.04)	Gamma(1600,40000)	0.04	[0.04;0.04]
π 1,5	Transition parameter from state 1 to state 5	Beta(156.171, 312186.6)	Gamma(156.25,312500)	<0.01	[<0.01;<0.01]
η	Probability of STI diagnosis	Beta-Binomial model	Equivalent to BMM	0.90	[0.88;0.92]
σ	Screening probability	Beta-Binomial model	Equivalent to BMM	0.90	[0.87;0.92]
α	Vaccine coverage parameter	Beta-Binomial model	Equivalent to BMM	0.90	[0.87;0.92]
γ	Vaccine efficacy parameter	Beta-Binomial model	Equivalent to BMM	0.90	[0.87;0.92]
c screen	Unit cost of screening in £	Lognormal(2.996, 0.693)	Equivalent to BMM	25.39	[5.19;77.53]
c vac	Unit cost of vaccination in £	Lognormal(5.011, 0.01)	Equivalent to BMM	150.02	[147.14;152.98]
c test	Unit cost of STI test in £	Lognormal(2.996, 0.03)	Equivalent to BMM	20.01	[18.83;21.19]
c blood	Unit cost of blood test in £	Lognormal(3.401, 0.03)	Equivalent to BMM	30	[28.26;31.79]
c treat	Unit cost of treatment in £	Lognormal(8.517, 0.015)	Equivalent to BMM	4999.78	[4853.56;5149.24]
c dis	Unit cost of disease treatment in £	Lognormal(9.210, 0.01)	Equivalent to BMM	9999.95	[9802.97;10198.10]
c gp	Unit cost of visit to general practitioner in £	Lognormal(3.912, 0.02)	Equivalent to BMM	50.01	[48.08;52.01]
u 2	Health utility of infected (min=0, max=1)	Beta(1469.3, 629.7)	Equivalent to BMM	0.70	[0.68;0.72]
u 3	Health utility of asymptomatic (min=0, max=1)	Beta(1439.4, 959.6)	Equivalent to BMM	0.60	[0.58;0.62]
u 4	Health utility of morbid (min=0, max=1)	Beta(629.7, 1469.3)	Equivalent to BMM	0.30	[0.28;0.32]

The values are fictional and were chosen so as to produce most realistic prevalence outcome and cost-effectiveness results

Depending on the medical context, correlations between posterior distributions are possible and informative priors can result in quite different posterior distributions, if they are influenced by the updating process of other priors. This might occur in parameters which are not independent, for example transition parameters to certain states calculated by means of hierarchical models. In this respect, for example the transition probability to the state Dead might be dependent on one to a less severe state such as Morbid. Potential correlation is automatically accounted for in the PSA.

The dODE is calibrated through a frequentist probabilistic calibration approach [50, 55], whereas the BODE and BMM are calibrated through Bayesian calibration approaches [56, 57]. The first involves the calculation of goodness-of-fit statistics. The advantage of the latter is that the model parameters can be inferred by fitting the models to data directly (in one step). Simulated time series data on the number of high-risk people in the states Susceptible, Infected, Asymptomatic and Morbid in the first five years of follow-up are used for calibration. A short observation time horizon of only a couple of years (five years in this case) is common in available time series data. These data are simulated by running the dODE under the status quo for a follow-up period of five years; the input parameter values are informed through the means of the parameters of the BODE. In many infectious diseases, only data on the number of infected individuals are available. The model code can easily be adapted to only calibrate the outcome on one of the states since every state is calibrated separately. Further details are given in Appendix 2.

The dODE is estimated using a combination of the R packages EpiModel [58] and deSolve [59]. As for the BMM and BODE, we estimate the model parameters using a MCMC procedure; we run two chains with a total of 1,000 simulations after convergence. In this setting, convergence is sufficiently achieved with a Potential Scale Reduction in all model parameters [49], and there are no issues with autocorrelation. We fit the BODE in WinBUGS through the ODE solver interface WBDiff [60]. The BMM is estimated using WinBUGS and JAGS [61], an alternative, established software to perform Gibbs sampling, in order to compare computational efficiency.

The models are run on a Dell Latitude E6320 (Intel Core i5-2520M, 2x4GB DDR3 RAM, 500GB SATA HDD (2.5", 7200rpm)). The computation times are 4,480.19 seconds (around 1 hour 15 minutes) and 6,587.67 seconds (around 1 hour 50 minutes) for the dODE and BODE, respectively. Interestingly, the BMM runs much faster in JAGS (149.81 seconds) than in WinBUGS (449.42 seconds). This difference is perhaps due to the way in which the two programmes handle logical nodes, which are instrumental to defining the state allocation algorithm (see (2)). All run times include model calibration which considerably increases the computational effort in the ODE systems, but not for the BMM. The relevant model codes are presented as Additional file 1BODE.txt, BMM(WinBUGS).txt and BMM(JAGS).txt.

As discussed in “Dynamic Bayesian Markov models” section, identifiability could also become an issue when the posteriors of β and ω are multimodal, which may happen if certain parameter combinations are more likely to favour the data on partner acquisition rates and STI transmission probabilities. As a consequence, the sampling algorithm could get stuck in one of the local modes, resulting in biased summary statistics. We generate different initial values for each Markov chain to reduce the risk that parts of the posterior distribution are not visited when producing the relevant samples. In addition, we evaluate whether the two chains show proper mixing through trace plots.

Results

Natural history of disease

Figure 2 shows the outcome on the natural history of the fictional chronic STI following calibration. Only the results on high-risk females are displayed; those on high-risk males are comparable. The BMM produces results which are comparable to the “gold standard” of the BODE. The model outcome on the states Susceptible and Morbid of the BMM and BODE is basically identical, whereas the outcome of the BODE shows slightly higher estimates on the number of infected and asymptomatic high-risk females. The outcome on the number of susceptible and morbid high-risk females is higher in the dODE; in contrast, the outcome on those in the states Infected and Asymptomatic is lower when compared to the two Bayesian models. The ranges of the 95% credible intervals of the BODE and BMM are similar, showing wider ranges in the BMM. The 97.5% quantiles of the scenario analysis of the dODE are lower than the upper bounds of the 95% CIs in the states Infected and Asymptomatic and higher in the states Susceptible and Morbid, whereas the 2.5% quantiles are considerably lower than the lower bounds of the CIs (apart from the state Susceptible).

Fig. 2

Calibration results on the number of high-risk females in the states following a systematic probabilistic calibration approach. The results of the Bayesian models are similar, with a slightly higher number of high-risk females in the states Infected and Asymptomatic estimated by the Bayesian ODE-based model. In contrast, the deterministic ODE-based model results in a lower estimate on the number of high-risk females in the states Infected and Asymptomatic; however, the outcome on the state Morbid is reversed

Cost-effectiveness analysis

We denote the unit costs and utilities as c and u, with indices s, t and i representing states s∈{1,…,S}, observation time points t and interventions i=1 (status quo) and i=2 (vaccination). We assume decreasing utility values for more severe states. Costs are induced by screening, vaccination, a visit at the general practitioner and diagnostic tests. Following a positive STI diagnosis, further diagnostic tests and treatment are necessary. For all these quantities, the distributional assumptions are presented in Table 1.

The overall costs per intervention are calculated as where n are the number of people in state s at time t when intervention i is applied and δ is the yearly discount rate. In both the continuous- and discrete-time approaches, the model output on the natural history of disease infection and progression is evaluated at pre-specified time points t∈{1,…,T}, where T represents the end of follow-up. We discount both costs and benefits at a fixed yearly rate δ=0.03, following ISPOR recommendations [62]. Similarly, the overall utilities are computed as

Overall costs and utilities define the monetary net benefit NB()=kU−C. The economic evaluation is performed by calculating suitable summaries such as the increment in mean cost Δ=C2−C1 and the increment in mean effectiveness Δ=U2−U1 between vaccination and the status-quo, or the incremental cost-effectiveness ratio

In the BMM and BODE, these quantities are estimated directly as function of the parameters, while in the dODE, we conduct a scenario analysis including the 2.5% and 97.5% quantiles of the ICER to evaluate the range of “plausible” results. A cut-off point of a willingness-to-payk of approximately £20,000 – £30,000 per QALY gained, adopted by NICE [63], is used as the benchmark of value for money.

As for PSA, it is usually based on: (i) the analysis of the cost-effectiveness plane, depicting the joint probability distribution of (Δ,Δ); (ii) the cost-effectiveness acceptability curve CEAC= Pr(kΔ−Δ>0), which shows the probability that the reference intervention is cost-effective as a function of the willingness to pay k; and (iii) the expected value of “perfect” information which quantifies the maximum amount of money that the decision-maker should be willing to invest (e.g. in a new study) in order to resolve parameter uncertainty and thus make a “better” decision. The Bayesian models can perform these analyses in a straightforward way, since these quantities are all functions of the model parameters and thus a full posterior distribution can be directly obtained.

The ICER of the dODE results in £ 7,203.416, ranging between £ 2,592.44 and £ 469,906 in the scenario analysis at the 2.5% and 97.5% quantiles. The ICERs of the Bayesian models are comparable to the dODE, resulting in £ 6,054.82 and £ 6,287.62 in the BODE and BMM, respectively.

Figure 3 displays the cost-effectiveness plane. Each point of the MCMC simulation lies within the grey sustainability area, indicating that STI vaccination is cost-effective at a threshold of £ 25,000 when compared to STI screening. STI vaccination is deemed to be both more expensive and more effective than STI screening since all points are located in the upper right quadrant of the graph. The corresponding ICERs of the BMM and BODE are displayed as red and blue dots, respectively.

Fig. 3

Cost-effectiveness planes of the Bayesian ODE system and Bayesian Markov model. The cost-effectiveness plane indicates that vaccination is both more expensive and more effective than the status quo. All points lie within the sustainability area of cost-effectiveness. The ICERs of £ 6,054.82 (blue dot, BODE) and £ 6,287.62 (red dot, BMM) indicate cost-effectiveness of STI vaccination in comparison to STI screening at a threshold of £ 25,000

Figure 4 shows the results of the CEACs and the population EVPIs of the two Bayesian models. The amount of uncertainty in the BMM is slightly larger than in the BODE; however, 80% cost-effectiveness is clearly reached at the break-even points of the ICERs of both models. The population EVPI of the BMM at around £ 500,000,000 is higher than in the BODE, where it reaches a value of around £ 400,000,000. The higher EVPI value of the BMM is a consequence of the slightly larger amount of uncertainty. Thus, the value of additional research is higher in the BMM.

Fig. 4

Cost-effectiveness acceptability curves and expected value of information of the Bayesian ODE system and Bayesian Markov model. The results of the BMM are displayed in grey, whereas those of the BODE are shown in black. The amount of parameter uncertainty is higher in the BMM. The CEACs in the left panel reach values of 80% at a willingness-to-pay corresponding to the ICERs. The EVPIs for the whole population at around £ 500,000,000 and £ 400,000,000 in the BMM and BODE, respectively, are shown in the right panel

Discussion

In this paper we have presented a comparison of modelling methods for the economic evaluation of interventions in infectious disease. We acknowledge that ODE-based models have several advantages and consider the Bayesian ODE structure as ideal to combine transmission modelling with economic evaluation. However, including a large number of states and probability distributions on a high number of model parameters induces computational issues, which effectively acts as a barrier to the application of complex economic modelling in this area. This possibly explains why the extensive application of PSA is limited in comparison to many other disease areas, in health economics.

From the technical point of view, constructing a BMM is relatively simple and does not require the use of specialised software — in fact our analysis has been performed using R and the Gibbs samplers WinBUGS and JAGS, which are often used by statisticians and health economic modellers and thus by reviewers and advisers for health technology assessment agencies. This may again facilitate the communication of complex modelling assumptions and thus the economic assessment of complex interventions such as those based on vaccination programmes.

We of course acknowledge that the run time of the BODE could be reduced, for example using optimised code implemented in C. This is also true of the BMM, which could be coded using C directly. We note, however that a Bayesian implementation of an ODE system is more computationally intensive than the corresponding deterministic version, due to the MCMC component needed to obtain the estimate for the posterior distributions of the relevant parameters. In addition, our proposal aims at giving a general framework: it is possible that in specific cases modellers will be able to write highly optimised code that reduces the computational time. In general terms, however, our experience suggests that many modellers working in health economic evaluation would not be necessarily familiar with specialised code (e.g. C).

Moreover, it is important that the code used to provide the evaluation of the transmission model is made available to reviewers (e.g. NICE appraisal committee), among whom a language such as WinBUGS is familiar, but other, more advanced programming languages are not. Finally, the possibility of estimating the transmission model and embedding in the economic model has the potential to improve the overall process (because the modellers developing the latter would be equipped with the technical expertise to develop the former too).

In addition to language, a more crucial point with respect to computational efficacy is given by the choice of the MCMC sampler. An alternative to a Gibbs sampler would for example be Hamiltonian Monte Carlo (HMC), which is a more efficient algorithm. Random walks are avoided using HMC since the gradient information is used. As a consequence, the algorithm is able to sample more efficiently from regions of high probability [64].

Conclusion

Our proposal of a dynamic Bayesian Markov model can be seen as an effective compromise between the ideal ODE-based models including suitable distributions on all model parameters and simpler structures of MMs that fail to account for time-dependent changes in prevalence and the effects of herd immunity. While providing a sparser temporal resolution in the way in which transmission is modelled, our methodology has the advantage of accounting for parameter uncertainty. This in turn means that standard economic analysis, including PSA, can be performed in a straightforward way. In addition, Markov models are a well established tool in health economics, which may facilitate the translation of the modellers’ work to the regulators and assessors.

In the fictional example presented in this paper, our BMM performs just as well as the BODE, with seizable computational savings. Model predictions should always be calibrated, given time series or prevalence data are available. Systematic calibration approaches usually induce a considerably high computational effort. We could show that including a direct calibration approach in the Bayesian models, the BMM runs 15 times faster than the BODE.

Appendix 1

Sex- and behaviour-specific force of infection

In the case study described in “Case study” section, the force of infection presented in (3) is adjusted to additionally account for sex and behaviour. The indices v,v′= (Male, Female) indicate the respective sex and its opposite. For example, a male is represented through the index M; the index of his female mixing partner is F. The sexual behaviour group is represented through the index b= (Low, High). In the ODE models, the transition rate ρ(t) from the state Susceptible to the state Infected depends on the covariates sex and behaviour. In the BMM, the sex- and behavioural specific transition probability π is estimated by transforming ρ(t) according to (4).

For simplicity, we exclusively account for random mixing; members of the population of sex v randomly select sexual partners of the opposite sex v′. Because of the impact of the covariates, the estimation of the overall prevalence in the sexual partners of sex v′ is a weighted average of the prevalence in both behaviour groups of sex v′. We show the corresponding equations for the continuous-time approach as functions of t; for a discrete-time approach, these are similar.

The time-specific probability of selecting a partner from the high-risk group, which we indicate as , depends on the partner acquisition rates and as well as on the population sizes and . The probability of selecting a partner from the low-risk group is represented by . The corresponding equations adapted from [28] only account for two sexual behaviour groups and are thus extended for sex, to give

We estimate the sex-, behavioural- and time-specific force of infection where is the weighted average of the STI population prevalence, which is estimated as a function of the probabilities of selecting a partner of the opposite sex from one of the two sexual behaviour groups and the time-, sex- and behavioural-specific population prevalence . The number of infectious people is estimated as those in the state Infected of the respective sex and behaviour group. In line with (3), the force of infection ρ(t) is a function of the STI transmission probability per partnership β, the partner acquisition rates ω, and the population prevalence .

Appendix 2

Model calibration

The dODE is calibrated through a frequentist probabilistic calibration approach [50, 55]. Suitable distributions are assigned to the model parameters which are relevant in terms of natural history of disease. As a next step, sets of 50,000 parameter combinations are sampled from the distributions through Monte Carlo sampling. We calculate the sum of squared errors between the outputs of the model runs (for each parameter set combination) and simulated data for the first five years of follow-up in four states as

The dODE is based on a continuous-time approach; however, the corresponding output is evaluated at yearly time intervals t∈{1,…,5}. The simulated data on the number of high-risk people in state s at time t are indicated by y(t), and f(t ∣ ) is the model output on high-risk people given the input parameter set . As a final step, the set corresponding to the output which results in the least sum of squares is selected; it is displayed in Table 2.

In contrast to the dODE, the BMM and BODE are calibrated through a Bayesian calibration approach [56, 57]. As described in “Model parameters and related distributions” section, the advantage is that model parameters can be inferred by fitting the models to data directly (in one step).

In a Bayesian framework, data can be considered in several ways. We include simulated data on a selection of parameters to update the priors into the corresponding posteriors. However, this only ensures that the corresponding parameters are informed by available evidence. Despite posterior sampling, it could still be possible that the predicted outcome implied by the Bayesian models was not comparable to high-quality data obtained from large data registries. In the BODE, the output on the number of people in the states over follow-up corresponds to the solutions of the ODEs. In the BMM, the output corresponds to the solutions of the state allocation algorithm, which is given by (2). In addition to posterior sampling, the outputs of both models are constantly updated through the simulated time series data. In fact, this is an additional process of Bayesian inference, updating the model outcome through additional data.

The updating takes place by assigning Poisson distributions to the data; the event rate λ is then represented by the solutions of the ODEs and the state allocation algorithm, respectively. The resulting model outcome is already calibrated, and no further steps such as the calculation of goodness-of-fit statistics are necessary. The evaluation of the natural history of disease and the cost-effectiveness analysis including PSA can be conducted straightforwardly.

The model codes of the BODE (in WinBUGS) and BMM (WinBUGS and JAGS versions) are provided as additional files BODE.txt, BMM(WinBUGS).txt and BMM(JAGS).txt. (ZIP 8 kb)

Table 2

Point estimates of the parameters of the deterministic ODE-based model obtained through a frequentist probabilistic calibration approach

Parameter	Description	Point estimate
ω MH	Partner acquisition rate (high-risk males)	8.3515
ω ML	Partner acquisition rate (low-risk males)	2.4526
ω FH	Partner acquisition rate (high-risk females)	8.3836
ω FL	Partner acquisition rate (low-risk females)	1.6085
χ	Proliferation parameter	0.0100
β	STI transmission probability per partnership	0.1639
π 2,3	Transition parameter from state 2 to state 3	0.7957
π 3,4	Transition parameter from state 3 to state 4	0.0891
π 4,5	Transition parameter from state 4 to state 5	0.0232
π 1,5	Transition parameter from state 1 to state 5	0.0005

The parameter set with the best fit to simulated data minimises the sum of squared errors