 Research article
 Open Access
 Open Peer Review
 Published:
A dynamic Bayesian Markov model for health economic evaluations of interventions in infectious disease
BMC Medical Research Methodology volume 18, Article number: 82 (2018)
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 ODEbased 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 ODEbased 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 prerequisite to their implementation, health interventions such as vaccines are thus increasingly subject to costeffectiveness 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 costeffectiveness [10]. However, there are currently no vaccinespecific guidelines for developing clinical or costeffectiveness 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 decisionmaking 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 ODEbased models, systems of equations can also be defined in discretetime, 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 discretetime 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 resampling 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 timedependent 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 HPVinduced 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 runtime 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 ODEbased 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 prespecified 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.
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 socalled “continuoustime approach”. Alternatively, it is possible to assume that transitions occur in discrete time where only one transition is possible within a predefined 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 plugin 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)=(n_{1}(t),…,n_{S}(t))^{′}, where n_{s}(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 ODEbased 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 ODEbased models which incorporate distributions on all parameters are rare exceptions in the literature on infectious disease transmission modelling [42].
Discretetime 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 timespecific 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 ODEbased models, MMs can be implemented for continuous time.
However, the vast majority of MMs in the health economic literature is based on a discretetime 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
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 n_{t}=(n_{1t},…,n_{St})^{′}, where n_{st} 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
MMs are relatively straightforward to implement and are commonly used to model the progression of noncommunicable conditions such as cardiovascular disease and cancer. Therefore, they are established in the health economic literature and wellknown 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 ODEbased 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 timedependent pathogen prevalence
$$\psi_{t} = \frac{I_{t}}{N_{t}}, $$where I_{t} represents the number of people in the state of infection and, assuming that state S indicates death,
$$N_{t}=\sum_{s=1}^{S1}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
Since ω is a rate, (3) also results in a transition rate. Assuming that λ_{t} remains constant within each time interval, the corresponding timedependent transition probability for the discretetime 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 \(\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 followup, 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 followup.
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 continuoustime 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, timedependent 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 ODEbased model would be extremely high due to i) numerical integration, e.g. through the RungeKutta solver [46]; ii) considerably smaller step sizes of ODE solvers when compared to the cycle length of discretetime 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 agespecific 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 nonidentifiability, 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 followup. 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 costeffectiveness 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 decisionmaking process. This can be postprocessed (e.g. using the R package BCEA [51]) to produce relevant summaries such as the costeffectiveness plane, the costeffectiveness acceptability curve and the analysis of the value of information (see “Costeffectiveness 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 lowrisk sexual behaviour. The duration of the followup 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 highrisk 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 sexspecific 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 healthcare interventions. We assume that in the statusquo, screening takes place at intervals of five years at a predefined 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 behaviouralspecific 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 costeffectiveness 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 minimallyinformative priors to the remaining parameters and update these through simulated individuallevel 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 PoissonGamma models. BetaBinomial 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,u_{s}}. 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 discretetime approach, the corresponding transition probabilities are modelled using Beta distributions. In contrast, the transition rates of the continuoustime 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.
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 goodnessoffit 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 highrisk people in the states Susceptible, Infected, Asymptomatic and Morbid in the first five years of followup 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 followup 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 i52520M, 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 highrisk females are displayed; those on highrisk 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 highrisk females. The outcome on the number of susceptible and morbid highrisk 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).
Costeffectiveness analysis
We denote the unit costs and utilities as c_{sti} and u_{s}, 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_{sti} 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 discretetime approaches, the model output on the natural history of disease infection and progression is evaluated at prespecified time points t∈{1,…,T}, where T represents the end of followup. 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_{i}(θ)=kU_{i}−C_{i}. The economic evaluation is performed by calculating suitable summaries such as the increment in mean cost Δ_{c}=C_{2}−C_{1} and the increment in mean effectiveness Δ_{e}=U_{2}−U_{1} between vaccination and the statusquo, or the incremental costeffectiveness 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 cutoff point of a willingnesstopayk 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 costeffectiveness plane, depicting the joint probability distribution of (Δ_{e},Δ_{c}); (ii) the costeffectiveness acceptability curve CEAC= Pr(kΔ_{e}−Δ_{c}>0), which shows the probability that the reference intervention is costeffective 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 decisionmaker 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 costeffectiveness plane. Each point of the MCMC simulation lies within the grey sustainability area, indicating that STI vaccination is costeffective 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.
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% costeffectiveness is clearly reached at the breakeven 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.
Discussion
In this paper we have presented a comparison of modelling methods for the economic evaluation of interventions in infectious disease. We acknowledge that ODEbased 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 ODEbased models including suitable distributions on all model parameters and simpler structures of MMs that fail to account for timedependent 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 behaviourspecific 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 continuoustime approach as functions of t; for a discretetime approach, these are similar.
The timespecific probability of selecting a partner from the highrisk 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 lowrisk 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
We estimate the sex, behavioural and timespecific force of infection
where
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 behaviouralspecific 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 followup in four states as
The dODE is based on a continuoustime approach; however, the corresponding output is evaluated at yearly time intervals t∈{1,…,5}. The simulated data on the number of highrisk people in state s at time t are indicated by y_{s}(t), and f_{s}(t ∣ θ) is the model output on highrisk 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 highquality data obtained from large data registries. In the BODE, the output on the number of people in the states over followup 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 goodnessoffit statistics are necessary. The evaluation of the natural history of disease and the costeffectiveness analysis including PSA can be conducted straightforwardly.
Abbreviations
 ANOVA:

Analysis of variance
 BCEA:

Bayesian costeffectiveness analysis
 BMM:

Bayesian Markov model
 BODE:

Bayesian ODEbased model
 CEA:

Costeffectiveness analysis
 CEAC:

Costeffectiveness acceptability curve
 CI:

Credible interval
 DDR:

Double data rate
 dODE:

deterministic ODEbased 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 costeffectiveness 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:

Qualityadjusted 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)
References
 1
Orenstein W. Eradicating polio: how the world’s pediatricians can help stop this crippling illness forever. Pediatrics. 2015; 135(1):196–202.
 2
Witty C. Milroy Lecture: eradication of disease: hype, hope and reality. Clin Med. 2014; 14(4):419–21.
 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.
 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 CostEffectiveness 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.
 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.
 9
Baio G. Bayesian Methods in Health Economics. Boca Raton: Chapman & Hall, CRC Biostatistics Series; 2013.
 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.
 11
Anderson R, May R. Vaccination and herd immunity to infectious diseases. Nature. 1985; 318:323–29.
 12
Pitman R, Fisman D, Zaric G, Postma M, Kretzschmar M, Edmunds J, Brisson M. Dynamic Transmission Modeling: A Report of the ISPORSMDM Modeling Good Research Practices Task Force5. Value Health. 2012; 15:828–34.
 13
Edmunds J, Medley G, Nokes D. Evaluating The CostEffectiveness Of Vaccination Programmes: A Dynamic Perspective. Stat Med. 1999; 18:3263–82.
 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.
 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.
 16
Jit M, Choi Y, Edmunds W. Economic Evaluation of Human Papillomavirus Vaccination in the United Kingdom. BMJ. 2008; 337:1–12.
 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.
 18
Khazeni N, Hutton D, Garber A, Hupert N, Owens D. Effectiveness and CostEffectiveness of Vaccination against Pandemic (H1N1) 2009. Ann Intern Med. 2009; 151:829–39.
 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.
 20
Juusola J, Brandeau M, Long E, Owens D, Bendavid E. The costeffectiveness of symptombased testing and routine screening for acute HIV infection in men who have sex with men in the USA. AIDS. 2011; 25:1779–87.
 21
Xiu D, Karniadakis G. The WienerAskey Polynomial Chaos for Stochastic Differential Equations. SIAM J Sci Comput. 2002; 24:619–44.
 22
Sudret B. Global sensitivity analysis using polynomial chaos expansions. Reliab Eng Syst Saf. 2008; 93:964–79.
 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.
 24
Saltelli A, Tarantola S, Chan K. A Quantitative ModelIndependent Method for Global Sensitivity Analysis of Model Output. Technometrics. 1999; 41:39–56.
 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.
 26
Long E, Stavert R. Portfolios of Biomedical HIV Interventions in South Africa: A CostEffectiveness Analysis. J Gen Intern Med. 2013; 28(10):1294–301.
 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 decisionmaking. Health Technol Assess. 2009; 13(29):1–100.
 28
Vynnycky E, White R. An Introduction to Infectious Disease Modelling. New York: Oxford University Press; 2010.
 29
Brisson M, Edmunds W. Economic Evaluation of Vaccination Programs: The Impact of HerdImmunity. Med Dec Making. 2003; 23:76–82.
 30
Ross E, Cinti S, Hutton D. A CostEffective, Clinically Actionable Strategy for Targeting HIV Preexposure Prophylaxis to HighRisk Men Who Have Sex With Men. J Acquir Immune Defic Syndr. 2016; 72:61–7.
 31
Cooper B, Lipsitch M. The analysis of hospital infection data using hidden Markov models. Biostatistics. 2004; 5(2):223–37.
 32
Gibson GJ, Renshaw E. Likelihood estimation for stochastic compartmental models using Markov chain methods. Stat Comput. 2001; 11(4):347–58.
 33
Forrester M, Pettitt A. Use of Stochastic Epidemic Modeling to Quantify Transmission Rates of Colonization With MethicillinResistant Staphylococcus aureus in an Intensive Care Unit. Infect Control Hosp Epidemiol. 2005; 26(7):598–606.
 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.
 35
Gibson GJ, Austin EJ. Fitting and testing spatiotemporal stochastic models with application in plant epidemiology. Plant Pathol. 1996; 45(2):172–84.
 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. CostEffectiveness Analysis of Universal Human Papillomavirus Vaccination Using a Dynamic Bayesian Methodology: The BEST II Study. Value Health. 2015; 18(8):956–68.
 37
Spiegelhalter D, Best N. Bayesian Approaches to Multiple Sources of Evidence and Uncertainty in Complex CostEffectiveness Modelling. Stat Med. 2003; 22:3687–709.
 38
Alizon S, Magnus C. Modelling the Course of an HIV Infection: Insights from Ecology and Evolution. Viruses. 2012; 4:1984–2013.
 39
Welton N, Sutton A, Cooper N, Abrams K, Ades A. Evidence Synthesis for Decision Making in Healthcare. Chichester: John Wiley & Sons, Ltd.; 2012.
 40
Unit MB. WinBUGS 14. 2016. Available from: http://www.mrcbsu.cam.ac.uk/software/bugs/thebugsprojectwinbugs/. Accessed 02 Jan 2017.
 41
Stan Modeling Language. User’s Guide and Reference Manual, Version 2.14.0. 2016. Available from: http://mcstan.org/. Accessed 02 Jan 2017.
 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.
 43
van Rosmalen J, Toy M, O’Mahony J. A Mathematical Approach for Evaluating Markov Models in Continuous Time without DiscreteEvent Simulation. Med Dec Making. 2013; 33:767–79.
 44
Hudgens M, Li C, Fine J. Parametric Likelihood Inference for Interval Censored Competing Risks Data. Biometrics. 2014; 70:1–9.
 45
Andersen P, Abildstrom S, Rosthoj S. Competing risks as a multistate model. Stat Methods Med Res. 2002; 11:203–15.
 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.
 47
Cooper N, Abrams K, Sutton A, Turner D, Lambert P. A Bayesian approach to Markov modelling in costeffectiveness analyses: application to taxane use in advanced breast cancer. J R Stat Soc Ser A. 2003; 166:389–405.
 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.
 49
Gelman A, Carlin J, Stern H, Rubin D. Bayesian Data Analysis. London: Chapman & Hall/CRC; 2004.
 50
Vanni T, Karnon J, Madan J, White R, Edmunds W, Foss A, Legood R. Calibrating models in economic evaluation: a sevenstep approach. Pharmacoeconomics. 2011; 29:35–49.
 51
Package BCEA. Bayesian costeffectiveness analysis. 2015. Available from: http://cran.rproject.org/web/packages/BCEA/BCEA.pdf. Accessed 02 Jan 2017.
 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.
 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.
 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.
 55
Van de Velde N. Modeling Human Papillomavirus Vaccine Effectiveness: Quantifying the Impact of Parameter Uncertainty. Am J Epidemiol. 2007; 165:762–75.
 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.
 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.
 58
Jenness S. Package EpiModel. 2015. Available from: https://cran.rproject.org/web/packages/EpiModel/EpiModel.pdf. Accessed 02 Jan 2017.
 59
Soetaert K, Petzoldt T, Woodrow Setzer R. Package deSolve: Solving Initial Value Differential Equations in R. 2015. Available from: https://cran.rproject.org/web/packages/deSolve/vignettes/deSolve.pdf. Accessed 02 Jan 2017.
 60
Lunn D. WinBUGS Differential Interface (WBDiff). 2004. Available from: http://winbugsdevelopment.mrcbsu.cam.ac.uk/wbdiff.html. Accessed 02 Jan 2017.
 61
Plummer M. JAGS Version 4.0.0 user manual. 2015. Available from: http://sourceforge.net/projects/mcmcjags/files/Manuals/4.x/jags/_user/_manual.pdf. Accessed 02 Jan 2017.
 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.
 64
Brooks S, Gelman A, Jones G, Meng XL. Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC, Boca Raton: Academic Press; 2011.
Acknowledgements
Dr Simón Lunagómez, lecturer at Lancaster University, provided thoughtprovoking 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.
Author information
Affiliations
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.
Corresponding author
Correspondence to Katrin Haeussler.
Ethics declarations
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.
Additional file
Additional file 1
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)
Rights and permissions
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.
About this article
Received
Accepted
Published
DOI
Keywords
 Infectious disease
 Herd immunity
 Dynamic Markov model
 Bayesian framework
 Costeffectiveness analysis
 Health economic evaluation
 Probabilistic sensitivity analysis