Estimating Expected Value of Sample Information for Incomplete Data Models Using Bayesian Approximation

Abstract

Expected value of sample information (EVSI) involves simulating data collection, Bayesian updating, and reexamining decisions. Bayesian updating in incomplete data models typically requires Markov chain Monte Carlo (MCMC). This article describes a revision to a form of Bayesian Laplace approximation for EVSI computation to support decisions in incomplete data models. The authors develop the approximation, setting out the mathematics for the likelihood and log posterior density function, which are necessary for the method. They compare the accuracy of EVSI estimates in a case study cost-effectiveness model using first- and second-order versions of their approximation formula and traditional Monte Carlo. Computational efficiency gains depend on the complexity of the net benefit functions, the number of inner-level Monte Carlo samples used, and the requirement or otherwise for MCMC methods to produce the posterior distributions. This methodology provides a new and valuable approach for EVSI computation in health economic decision models and potential wider benefits in many fields requiring Bayesian approximation.

Keywords

Bayesian decision theory value of information uncertainty sample size clinical trial design stochastic CEA approximate Bayesian inference incomplete data models

Introduction

Expected Value of Sample Information

Expected value of sample information (EVSI) quantifies the expected value to the decision maker of obtaining sample information before making a decision.¹ In health economics, the use of value-of-information methods in general is an active area of methodological development, with authors investigating and promoting their use in sensitivity analysis and to quantify the value of research.^2–11 EVSI is promoted for determining optimum sample sizes and allocation rates in health and clinical studies.^12–18

The simplest form of medical research design is a clinical trial of 2 treatments with a single end point and equal sample size in each arm of the trial. In this article, we focus on computing EVSI for slightly more complex research designs whereby data might also be collected on quality of life (utilities), costs, and longer term effectiveness, often with different sample sizes for each different component of the data collection exercise. We describe a novel form of Bayesian Laplace approximation for EVSI computation to support decisions in incomplete data models, where a different sample size is proposed for different portions of the data. One immediate practical value of this novel approximation lies in reducing the number of computations and the time required when applied to the calculation of the EVSI from decision theory. The case study models demonstrate the feasibility of the new approach, showing that EVSI-approximated results are similar to the standard nested Monte Carlo sampling method but are achieved with significant computation time reductions.

For the mathematical description of EVSI, we assume a decision model with uncertain parameters θ, a joint prior probability distribution p(θ), and a choice between intervention D = {D₁, D₂, …, D_N} with net benefit functions NB(D_i, θ). A new research study would provide data $X_{θ_{I}}$ on parameters of interest θ_I, giving a posterior density via Bayesian updating $p (θ {| X}_{θ_{I}})$ and enabling a revised decision given the data. EVSI is the difference between the expected value of a decision made after the proposed research and the expected value of a decision made with only current information^16,19—that is,

EVSI = E_{X_{θ_{I}}} [\max_{D} {E_{(θ_{I}^{C}, θ_{I} | X_{θ_{I}})} [N B (D, θ_{I}, θ_{I}^{c}]}] - \max_{D} {E_{θ} [N B (D, θ)]},

where $θ = {θ_{I}, θ_{I}^{c}}$ and E_z[f(z)] = ∫[f(z)p(z)]dz is the expectation of function f averaged over values taken by random variable Z.

Current algorithms to compute equation (1) recommend nested Monte Carlo sampling combined with Bayesian updating.^16,19 First, we use Monte Carlo to produce a sample of the parameters of interest θ_Isample. Next we use Monte Carlo to simulate a data set of a specified sample size and design $X_{θ_{I s a m p l e}}$ by sampling data, assuming that the true values of the model parameters of interest are θ_Isample. The simulated data set provides new evidence on the model parameters, and next we use Bayesian updating to compute or estimate the posterior probability distribution $p (θ {| X}_{θ_{I s a m p l e}})$ . Finally, we estimate the expected net benefit produced by a revised decision given the data. To do this, we undertake a revised probabilistic sensitivity analysis using Monte Carlo integration conditional on the data set, sampling parameters from the posterior $p (θ {| X}_{θ_{I s a m p l e}})$ , to evaluate the inner expectations for each intervention in the first term of (1). This whole process is repeated many times, and the outer Monte Carlo sampling process is used to produce different simulated data sets (we denote the number of simulated data sets K), followed for each data set by a Bayesian update and then a nested inner Monte Carlo sampling process (we denote the number of inner-loop Monte Carlo samples J) to estimate the conditional expected net benefits.

The Bayesian updating process can also be complex. To date, EVSI studies^{12,14,16,19,20} have focused on prior probability distributions and complete simulated data that are conjugate, simplifying the computation using analytic formulas for the posterior probability distribution of the model parameters. Computation remains significant because we have many sampled data sets, each of which requires Bayesian updates and Monte Carlo simulation to compute the posterior expectation. If the prior and the data are nonconjugate (Weibull proportional hazards models in Brennan and Kharroubi²⁰) or if the data set is incomplete (as is the case in this article)—that is, different sample sizes for each different component of the data collection exercise—then the dominant approach currently is to use Markov chain Monte Carlo (MCMC) methods²¹ to undertake Bayesian updating and produce simulated samples from the posterior probability distribution (e.g., an application in WinBUGS http://www.mrc-bsu.cam.ac.uk/bugs/welcome.shtml or OpenBUGS http://www.openbugs.info/w/). Using MCMC methods, this is a substantial computational expense itself and must be repeated for each simulated data set, which can result in very substantial computation times.

EVSI Using a Bayesian Approximation Methodology

In 2007, we published an article setting out a novel Bayesian approximation that, it turned out, could be used to speed up EVSI computation by replacing the Bayesian updating (second) and revised probabilistic sensitivity analysis (third) steps above with a single approximation function.²² This work, building on that of Sweeting and Kharroubi,²³ provides an approximation formula for the posterior expectation of a real valued function v(θ) of a d-dimensional vector of parameters (θ) given particular sample data X—namely,

E {v (θ) | X} ≅ v (\hat{θ}) + \sum_{i = 1}^{d} (α_{i}^{-} v (θ_{i}^{-}) + α_{i}^{+} v (θ_{i}^{+}) - v (\hat{θ})) .

The formula (2) exhibits the posterior expectation of v(θ) as a first term $v (\hat{θ})$ , plus a correction term that evaluates the function v(θ) at a series of sets of 3 points in d-dimensional space, $θ_{i}^{+}, θ_{i}^{-}$ , and $\hat{θ}$ , with weights $α_{i}^{+}$ and $α_{i}^{-}$ . Full details on the specification of these points and weights, as well as the methods used to calculate them for a given p(θ) and X, are provided in Brennan and Kharroubi.²² It is important only to note here that the calculation (undertaken by the computer programming in [R]) computes partial derivatives of the log posterior density function, inverts the information matrix, and undertakes some numerical optimization on the log-likelihood function. The approximation can therefore only be applied in circumstances where the joint posterior density is a smooth differentiable function. Note also that neither the points $θ_{i}^{+}, θ_{i}^{-}$ , and $\hat{θ}$ nor the weights $α_{i}^{+}$ and $α_{i}^{-}$ depend on the function v; they depend only on the posterior density function of θ.

The inner expectation of the first term of the equation for EVSI (formula (1)) can now be approximated using formula (2). We simply let v(θ) = NB(t, θ), and then the inner expectation of the first term of (1) is

N B (t, \hat{θ}) + \sum_{i = 1}^{d} (α_{i}^{-} N B (t, θ_{i}^{-}) + α_{i}^{+} N B (t, θ_{i}^{+}) - N B (t, \hat{θ})),

and so the Bayesian approximation formula for EVSI is

E V S I ≅ E_{X_{θ_{I}}} \max_{t} {N B (t, \hat{θ}) + \sum_{i = 1}^{d} (α_{i}^{-} N B (t, θ_{i}^{-}) + α_{i}^{+} N B (t, θ_{i}^{+}) - N B (t, \hat{θ}))} - \max_{t} E_{θ} N B (t, θ) .

Note that the number of evaluations of the net benefit function required to estimate EVSI using the approximation formula appears as 3d + 1 times the number of treatment strategies under consideration, although the last term within the $N B (t, \hat{θ})$ is the same as the first term; thus, recording this value rather than recomputing it each time, the true number of evaluations of the net benefit function required to estimate EVSI is 2d + 1 times the number of treatment strategies.

In our previous methodological development²² and a more complex case study on Weibull proportional hazards models,²⁰ the assumption has been that the proposed data collection exercise would collect the same amount of data on each of the parameters (i.e., the data matrix is rectangular with the same data being collected on all of the n individuals in the sample). In this article, we develop a refinement and apply the new methodology to a case study with “unbalanced” data, comparing the accuracy of the results and the efficiency of computation for the Brennan and Kharroubi (B&K) approximation v. the traditional 2-level approach incorporating MCMC.

Extending the Approximation Method to Unbalanced (Incomplete) Data

Balanced Case

The theoretical development for the approximation and the applications undertaken in the literature to date have assumed that the same amount of data (i.e., n samples) will be collected on each of the parameters of interest. We call this the “balanced” case. Applications have included those with multivariate normal prior and data in which conjugate distributions enable the use of analytic formulas for Bayesian updating, as well as the context of the Weibull distribution for survival data, a situation where distributions are nonconjugate and MCMC is required.²⁰ Assuming the same amount of data will be collected on each of the parameters of interest is not problematic for many real case studies such as randomized controlled trials with equal sample sizes in each arm (e.g., n = 300 for treatment 1 and n = 300 treatment 2).

However, it is obviously limiting for many more complex or “unbalanced” proposed data collection exercises (e.g., collect a subset of data on utilities for n = 100 alongside a trial of n = 300). From the perspective of economic evaluation, when synthesis of evidence from a variety of sources to inform model parameters is entirely the norm, then unbalanced designs and indeed different sources of evidence for different components of the evaluation are very common. Decision makers considering investment in further data collection might consider options including further randomized clinical trials, longer term observational data collection on effectiveness or safety, cross-sectional studies on health-related quality of life, and many more. To be generally applicable to real-world case studies, the approximation method needs to be further developed to deal with these situations. We coin the term unbalanced (see “Unbalanced Case”) because there is a different sample size for different parameters—for example, collecting n₁ samples on some parameters and n₂ on some others, producing a data matrix that we can consider either as no longer rectangular or as having missing data. The consequence of unbalanced or missing data for the approximation theory as developed to date is that it causes a problem in defining the likelihood function and hence computing θ^± and α^±.

In the development that follows, we extend the theory to the unbalanced case, illustrating the mathematical development, with the case of a prior for the model parameters p(θ) taking the multivariate normal distribution and the unbalanced data sets also having a multivariate normal distribution. The resulting extension to the approximation algorithm is not limited to multivariate normal, however, and is generally applicable to all smooth differentiable joint probability distributions.

Multivariate Normal Sample

Here we review the Bayesian updating process in the case of multivariate normal. We focus on prior probability density functions that are multivariate normal and simulated data (same amount of data [i.e., n samples] will be collected on each parameter of interest) that are conjugate, enabling us to use this conjugacy to obtain analytically the posterior density.

Let θ be a vector of d unknown parameters. Suppose that a prior probability density p(θ) for θ follows a multivariate normal such that θ ~ N_d(θ₀, Σ₀). Assume further that we are collecting n independent and identically distributed samples on θ that follow a multivariate normal such that X₁, …, X_n ~ N_d(θ, Σ), where Σ is known. We find using Bayes’ theorem that the posterior distribution of θ is

p (θ | X_{1}, ..., X_{n}) \propto \exp {- (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) / 2 - n (θ - \bar{X})^{'} Σ^{- 1} | (θ - \bar{X}) | /2} \propto \exp {- (θ - θ *)^{'} {(Σ *)}^{- 1} (θ - θ *) / 2}

where, after some algebra,

\begin{array}{l} Σ * = {(Σ_{0}^{- 1} + n Σ^{- 1})}^{- 1}, \\ θ^{*} = (Σ * (Σ_{0}^{- 1} θ_{0} + n Σ^{- 1} \bar{X}) . \end{array}

The posterior distribution of θ is therefore N_d(θ^*, Σ^*). See Appendix A for more details.

In the general case of a multivariate normal sample, both θ and Σ are unknown. The normal-inverse-Wishart distribution NIW(A, d, a, c), with A, d, a, and c known, is the natural conjugate prior for this problem, and the posterior distribution is NIW(A^*, d^*, a^*, c^*), where

\begin{array}{l} a * = {(c + n)}^{- 1} (c a + n \bar{X}), \\ c * = (c + n), \\ a * = A + n S + c n {(c + n)}^{- 1} (a - \bar{X}) (a - \bar{X})^{'} . \end{array}

Again see Appendix A for more details.

Unbalanced Case

In this section, we show the unbalanced case when there is a different sample size for different parameters producing a data matrix having missing data. We will show that the Bayesian updating process can be complex in this case, and so sequential use of Bayes’ theorem is needed to tackle this problem.

Sequential Use of Bayes’ Theorem

Let $θ = (θ_{1}, \dots, θ_{p}, θ_{p + 1}, \dots, θ_{i}, θ_{i + 1}, \dots, θ_{d}) = (θ_{i} θ_{i}^{c})$ , where θ_I is the vector of the first i parameters of θ, and $θ_{i}^{c}$ is its complement. Suppose that a prior probability density p(θ) for θ is available. Assume further that we are collecting n independent and identically distributed samples on θ_p = (θ₁, …, θ_p), n₁ independent and identically distributed samples on (θ_{p + 1}, …, θ_i) where n₁ < n (say), and no samples on (θ_{i + 1}, …, θ_d). Then, the data matrix D of interest takes the following form:

(\begin{matrix} X_{11} & X_{12} & \dots & X_{1 θ_{p}} & X_{1 θ_{p + 1}} & \dots & X_{1 i} & N A & \dots & N A \\ X_{21} & X_{22} & \dots & X_{1 θ_{p}} & X_{1 θ_{p + 1}} & \dots & X_{2 i} & N A & \dots & N A \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ X_{n_{1} 1} & X_{n_{1} 2} & \dots & X_{n_{1} θ_{p}} & X_{n_{1} θ_{p + 1}} & \dots & X_{n_{1} i} & N A & \dots & N A \\ X_{(n_{1} + 1) 1} & X_{(n_{1} + 1) 2} & \dots & X_{(n_{1} + 1) θ_{p}} & N A & \dots & N A & N A & \dots & N A \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ X_{n 1} & X_{n 2} & \dots & X_{n θ_{p}} & N A & \dots & N A & N A & \dots & N A \end{matrix}) .

For sequential use of Bayes’ theorem, the data matrix D can be divided into 2 parts in some way: D = (D₁, D₂), where D₁ contains the n₁ observations on θ_I—namely,

D_{1} = (\begin{matrix} X_{11} & X_{21} & \dots & X_{1 θ_{p}} & X_{1 θ_{p + 1}} & \dots & X_{1 i} \\ X_{21} & X_{22} & \dots & X_{1 θ_{p}} & X_{1 θ_{p + 1}} & \dots & X_{2 i} \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ & ⋱ & ⋮ \\ X_{n_{1} 1} & X_{n_{1} 2} & \dots & X_{n_{1} θ_{p}} & X_{n_{1} θ_{p + 1}} & \dots & X_{n_{1} i} \end{matrix})

and D₂ contains the n₂ = (n − n₁) observations on θ_p—namely,

D_{2} = (\begin{matrix} X_{(n_{1} + 1) 1} & X_{(n_{1} + 1) 2} & \dots & X_{(n_{1} + 1) θ_{p}} \\ ⋮ & ⋮ & ⋱ & ⋱ \\ X_{n 1} & X_{n 2} & \dots & X_{n θ_{p}} \end{matrix}) .

Now we can apply Bayes’ theorem in 2 stages to get p(θ|D), the associated posterior density of θ given D:

1. The prior density p(θ) may first be updated to a posterior density p(θ|D₁) using data D₁. Thus,

p (θ | D_{1}) = p (θ) p (D_{1} | θ) / p (D_{1}),

where p(D₁|θ) is the likelihood function for θ on data D₁, and p(D₁) is the marginal density of D₁, given as usual by

p (D_{1}) = \int p (D_{1}, θ) d θ = \int p (θ) p (D_{1} | θ) d θ .

2. Next, we can update this to the final posterior density p(θ|D₁, D₂) using data D₂. Because D₁ has now been observed, the prior density for this application is the posterior p(θ|D₁) in (8). Remember that the terms prior and posterior must always be interpreted relative to some data. p(θ|D₁) is the density of θ posterior to incorporating data D₁ but prior to D₂. Similarly, the likelihood in this application of Bayes’ theorem is p(D₂|θ, D₁). Then,

p (θ | D_{1}, D_{2}) = \int p (θ | D_{1}) p (D_{2} | θ, D_{1}) / p (D_{2} / D_{1}),

where

p (D_{2} | D_{1}) = \int p (D_{2}, θ | D_{1}) d θ = \int p (θ | D_{1}) p (D_{2} / θ, D_{1}) d θ .

Substituting for (8) in (10), we finally get the associated posterior density of θ given the full data D = (D₁, D₂) in the following form:

p (θ | D_{1}, D_{2}) = p (θ) p (D_{1} | θ) p (D_{2} | θ, D_{1}) / p (D_{1}, D_{2}),

where, by (9) and (11),

\begin{array}{l} p (D_{1}, D_{2}) = p (D_{1}) p (D_{2} | D_{1}) \\ = p (D_{1}) \int p (θ | D_{1}) p (D_{2} | θ, D_{1}) d θ \\ = \int p (θ) p (D_{1} | θ) p (D_{2} | θ, D_{1}) d θ . \end{array}

Note that to compute the posterior mode $\hat{θ}$ in (12), for the purpose of calculating all the components of (3), p(θ|D) needs only to be available up to the normalizing constant—that is,

p (θ | D) \propto p (θ) p (D_{1} | θ) p (D_{2} | θ, D_{1}),

and so the constant of proportionality (13) does not need to be computed in our analysis.

Naturally, the data may be divided further and the posterior density reached by a series of applications of Bayes’ theorem. This may be appropriate when the vector of parameters θ is divided further. Then, p(θ|D₁, D₂, …, D_n) expresses all that is known about θ after observing data D₁, D₂, …, D_n. When D_{n + 1} becomes available, Bayes’ theorem produces p(θ|D₁, D₂, … D_n, D_{n + 1}), expressing all that is now known about θ.

It is important to remember that at each stage, the likelihood is conditional on all data incorporated so far, as in the use of p(D₂|θ, D₁) in (10). It is common in many models that the data collected are independent and identically distributed. This implies that D₁ and D₂ will be independent given θ, so that p(D₂|θ, D₁) reduces to p(D₂|θ) in (10) and so in (14).

Multivariate Normal Sample (Unbalanced Case)

Here we get the associated posterior density function of θ given D in the case of multivariate normal density. As already mentioned, the prior density p(θ) may first be updated to a posterior density p(θ|D₁) using data D₁. Thus, it follows from (8) that

p (θ | D_{1}) \propto p (θ) p (D_{1} | θ),

where p(D₁|θ) is the likelihood function for θ on data D₁, which, using the properties of multivariate normal distributions in Appendix B, follows a multivariate normal N_I(θ_p, Σ_II). Thus,

p (X_{1}, \dots, X_{n_{1}} | θ) \propto \exp {- n_{1} {(θ_{I} - {\bar{X}}_{_{D_{1}}})}^{'} Σ_{II}^{- 1} (θ_{I} - {\bar{X}}_{D_{1}}) / 2} .

This implies that

p (θ | D_{1}) \propto \exp {- (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) / 2 - n_{1} {(θ_{I} - {\bar{X}}_{D_{1}})}^{'} Σ_{II}^{- 1} (θ_{I} - {\bar{X}}_{D_{1}}) / 2} .

Next, we can update this to the final posterior density p(θ|D₁, D₂) using data D₂. Because D₁ has now been observed, the prior density for this application is the posterior p(θ|D₁). The likelihood in this application of Bayes’ theorem is p(D₂|θ, D₁), the only term left to get (14). Note here that D₁ and D₂ are independent given θ, so that p(D₂|θ, D₁) reduces to p(D₂|θ). It also follows from the properties of multivariate normal distributions that the likelihood function for θ on data D₂ follows a multivariate normal N_p(θ_p, Σ_pp)—that is,

p (X_{1}, \dots, X_{n_{2}} | θ) \propto \exp {- n_{2} (θ_{p} - {\bar{X}}_{D_{2}})^{'} Σ_{pp}^{- 1} (θ_{p} - {\bar{X}}_{D_{2}}) / 2} .

This implies that the final posterior density function of θ given D is

p (θ | D_{1}) \propto \exp {- (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) / 2 - n_{1} {(θ_{I} - {\bar{X}}_{D_{1}})}^{'} Σ_{II}^{- 1} (θ_{I} I - {\bar{X}}_{D_{1}}) / 2 - n_{2} {(θ_{p} - {\bar{X}}_{D_{2}})}^{'} Σ_{pp}^{- 1} (θ_{p} - {\bar{X}}_{D_{2}}) / 2}

and so the log posterior density is readily available. Note that if n₁ = n and hence n₂ = n − n₁ = 0, the posterior distribution (4) in the balanced case is also available. Having obtained (15), we can then maximize to find the posterior mode $\hat{θ}$ and also use to compute all of the components $θ_{i}^{+}, θ_{i}^{-}, α_{i}^{+}$ , and $α_{i}^{-}$ of equation (3). This is achieved using an iterative numerical optimization process such as the Newton-Raphson technique to estimate $\hat{θ}$ . In our case study applications that follow, we used [R] software with the optim or nlm functions to minimize the negative log-posterior function. We obtain partial derivatives of the log-posterior density function with respect to each model variable analytically, although of course, if the functions are too complex to be easily differentiable, then numerical approximations will be needed for this step.²²

Case Study and Analysis

Case Study Cost-Effectiveness Model

The case study model is a simple decision tree comparing 2 strategies: treatment with drug T0 v. treatment with drug T1.²² Table 1 shows the 19 uncertain model parameters, with prior mean values shown for T0 (column a), T1 (column b), and hence the incremental analysis (column c). Costs include cost of drug and cost of hospitalizations (product of % hospitalized, days in hospital, and cost per day). Benefits (quality-adjusted life years [QALYs]) come from responders receiving a utility improvement for a specified duration, and some patients have side effects with a utility decrement for a specified duration. The illustrative threshold cost per QALY is set at λ_w = $10, 000 (this is a purely illustrative model—the value of λ_w was arbitrarily chosen and is lower than the threshold of many Western countries). Uncertain parameters have multivariate normal distributions, with standard deviations (columns d and e). There are correlations between several of the model parameters. Parameters θ⁵,θ⁷,θ¹⁴, and θ¹⁶ are each correlated (correlation coefficient = 0.6), and parameters θ⁶ and θ¹⁵ are independent of these but correlated with each other (again 0.6 correlation). Each parameter can be informed by collection of further data on individual patients, and it is assumed that patient-level variance is known for each parameter (columns f and g). The net benefit function for each treatment takes a sum-product form—that is,

Table 1

Summary of Illustrative Model

	a	b	c	d	e	f	g
	Parameter Mean Values Given			Uncertainty in Means		Patient-Level Variability
	Existing Evidence			Standard Deviations		Standard Deviations
Parameters	T0	T1	(T1–T0)	T0	T1	T0	T1
Cost of drug (θ¹, θ¹¹), $	1000	1500	500	1	1	500	500
% Admissions (θ², θ¹²)	10	8	−2	2	2	25	25
Days in hospital (θ³, θ¹³)	5.20	6.10	0.90	1.00	1.00	4.00	4.00
Cost per day (θ⁴), $	400	400	0	200	200	200	200
% Responding (θ⁵, θ¹⁴)	70	80	10	10	10	20	20
Utility change if respond (θ⁶, θ¹⁵)	0.30	0.30	0	0.10	0.05	0.20	0.20
Duration of response, y (θ⁷, θ¹⁶)	3.0	3.0	0	0.5	1.0	1.0	2.0
% Side effects (θ⁸, θ¹⁷)	25	20	−5	10	5	20	10
Change in utility if side effect (θ⁹, θ¹⁸)	−0.10	−0.10	0.00	0.02	0.02	0.10	0.10
Duration of side effect (years) (θ¹⁰, θ¹⁹)	0.50	0.50	0	0.20	0.20	0.80	0.80
Central estimate results
Net benefit given threshold (λ = $10,000), $	5057.2	5584.7	527.5
Overall expected value of perfect information per patient (estimated using 1,000,000 Monte Carlo samples), $			844.4

\begin{array}{l} N B (T 0) = λ_{w} (θ^{5} θ^{6} θ^{7} + θ^{8} θ^{9} θ^{10}) - (θ^{1} + θ^{2} θ^{3} θ^{4}), \\ N B (T 1) = λ_{w} (θ^{14} θ^{15} θ^{16} + θ^{17} θ^{18} θ^{19}) - (θ^{11} + θ^{12} θ^{13} θ^{4}) . \end{array}

There is considerable decision uncertainty in this illustrative model. The overall expected value of perfect information (EVPI) is estimated at $844.4 per patient, using 1,000,000 MC simulated samples in a probabilistic sensitivity analysis. Note when we examine the EVSI results later, we index them to the overall EVPI (i.e., index of 1.00 is equivalent to $844.4).

Analysis Plan

The objectives of the analysis are to examine the accuracy and efficiency of the B&K algorithm compared to the standard approach of using MCMC in OpenBUGS, in the context of unbalanced data collection.

To undertake unbalanced data collection exercises, we defined a separate sample size for 3 different parameter subsets:

n₁ is the sample size for a randomized controlled clinical trial measuring only short-term effectiveness response rate parameters (parameters θ⁵ and θ¹⁴).

n₂ is the sample size for an observational study on utility parameters only (parameters θ⁶ and θ¹⁵).

n₃ is the sample size for an observational study of the duration of long-term response to therapy (parameters θ⁷ and θ¹⁶).

We defined 6 unbalanced data collection exercises with values for (n₁, n₂, n₃) equal to (10, 50, 200), (10, 200, 50), (50, 10, 200), (50, 200, 10), (200, 10, 50), and (200, 50, 10). This allowed us to explore the effect of a range of sample sizes across each of the parameter subset domains.

The first stage of the analysis was to determine the required number of posterior distribution “inner-loop” samples, J, for some given level of accuracy when estimating the expected net benefit via MCMC in OpenBUGS. We considered that a reasonable level of accuracy for estimating an expected net benefit conditional on sample data would be achieved when the Monte Carlo error (i.e., the standard error of the expectation) was less than 0.5% of the expectation.

The second stage of the analysis was to compare the performance of the B&K algorithm with the MCMC-based algorithm with J inner-loop samples. For each data collection exercise, we generated K = 1000 sample data sets. Conditional on each of these, we calculated the expected net benefits of the 2 treatments using both B&K and MCMC methods. We compared the maximum posterior expected net benefit estimated by the 2 methods for each of the K samples. Finally, we compared the runtime for the B&K algorithm with the MCMC-based algorithm.

The third stage of the analysis was to determine the number of outer-loop sample data sets, K, that are required to produce a stable estimate of the EVSI for each data collection exercise. For this purpose, we used only the B&K method. We considered that a reasonable level of accuracy for estimating the EVSI would be achieved when the standard error of the EVSI estimate was less than 0.5% of the mean EVSI. Once we determined an adequate value of K, we calculated the EVSI for the 6 exercises, along with the EVSI index, defined as the ratio of the EVSI to the overall EVPI.

The final stage of the analysis was to expand the range of sample sizes on which we examine the B&K results to show EVSI estimates for all 216 (6³) possible combinations for (n₁, n₂, n₃) being n = 0, 10, 50, 100, 200, and 500. In each case, we sampled the largest data set (i.e., n1 = n2 = n3 = 500) K = 1000 times, and then, to sample smaller data sets, we discarded excess data. This enabled us to reduce the additional noise that would have resulted from resampling.

The economic model and B&K approximation method were written in R 2.11.1 (http://www.r-project.org/). The inner MCMC-based sampling was carried out in OpenBUGS 3.1.2 (http://www.openbugs.info/w/). We conducted all analyses on the Sheffield University Linux high-performance computer cluster (“Iceberg”). This allowed us to run the K outer runs of the model in parallel on separate processors with a significant gain in speed. Each processor on the cluster was clocked at 2.55 GHz, and we report all timings in either processor-seconds or processor-hours.

The models and the programs to undertake both the traditional 2-level nested MCMC and proposed B&K approximation EVSI computations are available on the CHEBS website (http://www.shef.ac.uk/chebs/software).

Case Study Results

Determining an Adequate Inner-Loop Size (J) for the MCMC Method

In our case study, the parameters in $θ_{i}^{c}$ are pairwise independent and independent of all parameters in θ_I. Because our model is of sum-product form, we were therefore able to plug in expectations for $θ_{i}^{c}$ into (1). This had the dual effect of increasing the speed of the MCMC sampling in OpenBUGS while reducing Monte Carlo sampling error.

Table 2 shows the Monte Carlo error as a proportion of the expected net benefit (for treatments T0 and T1) at a range of values of J for the 6 data collection exercises. The time taken for a single outer-loop sample K = 1 ranged from 3 processor-seconds for J = 1000 inner samples to 92 processor-seconds for J = 100,000 samples. Convergence, as assessed by the Brooks Gelman Rubin method, was achieved at around 5000 samples, so in each case, a burn-in of 10,000 samples was discarded before collecting the J inner-loop samples. The Monte Carlo error for each run was calculated within OpenBUGS using the “batch method.”²⁴ A Monte Carlo error of less than 0.5% of expected net benefit was achieved at a sample size of J = 100, 000 in all 6 data collection exercises. We therefore chose to use J = 100, 000 for the second stage of the analysis.

Comparison of B&K v. MCMC (J = 100,000)

The second stage of the analysis assessed the accuracy and efficiency of the B&K method to compute EVSI for 6 proposed data collection exercises compared to using the 2-level MCMC approach with K = 1000 outer loops and J = 100, 000 inner loops. Results are shown in Table 3. The runtime for the 6 data collection exercises is of the order of 24 processor-hours for the MCMC algorithm v. just under 1 processor-hour for the B&K algorithm. The B&K algorithm is therefore considerably more efficient, approximately 25 to 30 times faster than the MCMC-based algorithm with J = 100, 000.

Table 2

Monte Carlo Error as a Percentage of the Expected Net Benefit for Treatments T0 and T1 for Data Collection Exercises 1 to 6 for Different Values of J (K = 1)

		Data Collection Exercise
		1		2		3		4		5		6
J	Runtime, s	T0	T1	T0	T1	T0	T1	T0	T1	T0	T1	T0	T1
1,000	9.5	1.40	1.13	0.88	1.19	0.59	0.57	1.12	2.92	1.01	1.45	1.11	2.16
10,000	18.5	0.86	0.59	0.35	0.34	0.25	0.21	0.42	1.09	0.40	0.65	0.69	1.47
100,000	106.8	0.22	0.15	0.09	0.10	0.09	0.07	0.14	0.34	0.13	0.21	0.25	0.46
1,000,000	882.5	0.07	0.06	0.03	0.04	0.02	0.02	0.04	0.12	0.04	0.05	0.08	0.15

Table 3

Comparison of Speed and Agreement between the B&K and MCMC Methods for Computing EVSI for 6 Proposed Data Collection Exercises (K = 1000)

	Data Collection Exercise
	1	2	3	4	5	6
n ₁	10	10	50	50	200	200
n ₂	50	200	10	200	10	50
n ₃	200	50	200	10	50	10
EVSI estimates
B&K EVSI estimate	858.1	845.6	826.5	778.9	775.6	762.7
MCMC (J = 100,000) EVSI estimate	857.5	845.6	825.6	779.4	772.7	756.6
B&K indexed EVSI	1.016	1.001	0.979	0.922	0.918	0.903
MCMC indexed EVSI	1.015	1.001	0.978	0.923	0.915	0.896
Difference in indexed EVSI (B&K – MCMC)	0.001	0.000	0.001	−0.001	0.003	0.007
Measure of agreement between methods Σ_i[(BK_i − MCMC_i)/{(BK_i + MCMC_i)/2}]/K	0.011	−0.001	0.015	−0.009	0.045	0.096
Time (processor hours)
B&K	0.92	0.92	0.92	0.92	0.92	0.92
MCMC (J = 100,000)	26.3	24.2	26.4	24.2	26.4	26.4
Runtime ratio (B&K/MCMC), %	3.5	3.8	3.5	3.8	3.5	3.5

B&K, Brennan and Kharroubi; MCMC, Markov chain Monte Carlo; EVSI, expected value of sample information.

Figure 1 shows a comparison of the maximum expected net benefits between the MCMC and B&K algorithms for the 6 data collection exercises. All points lie on or very close to the line of equality. The difference in maximum expected net benefit between the 2 methods, expressed as a percentage of the mean value of the 2 methods, is on the order of 0.001% to 0.1%.

Determining an Adequate Outer-Loop Size (K)

For each data collection exercise, we determined the variance of the maximum expected net benefit from the results of the K = 1000 outer runs in stage 2. From this, we were able to calculate the number of runs required to achieve a standard error for the EVSI that was less than 0.5% of the mean EVSI. This was achieved at an outer-loop size of K = 100, 000 for all 6 data collection exercises. EVSI estimates calculated via the B&K algorithm with K = 100, 000 sample data sets along with the processor time required are shown in Table 4. The overall EVPI, calculated using a sample size of 1,000,000 to ensure stability, was $844.4. Values for the EVSI index (the ratio of the EVSI to the overall EVPI) are also shown in Table 4.

Table 4

EVSI Estimates and Timings in Processor Hours for the 6 Exercises Calculated via the B&K Method with K = 100,000 Sampled Data Sets (EVSI Index Is EVSI/EVPI)

	Data Collection Exercise
	1	2	3	4	5	6
n ₁	10	10	50	50	200	200
n ₂	50	200	10	200	10	50
n ₃	200	50	200	10	50	10
Time (processor hours)	86.7	89.1	84.9	91.9	89.1	89.6
EVSI (SE), $	807.7 (9.43)	802.1 (9.44)	764.0 (9.11)	749.4 (9.18)	743.6 (8.98)	755.3 (9.17)
EVSI index	0.96	0.95	0.90	0.89	0.88	0.89

B&K, Brennan and Kharroubi; EVSI, expected value of sample information, EVPI, expected value of perfect information.

Results for 216 Different Sample Size Designs

Figure 2 demonstrates that it is feasible to compute B&K results for large numbers of alternative unbalanced research designs. The EVSI results are estimated for the 216 (6³) possible sample size designs for the 3 parameter subgroups examined in our case study. These results have been run using K = 1000. Given that the EVSI estimate was calculated to be within the threshold of 0.5% of the mean using K = 100, 000, we can consider that these estimates using 100 times fewer samples should have a standard error 10 times larger (i.e., 5% of the mean EVSI).

Note that beyond a sample size of around n = 100 on any dimension (n₁, n₂, n₃), the additional value from collecting further data is very small, but each dimension also is important, so a large sample size for n₂ and n₃ but no further data collection on the first dimension (n₁ = 0) has a relatively low indexed EVSI and would still leave considerable decision uncertainty.

Discussion

In our previous work, we developed and tested a method for computing posterior expected net benefits given sample data, which enabled much more efficient computation of EVSI, but this method was only applicable to balanced data sets.^20,22 In this article, we have further developed the method by thinking through the mathematics of unbalanced data sets (where a different sample size is proposed for different portions of the data), developing a revised approach and applying this method to an illustrative case study. Our intended approach was to explore using a sequence of nested uses of our original method applied to balanced subsets of the data. However, in developing the mathematics, it has become clear that, even for unbalanced data sets, a single use of the method applied to a log posterior density function that accounts for data incompleteness will provide an estimate of the posterior expected net benefit. In general, the complexities of unbalanced or incomplete data models mean that the posterior distribution is not available analytically. Even distributions that are conjugate for balanced data sets (e.g., multivariate normal prior and data) become nonconjugate in incomplete data models; therefore, to compute posterior expected net benefit without our method requires the use of MCMC to provide a posterior sampling distribution. The case study therefore examines the accuracy and efficiency of estimating EVSI using our revised B&K method compared with using 2-level Monte Carlo sampling incorporating MCMC within the inner loop for each simulated data set. The results suggest that the MCMC method can require a very large number of inner-loop samples, and hence a very large number of evaluations of the expected net benefit function, to provide an accurate estimate of posterior expected net benefit. In contrast, our method provides an estimate using just 2d + 1 times the number of treatment strategies, where d is the number of uncertain model parameters (d = 19 in our case study, so 39 evaluations as opposed to thousands or hundreds of thousands). The accuracy of the revised B&K method in the case study appears at least as good as using very large numbers of MCMC inner iterations, and efficiency is shown by computing EVSI estimates around 30 times faster than when using J = 100,000 inner loops.

The approximation approach is valid for any net benefit function provided there is a finite number of decision strategy options. It works for analytic models whereby there is a system of assumptions and model parameters producing a large analytic formula for net benefit. It also works for stochastic models (e.g., estimated net benefit is the result of sampling, say, 10,000 patients receiving different treatments in an individual-level simulation model). The approximation method has much greater computational savings for computationally expensive net benefit functions because the functions are computed many fewer times. The approximation approach is also valid for any probability distributions used to characterize parameter uncertainty as long as they are smooth, differentiable functions. The essential requirement is that we can determine the various components of the approximation via numerical optimization techniques, and so we need to be able to differentiate the log posterior density function and invert the information matrix. This smooth differentiable criterion will cover the vast majority of characterized uncertainty in health economic models but importantly does exclude both discrete parametric distributions and empirical distributions such as nonparametric histograms.

Although this revised approximation method extends the scope for computation of EVSI within a reasonable time, there remain methodological issues for which further research would be potentially useful. First, it would be useful to examine empirically how many inner samples are needed for accurate estimation of maximum posterior expected net benefit in a wider set of case studies. Methodological work has been done on the number of inner samples needed in partial EVPI estimation,²⁵ which could be extended to EVSI. A second question concerns how accurate the estimated EVSI answer needs to be. Sometimes the analyst wants to “get a feel” for how important current parameter uncertainty is and how amenable to reduction via realistically sized data collection exercises the uncertainty might be. On other occasions, the aim is to perform quantified tradeoffs against costs of research, which requires a model of the costs of proposed data collection, including the “decision prevalence” (i.e., the number of people affected over the “lifetime” of the decision choice). The accuracy needed for EVSI estimates is therefore in part dependent on the costs of data collection and decision prevalence numbers.

There remain wider methodological issues for EVSI, including computation, not just for the EVSI of one particular study but how to go about estimating the best value from a range of research study options with limited research funds. Also, because there is no single global decision maker, the value of any research undertaken may well have wider than national impact, meaning that application of the approach in particular jurisdictions may be underestimating the value obtained globally. Most EVSI studies have taken a societal decision maker perspective, trading off research costs and societal benefits, but because data collection is undertaken by pharmaceutical or medical technology developers also, more case studies using a commercial net benefit as opposed to a societal net benefit function would be interesting. Finally, it would be useful to explore the relationship between these broad estimates of EVSI for decision problems and the examination of sequences of research studies with nested decision rules for progression to the next study that can be done using Bayesian clinical trial simulation approaches.²⁶

In conclusion, we have developed a revision to our Bayesian approximation formula for posterior expectations of real valued functions given observed data E(v(θ) |X), which enables its use for estimation of EVSI with proposed data collection exercises that are “unbalanced.” This is particularly important because the tradeoffs or synergistic value between more short-term trials and larger/longer observational studies are a common issue, especially in research in chronic diseases. This method has been applied to a relatively simple case study decision model and shown to be accurate and much more efficient in comparison to 2-level Monte Carlo estimation with nested MCMC to compute posterior expected net benefits. The method is generalizable to any net benefit function, including stochastic simulation models and any smooth mathematically defined probability distribution. Greater computation time savings are likely for more complex, computationally expensive decision models.

Footnotes

Appendix A

The treatment here in this appendix is based on material from O’Hagan and Forster.²⁷ Let θ be a vector of d unknown parameters. Suppose that a prior probability density p(θ) for θ follows a multivariate normal such that θ ~ N_d(θ₀, Σ₀). This implies that

p (θ) \propto \exp {- (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) / 2} .

Assume further that we are collecting n independent and identically distributed samples on θ that follow a multivariate normal such that X₁, …, X_n ~ N_d(θ, Σ). The likelihood function is

\begin{matrix} p (X_{1}, \dots, X_{n} | θ, Σ) \propto \prod_{i = 1}^{n} | Σ |^{- 1 / 2} \exp {- {(X_{i} - θ)}^{'} Σ^{- 1} (X_{i} - θ) / 2} \\ = | Σ |^{- n / 2} \exp {- \sum_{i = 1}^{n} {(X_{i} - θ)}^{'} Σ^{- 1} (X_{i} - θ) / 2} . \end{matrix}

Expanding the quadratic forms, we have

\begin{array}{l} \sum_{i = 1}^{n} (X_{i} - θ)^{'} Σ^{- 1} (X_{i} - θ) = n θ^{'} Σ^{- 1} θ - n {\bar{X}}^{'} Σ^{- 1} θ - n θ^{'} Σ^{- 1} \bar{X} \\ + \sum_{i = 1}^{n} X_{i}^{'} Σ^{- 1} X_{i} = n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) + r \end{array}

where $\bar{X} = n^{- 1} \sum_{i} X_{i}$ is the sample mean vector,

r = \sum_{i = 1}^{n} X_{i}^{'} Σ^{- 1} X_{i} - n {\bar{X}}^{'} Σ^{- 1} \bar{X} = \sum_{1 = 1}^{n} (X_{i} - \bar{X})^{'} Σ^{- 1} (X_{i} - \bar{X}) .

Because each term $(X_{i} - \bar{X})^{'} Σ^{- 1} (X_{i} - \bar{X})^{'}$ is a scalar, it equals the trace itself. Hence,

(X_{i} - \bar{X})^{'} Σ^{- 1} (X_{i} - \bar{X}) = t r Σ^{- 1} (X_{i} - \bar{X}) (X_{i} - \bar{X})^{'},

and so,

r = n t r Σ^{- 1} S,

where $S = n^{- 1} \sum_{i = 1}^{n} (X_{i} - \bar{X}) (X_{i} - \bar{X})^{'}$ is the sample covariance matrix. Therefore,

\begin{array}{l} p (X_{1}, \dots, X_{n} | θ, Σ) \propto | Σ |^{- n / 2} \\ \exp {- n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) / 2 - n t r Σ^{- 1} S / 2} \end{array}

Now if Σ is known, 2 terms can be dropped from the term above, leaving

p (X_{1}, \dots, X_{n} | θ) \propto \exp {- n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) / 2} .

It is very easy to work with the general family of all multivariate normal distributions. Given θ ~ N_d(θ₀, Σ₀), we find using Bayes’ theorem

p (θ | X_{1}, \dots, X_{n} | θ) \propto \exp {- (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) / 2 - n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) / 2} .

Now complete the square by

\begin{array}{l} (θ - θ_{0})^{'} Σ_{0}^{- 1} (θ - θ_{0}) - n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) = θ^{'} (Σ_{0}^{- 1} + n Σ^{- 1}) θ - (θ_{0}^{'} Σ_{0}^{- 1} + n {\bar{X}}^{'} Σ^{- 1}) θ \\ - θ^{'} (Σ_{0}^{- 1} θ_{0} + n Σ^{- 1} \bar{X}) + (θ_{0}^{'} Σ_{0}^{- 1} θ_{0} + n {\bar{X}}^{'} Σ^{- 1} \bar{X}) \\ = (θ - θ^{*})^{'} {(Σ *)}^{- 1} (θ - θ^{*}) | + R, \end{array}

where

Σ * = {(Σ_{0}^{- 1} + n Σ^{- 1})}^{- 1}, θ^{*} = Σ * (Σ_{0}^{- 1} θ_{0} + n Σ^{- 1} \bar{X}),

and R is a constant. The posterior distribution of θ is therefore N_d(θ^*, Σ^*). As is the case with univariate case, the posterior mean is a (matrix)-weighted average of the mean θ₀ and the sample mean $\bar{X}$ . The weights are respectively the prior precision matrix Σ₀⁻¹ and the analogous data precision nΣ⁻¹ (being the inverse of the covariance matrix of $\bar{X}$ ). The posterior precision matrix is the sum of Σ₀⁻ and nΣ⁻, and the posterior covariance matrix is smaller than either the prior covariance matrix Σ₀ or the data analogue nΣ⁻¹, in the sense that both Σ₀ − Σ^* and nΣ⁻¹ − Σ^* are positive definite.

In the general case of a multivariate sample with both θ and Σ unknown, the likelihood is (16). The conjugate prior distribution

\begin{array}{l} p (θ, Σ) \propto | Σ |^{- (d + p + 2) / 2} \\ \exp {- c (θ - a)^{'} Σ^{- 1} (θ - a) / 2 - (t r Σ^{- 1} A) / 2} \end{array}

is known as the normal-inverse-Wishart distribution NIW(A, d, a, c). Its hyperparameters are positive scalars d and c, a positive symmetric matrix A, and a vector a. Then combining this with (16), we have the posterior distribution

p (θ, Σ) \propto | Σ |^{- (d + p + 2 + n) / 2} \exp (Q / 2),

where

Q = t r Σ^{- 1} (A + n S) + c (θ - a)^{'} Σ^{- 1} (θ - a) + n (θ - \bar{X})^{'} Σ^{- 1} (θ - \bar{X}) .

Expanding the 2 quadratic forms in Q and collecting terms gives

Q = t r Σ^{- 1} A * + c^{*} (θ - a *)^{'} Σ^{- 1} (θ - a *) + n,

where

\begin{array}{l} a * = {(c + n)}^{- 1} (c a + n \bar{X}), c * = (c + n), \\ a * = A + n S + c n {(c + n)}^{- 1} (a - \bar{X}) (a - \bar{X})^{'} . \end{array}

Therefore, the posterior distribution is NIW(A^*, d^*, a^*, c^*), where d^* = d + n (O’Hagan and Forster²⁷).

Appendix B

Here we review the property of partitioning multivariate normal distributions that will form the basis for the results in the “Multivariate Normal Sample (Unbalanced Case)” section. Let X = (X₁, …, X_d) ~ N_d(θ, Σ). If X is partitioned into 2 sets of variables such that X₁ = (X₁, …, X_p) and X₂ = (X_{p + 1}, …, X_d), we may also partition θ and Σ as follows:

θ = [\begin{matrix} θ_{1} \\ θ_{2} \end{matrix}], Σ = [\begin{matrix} Σ_{11} & Σ_{12} \\ Σ_{21} & Σ_{22} \end{matrix}] .

Then it follows from the properties of multivariate normal distributions that X₁ will have a p-variate normal distribution (i.e., X₁ ~ N_p(θ₁, Σ₁₁)) and X₂ will have a (d − p)-variate normal distribution (i.e., X₂ ~ N_{d − p}(θ₂, Σ₂₂)).

MS is supported by a UK Medical Research Council Health Services Research/Health of the Public research training fellowship (grant number G0601721).

References

Raiffa

Schlaiffer

. Applied Statistical Decision Theory. New York: Wiley Interscience; 1967.

Claxton

Posnett

. An economic approach to clinical trial design and research priority-setting. J Health Econ. 1996;5:513–24.

Claxton

. The irrelevance of inference: a decision-making approach to the stochastic evaluation of health care technologies. J Health Econ. 1999;18:341–64.

Thompson

Graham

. Going beyond the single number: using probabilistic risk assessment to improve risk management. Hum Ecol Risk Assess. 1996;2(4):1008–34.

Felli

Hazen

. Sensitivity analysis and the expected value of perfect information. Med Decis Making. 1998;18:95-109.

Meltzer

. Addressing uncertainty in medical cost-effectiveness analysis: implications of expected utility maximization for methods to perform sensitivity analysis and the use of cost-effectiveness analysis to set priorities for medical research. J Health Econ. 2001;20:109–29.

Brennan

Kharroubi

Chilcott

O’Hagan

. A two level Monte Carlo approach to calculation expected value of sample information: how to value a research design. Presented at the 24th Annual Meeting of the Society for Medical Decision Making, October 23, 2002, Baltimore, MD.

Coyle

Buxton

O’Brien

. Measures of importance for economic analysis based on decision modeling. J Clin Epidemiol. 2003;56:989–97.

Tappenden

Chilcott

Eggington

Oakley

McCabe

. Methods for expected value of information analysis in complex health economic models: developments on the health economics of interferon-beta and glatiramer acetate for multiple sclerosis. Health Technol Assess. 2004;8(27):iii, 1–78.

10.

Yokota

Thompson

. Value of information literature analysis (VOILA): a review of applications in health risk management. Med Decis Making. 2004;24(3):287–98.

11.

Yokota

Thompson

. Value of information (VOI) analysis in environmental health risk management (EHRM). Risk Anal. 2004;24(3):635–50.

12.

Claxton

Neumann

Araki

Weinstein

. Bayesian value-of-information analysis: an application to a policy model of Alzheimer’s disease. Int J Technol Assess. 2001;17:38–55.

13.

Claxton

Thompson

. A dynamic programming approach to efficient clinical trial design. J Health Econ. 2001;20:797–822.

14.

Chilcott

Brennan

Booth

Karnon

Tappenden

. The role of modelling in prioritising and planning clinical trials. Health Technol Assess. 2003;7(23):iii, 1–125.

15.

Brennan

Kharroubi

Chilcott

O’Hagan

. A two level Monte Carlo approach to calculating expected value of perfect information: Resolution of the uncertainty in methods. Working paper, University of Sheffield, England, 2002.

16.

Ades

Claxton

. Expected value of sample information calculations in medical decision modeling. Med Decis Making. 2004;24:207–27.

17.

Eckermann

Willan

. Expected value of information and decision making in HTA. Health Econ. 2007;16:195–209.

18.

Eckermann

Willan

. The option value of delay in health technology assessment. Med Decis Making. 2008;28:300–5.

19.

Brennan

Chilcott

Kharroubi

O’Hagan

. A two level Monte Carlo approach to calculating expected value of sample information: how to value a research design. Working paper, University of Sheffield, England, 2002.

20.

Brennan

Kharroubi

. Expected value of sample information for Weibull survival data. Health Econ. 2008;16:1205–25.

21.

Spiegelhalter

Myles

Jones

Abrams

. Bayesian methods in health technology assessment: a review. Health Technol Assess. 2000;4:1–130.

22.

Brennan

Kharroubi

. Efficient computation of partial expected value of sample information using Bayesian approximation. J Health Econ. 2007;26:122–48.

23.

Sweeting

Kharroubi

. Some new formulae for posterior expectations and Bartlett corrections. Test. 2003;12:497–521.

24.

Roberts

. Markov chain concepts related to sampling algorithms. In: Gilks

Richardson

Spiegelhalter

eds. Markov Chain Monte Carlo in Practice: Interdisciplinary Statistics. London: Chapman & Hall; 1995. p 45–54.

25.

Oakley

Brennan

Tappenden

Chilcott

. Sample sizes for Monte Carlo partial EVPI calculations. Research Report No. 568/06, Department of Probability and Statistics, University of Sheffield, 2009.

26.

Nixon

O’Hagan

Oakley

. The Rheumatoid Arthritis Drug Development Model: a case study in Bayesian clinical trial simulation. Pharm Stat. 2009;8:371–89.

27.

O’Hagan

Forster

. Kendall’s Advanced Theory of Statistics: Bayesian Inference, Vol 2B, 2nd ed. London: Arnold; 2004.