Skip to content

Advertisement

  • Research article
  • Open Access
  • Open Peer Review

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

BMC Medical Research Methodology201818:82

https://doi.org/10.1186/s12874-018-0541-7

  • Received: 17 August 2017
  • Accepted: 12 July 2018
  • Published:
Open Peer Review reports

Abstract

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.

Keywords

  • Infectious disease
  • Herd immunity
  • Dynamic Markov model
  • Bayesian framework
  • Cost-effectiveness analysis
  • Health economic evaluation
  • Probabilistic sensitivity analysis

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 [1214].

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 [1517] and [1820], 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 [3335] 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 \(\mathcal {S}\) of mutually exclusive and exhaustive states describing disease infection and progression. We indicate the elements of \(\mathcal {S}\) 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 ϕr,s, where \(r,s\in \mathcal {S}\) 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
Fig. 1

Model structure of a hypothetical chronic sexually transmitted infection. The arrows represent the possible transitions. These are governed by the parameters ϕr,s with indices \(r,s\in \mathcal {S}\) 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 \(\mathcal {I}_{t} = [t, t + \kappa)\), 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 ρr,s(t), with \(r,s\in \mathcal {S}\) 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 n(t)=(n1(t),…,nS(t)), where ns(t) is the number of people in state s at time t. The corresponding ODE system is given by the set of equations

$$\begin{array}{*{20}l} \frac{dn_{1}(t)}{dt}&=\chi\!\left[n_{1}(t)+n_{2}(t) + n_{3}(t) + n_{4}(t)\right] - \rho_{1,2}(t)n_{1}(t)\\ &\quad-\rho_{1,5} n_{1}(t) \\ \frac{dn_{2}(t)}{dt}&=\rho_{1,2}(t)n_{1}(t)-\rho_{2,3} n_{2}(t)-\rho_{2,5} n_{2}(t) \\ \frac{dn_{3}(t)}{dt}&=\rho_{2,3} n_{2}(t)-\rho_{3,4} n_{3}(t)-\rho_{3,5} n_{3}(t) \\ \frac{dn_{4}(t)}{dt}&=\rho_{3,4} n_{3}(t)-\rho_{4,5} n_{4}(t) \\ \frac{dn_{5}(t)}{dt}&=\rho_{1,5} n_{1}(t)+\rho_{2,5} n_{2}(t)+\rho_{3,5} n_{3}(t) + \rho_{4,5} n_{4}(t). \end{array} $$
(1)

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 Π=(πr,s), 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
$$ \boldsymbol\Pi = \left(\begin{array}{ccccc} \pi_{1,1}& \pi_{1,2} & 0 & 0& \pi_{1,5}\\ 0& \pi_{2,2}& \pi_{2,3} & 0& \pi_{2,5}\\ 0& 0 & \pi_{3,3} & \pi_{3,4}& \pi_{3,5}\\ 0& 0 & 0 & \pi_{4,4} & \pi_{4,5}\\ 0& 0 & 0 & 0& 1\\ \end{array} \right), $$
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 nt=(n1t,…,nSt), where nst is the number of people in state s and at each time interval \(\mathcal {I}_{t}\), then transitions across the states from one time interval to the next are calculated as
$$ \boldsymbol{n}_{t+1}=\boldsymbol\Pi\boldsymbol{n}_{t}. $$
(2)

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
    $$\psi_{t} = \frac{I_{t}}{N_{t}}, $$
    where It represents the number of people in the state of infection and, assuming that state S indicates death,
    $$N_{t}=\sum_{s=1}^{S-1}n_{st}$$
    is the number of those alive at time interval \(\mathcal {I}_{t}\).
The force of infection is recalculated at each Markov cycle as
$$ \lambda_{t}=\beta \omega \psi_{t}. $$
(3)
Since ω is a rate, (3) also results in a transition rate. Assuming that λt remains constant within each time interval, the corresponding time-dependent transition probability for the discrete-time MM is estimated as
$$ \pi_{\text{{1,2,t}}} = 1-\exp^{-\lambda_{t}}. $$
(4)

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 \(\mathcal {I}_{t}\). 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
$$ \pi_{\text{{1,2,t}}}=\exp^{-\rho_{\text{{1,5}(t)}}}\left(1-\exp^{-\lambda_{t}}\right), $$
(5)

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 ψt 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 ψt and a fortiori also on λt, which is defined as a function of the three random parameters (β,ω,ψt). 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) [5254] 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 ωvb, 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 ωvb are available from a large data registry. To infer ωvb, 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 θ={ωvb,χ,β,πr,s,τ,ξ,γ,σ,c,us}. 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

ω M H

Partner acquisition rate (high-risk males)

Poisson-Gamma model

Equivalent to BMM

9.10

[8.77;9.29]

ω M L

Partner acquisition rate (low-risk males)

Poisson-Gamma model

Equivalent to BMM

2.98

[2.82;3.12]

ω F H

Partner acquisition rate (high-risk females)

Poisson-Gamma model

Equivalent to BMM

9.00

[8.71;9.26]

ω F L

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 \(\hat {R}<1.1\) 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
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 csti and us, 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
$$C_{i}=\sum\limits_{t=1}^{T}\sum\limits_{s=1}^{S}\frac{c_{si}n_{sti}}{(1+\delta)^{t-1}},$$
where nsti 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
$$U_{i}=\sum\limits_{t=1}^{T}\sum\limits_{s=1}^{S}\frac{u_{s}n_{sti}}{(1+\delta)^{t-1}}.$$
Overall costs and utilities define the monetary net benefit NBi(θ)=kUiCi. The economic evaluation is performed by calculating suitable summaries such as the increment in mean cost Δc=C2C1 and the increment in mean effectiveness Δe=U2U1 between vaccination and the status-quo, or the incremental cost-effectiveness ratio
$$\text{ICER}=\frac{\mathrm{E}[\Delta_{c}]}{\mathrm{E}[\Delta_{e}]}.$$

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 (Δe,Δc); (ii) the cost-effectiveness acceptability curve CEAC= Pr(kΔeΔc>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
$$\text{EVPI} = \mathrm{E}_{\boldsymbol{\theta}}\left[\max_{i}\text{NB}_{i}\left(\boldsymbol{\theta}\right)\right] - \max_{i}\mathrm{E}_{\boldsymbol{\theta}} \left[\text{NB}_{i}\left(\boldsymbol{\theta}\right)\right], $$
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
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
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 ρv,b,1,2(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 πv,b,1,2,t is estimated by transforming ρv,b,1,2(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 \(g_{v^{\prime }H}(t)\), depends on the partner acquisition rates \(\omega _{v^{\prime }H}\) and \(\omega _{v^{\prime }L}\) as well as on the population sizes \(N_{v^{\prime }H}(t)\) and \(N_{v^{\prime }L}(t)\). The probability of selecting a partner from the low-risk group is represented by \(g_{v^{\prime }L}(t)\). The corresponding equations adapted from [28] only account for two sexual behaviour groups and are thus extended for sex, to give
$$\begin{array}{@{}rcl@{}} g_{v^{\prime}H}(t)& = &\frac{\omega_{v^{\prime}H}N_{v^{\prime}H}(t)}{\omega_{v^{\prime}H}N_{v^{\prime}H}(t)+\omega_{v^{\prime}L}N_{v^{\prime}L}(t)} \\ g_{v^{\prime}L}(t)& = & 1-g_{v^{\prime}H}(t). \end{array} $$
We estimate the sex-, behavioural- and time-specific force of infection
$$ \rho_{v,b,1,2}(t)= \beta \omega_{vb}\overline{\psi}_{v^{\prime}}(t), $$
where
$$\overline{\psi}_{v^{\prime}}(t)=\left(g_{v^{\prime}H}(t)\frac{I_{v^{\prime}H}(t)}{N_{v^{\prime}H}(t)}+g_{v^{\prime}L}(t)\frac{I_{v^{\prime}L}(t)}{N_{v^{\prime}L}(t)}\right)$$
is the weighted average of the STI population prevalence, which is estimated as a function of the probabilities \(g_{v^{\prime }b}(t)\) 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 \(\frac {I_{v^{\prime }b}(t)}{N_{v^{\prime }b}(t)}\). The number of infectious people \(I_{v^{\prime }b}(t)\) is estimated as those in the state Infected of the respective sex and behaviour group. In line with (3), the force of infection ρv,b,1,2(t) is a function of the STI transmission probability per partnership β, the partner acquisition rates ωvb, and the population prevalence \(\overline {\psi }_{v^{\prime }}(t)\).

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
$$\mathcal{Q}(\boldsymbol{\theta})=\sum\limits_{t=1}^{5}\sum\limits_{s=1}^{4}\left[y_{s}(t)-f_{s}\left(t\mid\boldsymbol{\theta}\right)\right]^{2}. $$

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 ys(t), and fs(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.

Table 2

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

Parameter

Description

Point estimate

ω M H

Partner acquisition rate (high-risk males)

8.3515

ω M L

Partner acquisition rate (low-risk males)

2.4526

ω F H

Partner acquisition rate (high-risk females)

8.3836

ω F L

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

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.

Abbreviations

ANOVA: 

Analysis of variance

BCEA: 

Bayesian cost-effectiveness analysis

BMM: 

Bayesian Markov model

BODE: 

Bayesian ODE-based model

CEA: 

Cost-effectiveness analysis

CEAC: 

Cost-effectiveness acceptability curve

CI: 

Credible interval

DDR: 

Double data rate

dODE: 

deterministic ODE-based model

EVPI: 

Expected value of perfect information

EVPPI: 

Expected value of partial perfect information

EVSI: 

Expected value of sample information

FAST: 

Fourier amplitude sensitivity test

GB: 

Gigabyte

HDD: 

Hard disk drive

HIV: 

Human immunodeficiency virus

HMC: 

Hamiltonian Monte Carlo

HPV: 

Human papillomavirus

ICER: 

Incremental cost-effectiveness ratio

IPP: 

Independent inspection process

ISPOR: 

International society for pharmacoeconomics and outcomes research

JAGS: 

Just another gibbs sampler

MCMC: 

Markov chain Monte Carlo

MM: 

Markov model

NB: 

Monetary net benefit

NICE: 

National institute for health and care excellence

ODE: 

Ordinary differential equation

PSA: 

Probabilistic sensitivity analysis

QALY: 

Quality-adjusted life year

RAM: 

Random access memory

SATA: 

Serial advanced technology attachment

STI: 

Sexually transmitted infection

VoI: 

Value of information

WBDiff: 

WinBUGS differential interface

WinBUGS: 

Bayesian inference using Gibbs sampling (for windows)

Declarations

Acknowledgements

Dr Simón Lunagómez, lecturer at Lancaster University, provided thought-provoking and helpful advice in writing and revising the manuscript. His area of expertise is Bayesian modelling for network data with possible application to infectious disease. In addition, he is experienced in competing risk analysis.

Funding

Dr Gianluca Baio is partially supported as the recipient of a research grant sponsored by Mapi Group at University College London. The funding body was involved in the design of the study and in simulation, analysis and interpretation of data. The study, however, is based entirely on simulated data on a fictional chronic sexually transmitted infection.

Authors’ contributions

KH conducted the literature review on methodology in dynamic infectious disease transmission modelling. KH and GB contributed the approach of a dynamic Bayesian Markov model. KH wrote the corresponding programmes in JAGS and WinBUGS. GB and AH provided helpful advice in selecting the calibration approaches as well as improving the methodology and the writing. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Authors’ Affiliations

(1)
University College London, Department of Statistical Science, Torrington Place, London, WC1E 7JE, UK
(2)
ICON plc Clinical Research Organisation, Konrad-Zuse-Platz 11, München, 81829, Germany

References

  1. Orenstein W. Eradicating polio: how the world’s pediatricians can help stop this crippling illness forever. Pediatrics. 2015; 135(1):196–202.View ArticlePubMedGoogle Scholar
  2. Witty C. Milroy Lecture: eradication of disease: hype, hope and reality. Clin Med. 2014; 14(4):419–21.View ArticleGoogle Scholar
  3. Weiss R, Esparza J. The prevention and eradication of smallpox: a commentary on Sloane (1755). An account of inoculation. Phil Trans R Soc B. 2014; 370:1–11.Google Scholar
  4. WHO. Immunization, Vaccines and Biologicals. 2013. Available from http://www.who.int/immunization/monitoring/_surveillance/data/en/. Accessed 02 Jan 2017.
  5. WHO. Making Choices in Health: WHO Guide to Cost-Effectiveness Analysis. 2003. Available from: http://www.who.int/choice/publications/p/_2003/_generalised/_cea.pdf. Accessed 02 Jan 2017.
  6. Briggs A, Sculpher M, Claxton K. Decision Modelling for Health Economic Evaluation. Oxford: Oxford University Press; 2006.Google Scholar
  7. National Institute for Health and Care Excellence (NICE). Guide to the methods of technology appraisal 2013. London: NICE. 2013. Available from: http://www.nice.org.uk/article/pmg9. Accessed 02 Jan 2017.
  8. Baio G, Dawid A. Probabilistic Sensitivity Analysis in Health Economics. Stat Methods Med Res. 2011; 0:1–20.Google Scholar
  9. Baio G. Bayesian Methods in Health Economics. Boca Raton: Chapman & Hall, CRC Biostatistics Series; 2013.Google Scholar
  10. JCVI. Joint Committee on Vaccination and Immunisation. Code of Practice. 2013. Available from: https://www.gov.uk/government/uploads/system/uploads/attachment/_data/file/224864/JCVI/_Code/_of/_Practice/_revision/_2013/_-/_final.pdf. Accessed 02 Jan 2017.Google Scholar
  11. Anderson R, May R. Vaccination and herd immunity to infectious diseases. Nature. 1985; 318:323–29.View ArticlePubMedGoogle Scholar
  12. Pitman R, Fisman D, Zaric G, Postma M, Kretzschmar M, Edmunds J, Brisson M. Dynamic Transmission Modeling: A Report of the ISPOR-SMDM Modeling Good Research Practices Task Force-5. Value Health. 2012; 15:828–34.View ArticlePubMedGoogle Scholar
  13. Edmunds J, Medley G, Nokes D. Evaluating The Cost-Effectiveness Of Vaccination Programmes: A Dynamic Perspective. Stat Med. 1999; 18:3263–82.View ArticlePubMedGoogle Scholar
  14. Chong K, Zee B, Wang M. A statistical method utilizing information of imported cases to estimate the transmissibility for an influenza pandemic. BMC Med Res Methodol. 2017; 17(31):1–9.Google Scholar
  15. Blower S, Dowlatabadi H. Sensitivity and Uncertainty Analysis of Complex Models of Disease Transmission: An HIV Model, as an Example. Int Stat Rev. 1994; 62:229–43.View ArticleGoogle Scholar
  16. Jit M, Choi Y, Edmunds W. Economic Evaluation of Human Papillomavirus Vaccination in the United Kingdom. BMJ. 2008; 337:1–12.View ArticleGoogle Scholar
  17. Jit M, Chapman R, Hughes O, Choi Y. Comparing Bivalent and Quadrivalent Human Papillomavirus Vaccines: Economic Evaluation Based on Transmission Model. BMJ. 2011;:343:d5775.Google Scholar
  18. Khazeni N, Hutton D, Garber A, Hupert N, Owens D. Effectiveness and Cost-Effectiveness of Vaccination against Pandemic (H1N1) 2009. Ann Intern Med. 2009; 151:829–39.View ArticlePubMedPubMed CentralGoogle Scholar
  19. Alistar S, Owens D, Brandeau M. Effectiveness and Cost Effectiveness of Expanding Harm Reduction and Antiretroviral Therapy in a Mixed HIV Epidemic: A Modeling Analysis for Ukraine. PLoS Med. 2011; 8:1–15.View ArticleGoogle Scholar
  20. Juusola J, Brandeau M, Long E, Owens D, Bendavid E. The cost-effectiveness of symptom-based testing and routine screening for acute HIV infection in men who have sex with men in the USA. AIDS. 2011; 25:1779–87.View ArticlePubMedPubMed CentralGoogle Scholar
  21. Xiu D, Karniadakis G. The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations. SIAM J Sci Comput. 2002; 24:619–44.View ArticleGoogle Scholar
  22. Sudret B. Global sensitivity analysis using polynomial chaos expansions. Reliab Eng Syst Saf. 2008; 93:964–79.View ArticleGoogle Scholar
  23. Saltelli A, Sobol I. About the use of rank transformation in sensitivity analysis of model output. Reliab Eng Syst Saf. 1995; 50:225–39.View ArticleGoogle Scholar
  24. Saltelli A, Tarantola S, Chan K. A Quantitative Model-Independent Method for Global Sensitivity Analysis of Model Output. Technometrics. 1999; 41:39–56.View ArticleGoogle Scholar
  25. Zaric G, Bayoumi A, Brandeau M, Owens D. The Cost Effectiveness of Counseling Strategies to Improve Adherence to Highly Active Antiretroviral Therapy (HAART) Among Men Who Have Sex With Men. Med Dec Making. 2008; 28(3):359–76.View ArticleGoogle Scholar
  26. Long E, Stavert R. Portfolios of Biomedical HIV Interventions in South Africa: A Cost-Effectiveness Analysis. J Gen Intern Med. 2013; 28(10):1294–301.View ArticlePubMedPubMed CentralGoogle Scholar
  27. Andronis L, Barton P, Bryan S. Sensitivity analysis in economic evaluation: an audit of NICE current practice and a review of its use and value in decision-making. Health Technol Assess. 2009; 13(29):1–100.View ArticleGoogle Scholar
  28. Vynnycky E, White R. An Introduction to Infectious Disease Modelling. New York: Oxford University Press; 2010.Google Scholar
  29. Brisson M, Edmunds W. Economic Evaluation of Vaccination Programs: The Impact of Herd-Immunity. Med Dec Making. 2003; 23:76–82.View ArticleGoogle Scholar
  30. Ross E, Cinti S, Hutton D. A Cost-Effective, Clinically Actionable Strategy for Targeting HIV Preexposure Prophylaxis to High-Risk Men Who Have Sex With Men. J Acquir Immune Defic Syndr. 2016; 72:61–7.View ArticleGoogle Scholar
  31. Cooper B, Lipsitch M. The analysis of hospital infection data using hidden Markov models. Biostatistics. 2004; 5(2):223–37.View ArticlePubMedGoogle Scholar
  32. Gibson GJ, Renshaw E. Likelihood estimation for stochastic compartmental models using Markov chain methods. Stat Comput. 2001; 11(4):347–58.View ArticleGoogle Scholar
  33. Forrester M, Pettitt A. Use of Stochastic Epidemic Modeling to Quantify Transmission Rates of Colonization With Methicillin-Resistant Staphylococcus aureus in an Intensive Care Unit. Infect Control Hosp Epidemiol. 2005; 26(7):598–606.View ArticlePubMedGoogle Scholar
  34. Auranen K, Eichner M, Käyhty H, Takala A, Arjas E. A Hierarchical Bayesian Model to Predict the Duration of Immunity to Haemophilus influenzae Type b. Biometrics. 1999; 55(4):1306–13.View ArticlePubMedGoogle Scholar
  35. Gibson GJ, Austin EJ. Fitting and testing spatio-temporal stochastic models with application in plant epidemiology. Plant Pathol. 1996; 45(2):172–84.View ArticleGoogle Scholar
  36. Haeussler K, Marcellusi A, Mennini F, Favato G, Picardo M, Garganese G, Bononi M, Costa S, Scambia G, Zweifel P, Capone A, Baio G. Cost-Effectiveness Analysis of Universal Human Papillomavirus Vaccination Using a Dynamic Bayesian Methodology: The BEST II Study. Value Health. 2015; 18(8):956–68.View ArticlePubMedGoogle Scholar
  37. Spiegelhalter D, Best N. Bayesian Approaches to Multiple Sources of Evidence and Uncertainty in Complex Cost-Effectiveness Modelling. Stat Med. 2003; 22:3687–709.View ArticlePubMedGoogle Scholar
  38. Alizon S, Magnus C. Modelling the Course of an HIV Infection: Insights from Ecology and Evolution. Viruses. 2012; 4:1984–2013.View ArticlePubMedPubMed CentralGoogle Scholar
  39. Welton N, Sutton A, Cooper N, Abrams K, Ades A. Evidence Synthesis for Decision Making in Healthcare. Chichester: John Wiley & Sons, Ltd.; 2012.View ArticleGoogle Scholar
  40. Unit MB. WinBUGS 14. 2016. Available from: http://www.mrc-bsu.cam.ac.uk/software/bugs/the-bugs-project-winbugs/. Accessed 02 Jan 2017.
  41. Stan Modeling Language. User’s Guide and Reference Manual, Version 2.14.0. 2016. Available from: http://mc-stan.org/. Accessed 02 Jan 2017.Google Scholar
  42. Bilcke J, Chapman R, Atchison C, Cromer D, Johnson H, Willem L, Cox M, Edmunds W, Jit M. Quantifying Parameter and Structural Uncertainty of Dynamic Disease Transmission Models Using MCMC: An Application to Rotavirus Vaccination in England and Wales. Med Dec Making. 2015; 35:633–47.View ArticleGoogle Scholar
  43. van Rosmalen J, Toy M, O’Mahony J. A Mathematical Approach for Evaluating Markov Models in Continuous Time without Discrete-Event Simulation. Med Dec Making. 2013; 33:767–79.View ArticleGoogle Scholar
  44. Hudgens M, Li C, Fine J. Parametric Likelihood Inference for Interval Censored Competing Risks Data. Biometrics. 2014; 70:1–9.View ArticlePubMedPubMed CentralGoogle Scholar
  45. Andersen P, Abildstrom S, Rosthoj S. Competing risks as a multi-state model. Stat Methods Med Res. 2002; 11:203–15.View ArticlePubMedGoogle Scholar
  46. Dreyer T, Van Vuuren J. A comparison between continuous and discrete modelling of cables with bending stiffness. Appl Math Model. 1999; 23(7):527–41.View ArticleGoogle Scholar
  47. Cooper N, Abrams K, Sutton A, Turner D, Lambert P. A Bayesian approach to Markov modelling in cost-effectiveness analyses: application to taxane use in advanced breast cancer. J R Stat Soc Ser A. 2003; 166:389–405.View ArticleGoogle Scholar
  48. Korostil I, Peters G, Cornebise J, Regan D. Adaptive Markov Chain Monte Carlo Forward Projection for Statistical Analysis in Epidemic Modelling of Human Papillomavirus. Stat Med. 2012; 32(11):1917–53.View ArticlePubMedGoogle Scholar
  49. Gelman A, Carlin J, Stern H, Rubin D. Bayesian Data Analysis. London: Chapman & Hall/CRC; 2004.Google Scholar
  50. Vanni T, Karnon J, Madan J, White R, Edmunds W, Foss A, Legood R. Calibrating models in economic evaluation: a seven-step approach. Pharmacoeconomics. 2011; 29:35–49.View ArticlePubMedGoogle Scholar
  51. Package BCEA. Bayesian cost-effectiveness analysis. 2015. Available from: http://cran.r-project.org/web/packages/BCEA/BCEA.pdf. Accessed 02 Jan 2017.Google Scholar
  52. Heath A, Manolopoulou I, Baio G. Estimating Multiparameter Partial Expected Value of Perfect Information from a Probabilistic Sensitivity Analysis Sample: A Nonparametric Regression Approach. Med Dec Making. 2014; 34(3):311–26.View ArticleGoogle Scholar
  53. Heath A, Manolopoulou I, Baio G. Estimating the expected value of partial perfect information in health economic evaluations using integrated nested Laplace approximation. Stat Med. 2016; 35(23):4264–80.View ArticlePubMedPubMed CentralGoogle Scholar
  54. Heath A, Manolopoulou I, Baio G. Efficient Monte Carlo Estimation of the Expected Value of Sample Information Using Moment Matching. Med Dec Making. 2018; 38(2):163–73.View ArticleGoogle Scholar
  55. Van de Velde N. Modeling Human Papillomavirus Vaccine Effectiveness: Quantifying the Impact of Parameter Uncertainty. Am J Epidemiol. 2007; 165:762–75.View ArticlePubMedGoogle Scholar
  56. de Angelis D, Sweeting M, Ades A, Hickman M, Hope V, Ramsay M. An evidence synthesis approach to estimating Hepatitis C Prevalence in England and Wales. Stat Methods Med Res. 2009; 19:361–79.View ArticleGoogle Scholar
  57. Welton N, Ades A. A model of toxoplasmosis incidence in the UK: evidence synthesis and consistency of evidence. J R Stat Soc: Ser C: Appl Stat. 2005; 54:385–404.View ArticleGoogle Scholar
  58. Jenness S. Package EpiModel. 2015. Available from: https://cran.r-project.org/web/packages/EpiModel/EpiModel.pdf. Accessed 02 Jan 2017.Google Scholar
  59. Soetaert K, Petzoldt T, Woodrow Setzer R. Package deSolve: Solving Initial Value Differential Equations in R. 2015. Available from: https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf. Accessed 02 Jan 2017.Google Scholar
  60. Lunn D. WinBUGS Differential Interface (WBDiff). 2004. Available from: http://winbugs-development.mrc-bsu.cam.ac.uk/wbdiff.html. Accessed 02 Jan 2017.Google Scholar
  61. Plummer M. JAGS Version 4.0.0 user manual. 2015. Available from: http://sourceforge.net/projects/mcmc-jags/files/Manuals/4.x/jags/_user/_manual.pdf. Accessed 02 Jan 2017.Google Scholar
  62. ISPOR. ISPOR Pharmacoeconomics guidelines. 2012. Available from: http://www.ispor.org/peguidelines/index.asp. Accessed 02 Jan 2017.
  63. Rawlins M, Culyer A. National Institute for Clinical Excellence and its Value Judgements. BMJ. 2004; 329:224–27.View ArticlePubMedPubMed CentralGoogle Scholar
  64. Brooks S, Gelman A, Jones G, Meng X-L. Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC, Boca Raton: Academic Press; 2011.View ArticleGoogle Scholar

Copyright

Advertisement