Calculating the Expected Value of Sample Information using Efficient Nested Monte Carlo: A Tutorial

Objective: The Expected Value of Sample Information (EVSI) quantifies the economic benefit of reducing uncertainty in a health economic model by collecting additional information. This has the potential to improve the allocation of research budgets. Despite this, practical EVSI evaluations are limited, partly due to the computational cost of estimating this value using the"gold-standard"nested simulation methods. Recently, however, Heath et al developed an estimation procedure that reduces the number of simulations required for this"gold-standard"calculation. Up to this point, this new method has been presented in purely technical terms. Study Design: This study presents the practical application of this new method to aid its implementation. We use a worked example to illustrate the key steps of the EVSI estimation procedure before discussing its optimal implementation using a practical health economic model. Methods: The worked example is based on a three parameter linear health economic model. The more realistic model evaluates the cost-effectiveness of a new chemotherapy treatment which aims to reduce the number of side effects experienced by patients. We use a Markov Model structure to evaluate the health economic profile of experiencing side effects. Results: This EVSI estimation method offers accurate estimation within a feasible computation time, seconds compared to days, even for more complex model structures. The EVSI estimation is more accurate if a greater number of nested samples are used, even for a fixed computational cost. Conclusions: This new method reduces the computational cost of estimating the EVSI by nested simulation.


Introduction
The Expected Value of Sample Information (EVSI) [ 1] uses evidence about the cost and effectiveness of new treatments to determine the expected economic benefit of undertaking a proposed study [ 2]. The EVSI calculates this value by determining the extent to which the additional information from the study reduces the probability and expected loss of making inefficient treatment recommendations [ 3]. In general, an inefficient recommendation would spend health resources that would be used to improve patient outcomes elsewhere. Additionally, information from future studies has the potential to improve patient outcomes within the disease area under investigation.
As the EVSI calculates the value of a specific trial, it has the potential to be used to determine the optimal allocation of research funding. Despite this, the EVSI has rarely been used in practical scenarios [ 4]. This is partly due to the large computational effort required to calculate the EVSI using the "gold-standard" nested Monte Carlo method [ 2].
Recently, a number of methods have been developed to reduce the computation time for the EVSI [ 5,6,7,8,9,10]. While computationally efficient, these methods require the use of additional statistical techniques, e.g. specifying sufficient statistics [ 11], and may impose restrictions on the underlying economic model. Therefore, EVSI calculations have been repeatedly undertaken using Monte Carlo procedures, either using approximations to avoid nested simulations [ 12,13,14], or with high computational cost [ 15,16,17,18] Using nested Monte Carlo methods is advantageous as it can be used irrespective of the complexity of the economic model. Additionally, the EVSI is often understood in terms of this nested Monte Carlo estimation procedure [ 19,20,5,8] making the computation method easier to comprehend. To utilise these advantages, Heath et al developed a computationally efficient EVSI calculation method based solely on nested sampling [ 10] to calculate the EVSI in all practical scenarios, within a reasonable timeframe.
The validity of this method has been demonstrated elsewhere [ 10] and thus this paper is focused on presenting its practical implementation. Section 2 introduces the EVSI, key notation and the nested Monte Carlo method. Section 3 outlines the efficient Monte Carlo method using a simple "toy". Finally, Section 4 investigates the optimal choice of the key input for this efficient Monte Carlo method using an economic model evaluating a new Chemotherapy treatment [ 21].

Formal definition of the EVSI
In general, information has value as it reduces uncertainty in the inputs of a health economic model, denoted . Therefore, to calculate the EVSI, we must model the current level of uncertainty in using probability distributions that indicate the most plausible range for these inputs. This is performed as part of a Probabilistic Sensitivity Analysis (PSA) [ 22,23] where the parameters of these distributions are usually informed by literature reviews alongside clinical trials and expert opinion.
To value each of the treatment options, normally using a net benefit function [ 24], the model inputs are combined using an economic model. To ease our explanation, we assume throughout that this economic model only compares two treatment options -typically, an innovative treatment and the standard of care. The output of the economic model is then the incremental net benefit (INB), i.e. the difference between the net benefit of the innovative treatment and net benefit of the standard of care. In this setting, uncertainty in the model inputs implies that the INB is also uncertain. From this, we conclude that the innovative treatment is optimal given current uncertainty if the average INB is greater than 0 and the standard of care is optimal otherwise [ 22]. where the value of the optimal treatment is some monetary value, say ‫.ܥ$‬ To calculate the EVSI, we consider that additional information, denoted by , is going to be collected in a future study, e.g. a clinical trial or an observational study, to reduce the uncertainty in some of the model inputs. If had been collected, i.e. the future study had already been conducted, then the decision making process would be exactly as we have already described. The information in would be formally included in the PSA distributions and the value of the optimal decision given this additional information would be either 0, if the average INB is negative, or the average INB itself, equal to ‫ܨ$‬ in the bottom of Figure 1. Comparing ‫ܨ$‬ with ‫ܥ$‬ gives the Value of the Sample Information (VSI).
In general, we can consider the PSA distributions to be the prior for , which is combined with the data using Bayes Theorem to compute the posterior distribution for . This in turn implies a posterior distribution for the INB.
However, as the future study has not been run (and potentially will never be), we define a distribution over the possible, yet unknown study outcomes. The EVSI is then equal to the average VSI over all these possible future datasets. Mathematically, it can be expressed in terms of the INB as the posterior expectation of the INB for a specific sample . Typically, will only directly update a small number of model parameters, but in a Bayesian setting we still take expectation with respect to the joint posterior distribution ‫.)|(‬ To model the possible study outcomes, we specify the sampling distribution for the data, had they been collected. This will depend on some of the model inputs: for example, a binomial distribution models the number of people responding to a treatment conditional on the success rate of the drug. Combining this sampling distribution, conditional on the model inputs, with the PSA distributions for the model inputs gives a distribution over the potential future datasets. This process is represented by the red arrow in Figure 1. PSA is normally undertaken using a simulation approach [ 22,21,25] where ܵ Consequently, to calculate the EVSI by Monte Carlo requires ܵ × ܴ simulations, which can be relatively computationally expensive for standard choices of ܵ (normally around 1000) and ܴ (at least 600 [ 26]).
The method presented in this paper uses a nested Monte Carlo scheme but finds the posterior distribution of the model inputs for a small number of potential future datasets.
Specifically, the required number of posterior updates reduces from ܵ to ܳ < 100. This maintains the flexibility of the nested Monte Carlo method whilst drastically reducing its computational cost, irrespective of the model structure.

The moment matching method
The efficient nested Monte Carlo method [ 10] is based on "moment matching" and requires several elements. To begin, we must estimate the mean ߤ and variance ߪ ଶ of the INB, using the PSA simulations.

Uncertainty due to key parameters
As previously discussed, the sampling distribution of , typically, only depends directly on a small number of the model inputs. For example, within a full health economic model, a clinical trial would focus on the epidemiological parameters and not the economic disease burden. To formalise this, we assume that the model inputs split into two categories, (ࣘ , ), where the sampling distribution for is based solely on the inputs ࣘ (e.g. the epidemiological parameters) and all the remaining inputs are in . The moment matching method requires the distribution of the INB where uncertainty due to the model inputs in has been marginalised out. This is expressed mathematically as Notice that this quantity is required to calculate an alternative Value of Information measure known as the Expected Value of Partial Perfect Information (EVPPI) [ 6], which quantifies the economic value of learning the exact value of the model inputs ࣘ . As study information cannot give exact information, the EVPPI for ࣘ is an upper bound for the EVSI. If this upper bound is low then there is no value in a study targeting ࣘ and so the modeller can discount a trial targeting ࣘ before determining a sampling distribution for . Therefore, the EVPPI should always be calculated before proceeding to the EVSI [ 27], which means that the  To accurately estimate the EVSI, the variance of the posterior INB must be estimated for a number of potential datasets, say ܳ. Theoretically, ܳ should be greater than 30, as accurate estimation of the variance depends on the central limit theorem, but, the variance of the posterior INB is sufficiently accurate for values of ܳ close to 30 [ 10], provided the future datasets are simulated using the following procedure.
Loosely, we want to "space-out" the simulated future datasets to accurately estimate the EVSI. To achieve this, we find ܳ equally spaced values of ࣘ by determining the quantiles of the PSA simulations and use these to generate the future datasets. Practically, we proceed by ordering the PSA simulations for each of the model inputs in ࣘ and selecting the

Calculating the EVSI
To calculate the EVSI, we must combine all these different elements: Finally, all the above elements are used to rescale the INB ࣘ values: This gives ܵ rescaled INB ࣘ values, where ܵ is the number of PSA simulations, which are then used to estimate the EVSI,

Toy Example
To clarify the moment matching method, we estimate the EVSI in a very simple setting.
The effectiveness measure in this example is the probability of curing a disease. Information from a previous trial implies that the uncertainty in this effectiveness measure for two treatments under consideration can be modelled using the following Beta distributions ߨ ଵ~‫ܤ‬ ‫)4,3(ܽݐ݁‬ and ߨ ଶ~‫ܤ‬ ‫.)3,4(ܽݐ݁‬ A literature review determines that uncertainty in the incremental cost can be modelled using a Normal distribution ‫݉ݎܰ~ߜ‬ ݈ܽ (3,20). Therefore, in this example, = (ߨ ଵ , ߨ ଶ , ߜ). We then define the INB for this model as where 100 is the selected willingness to pay. Note that this threshold is used for illustrative purposes as the effectiveness measure is probability of a cure rather than Quality Adjusted Life Years (QALYs).
For this model, we report PSA simulations for ܵ= 10 in Table 1. Each parameter is simulated from its distribution and then each simulated row is used to calculate the INB using equation (1). The mean and variance of the INB are calculated as ߤ = −4.5 and ߪ ଶ = 722.
Using an EVPPI analysis, we found that reducing uncertainty in ߨ ଵ , the probability of a cure for treatment 1, has the highest value. Therefore, we design a trial where treatment 1 is given to 20 people and observe how many patients respond. This implies a Binomial sampling distribution for the future study: ‫݅ܤ~‬ ‫݉݊‬ ݅ ݈ܽ (20, ߨ ଵ ). The final element required for the moment matching method is the posterior variance estimated using nested sampling. For illustrative purposes, we set ܳ= 3. To begin, we order the observed ߨ ଵ values, which produces the following vector: These values are given in Table 1. Finally, the EVSI is estimated using moment matching, as EVSI = 1 10 (4 + 10 + 29 + 2) − max{0, −4.5} = 4.6 This example is illustrative of the moment matching method; in practice a large PSA simulation size should be used, ideally in excess of 1000.

Calculating the EVSI: a new chemotherapy drug
To demonstrate the moment matching method and investigate the optimal value for ܳ, we extend a model developed to compare a new chemotherapy treatment against the standard of care [ 21].

The optimal number of nested samples
To find the optimal number of nested samples, we fix the total number of simulations used to estimate the EVSI. Clearly, as with all simulation based methods, we increase the accuracy of the EVSI estimate by increasing the total number of simulations. However, in fixing the total computational cost we determine the relative importance of increasing ܳ versus increasing the number of simulations from the posterior distribution of the model inputs.
We use two sizes of nested simulation, 500 000, with a computation time of between 108 and 129 seconds for one EVSI estimate and 5 000, with a computation time of between 9 and 17 seconds. We then considered 8 different values of ܳ = 20, 30, … , 100, where 20 is below the recommended lower limit for ܳ [10]. As the EVSI is estimated by simulation, it is subject to random variance. Therefore, we calculated the EVSI 200 times for each combination of ܳ and simulation size. This allows us to calculate the variance of the EVSI estimate for each combination ( Figure 2) and the bias (Figure 3). To calculate the bias, we estimated the EVSI using nested simulation with ܵ = 100 000 and ܴ = 100 000 with a computational cost of around 60 days. Figure 2 demonstrates that larger ܳ values produce more precise EVSI estimates as the variance is smaller. This is true for both 500 000 and 5 000 simulations although the variance seems to plateau earlier for 5 000 implying that increasing ܳ is less important when the number of posterior simulations is not sufficient to estimate the posterior of . Therefore, computational effort should be used to increase ܳ and not ܴ, the number of simulations for each posterior.
Clearly, the variance of the estimated EVSI is lower with a greater number of simulation. Therefore, accurately estimating ߪ ଶ using nested simulations is useless if ߪ ଶ is poorly estimated.
Thus, the PSA simulations size ܵ must be relatively large. Additionally, the required accuracy of the estimates of ߪ ଶ and ߪ ଶ actually depends on the size of the EVSI. For a small EVSI, the posterior variance is close to the prior variance meaning that of ߪ ଶ − ߪ ଶ is small. As both variances are estimated by simulation, the sampling variation could be larger than the "true" difference, leading to inaccurate EVSI estimates. Therefore, the moment matching method should only be used for studies where the EVPPI for ࣘ is high.

Conclusion
In this paper, we have presented an efficient Monte Carlo estimation method for the EVSI. This method significantly reduces the number of nested simulations required to accurately estimate the EVSI. This maintains the flexibility and comprehensibility of the nested Monte Carlo method whist reducing computational time. The moment matching method was presented in theoretical terms alongside an illustrative worked example. Finally, we explored the optimal number of nested samples required to accurately estimate the EVSI. In general, the number of nested samples should be taken as large as possible, whilst maintaining a feasible computational cost and ensuring that the posterior distribution for the model inputs is adequately captured.
Another advantage of this efficient Monte Carlo method is that general purpose software has been developed to perform these EVSI calculations [ 34]. Early results also indicate that this method can be extended to estimate EVSI across different sample sizes with no additional computational cost meaning that the EVSI could be used as a tool for trial design within a realistic time frame.