An algorithm to generate correlated input-parameters to be used in probabilistic sensitivity analyses.

Background: Assessment of uncertainty in cost-effectiveness analyses (CEAs) is paramount for decision-making. Probabilistic sensitivity analysis (PSA) estimates uncertainty by varying all input parameters simultaneously within predefined ranges; however, PSA often ignores correlations between parameters. Objective: To implement an efficient algorithm that integrates parameter correlation in PSA. Study design: An algorithm based on Cholesky decomposition was developed to generate multivariate non-normal parameter distributions for the age-dependent incidence of herpes zoster (HZ). The algorithm was implemented in an HZ CEA model and evaluated for gamma and beta distributions. The incremental cost-effectiveness ratio (ICER) and the probability of being cost-effective at a given ICER threshold were calculated for different levels of correlation. Five thousand Monte Carlo simulations were carried out.
Results: Correlation coefficients between parameters sampled from the distribution generated by the algorithm matched the desired correlations for both distribution functions. With correlations set to 0.0, 0.5, and 0.9, 90% of the simulations showed ICERs below $25,000, $33,000, and $38,000 per quality-adjusted life-year (QALY), respectively, varying incidence only; and below $38,000, $48,000, and $58,000 per QALY, respectively, varying most parameters.
Conclusion: Parameter correlation may impact the uncertainty of CEA results. We implemented an efficient method for generating correlated non-normal distributions for use in PSA.

Background

Cost-effectiveness analyses (CEAs) are widely used to assess the cost-effectiveness ratio of a new health technology compared to an alternative intervention [1]. Assessment of uncertainty in CEAs is paramount since it informs decision makers whether results are robust under changes in assumptions and input parameters. Uncertainties in CEAs may stem from the modeling approach (structural uncertainty) or from uncertainty in input parameters, either because of imperfect knowledge of their true value or variability due to inherent heterogeneity, for example in patient characteristics. Parameter uncertainty may virtually affect all input parameters including demographic, epidemiological, clinical, economic, and quality of life parameters. Deterministic sensitivity analysis (DSA) and probabilistic sensitivity analysis (PSA) are often used to explore the impact of input parameters on CEA outcomes [2-5].

While DSA explores to which extent CEA results change when assuming extreme values for one or multiple parameters [6], PSA provides a global picture of uncertainty by sampling parameter values, typically using Monte Carlo (MC) simulations, from probabilistic distributions (e.g., normal, log normal, beta, gamma) [6]. In most instances, the correlation between input parameters is not incorporated in PSA [7]. Treating all input parameters as independent variables may lead to unrealistic combinations of parameter values and joint distributions that do not represent adequately the real-world situation. The importance of considering the correlation between different parameters when sampling their value for PSA has been raised in the literature [6-8]. Additionally, various guidelines on economic evaluations also recommend that the correlation among parameters should be incorporated, as it can affect both expected values and their degree of uncertainty (i.e., understating or overstating the true cost uncertainty) [6,9]. However, the techniques and algorithms for handling parameter correlation in a cost-effectiveness model, especially in the case of multivariate non-normal distributions, are less well documented. Nevertheless, methods of generating multivariate distributions from correlated non-normal marginal distributions have been described in other fields, such as behavioral science and financial risk analysis.

In this study, we used an algorithm based on Cholesky decomposition to generate correlated multivariate non-normal distributions (e.g., beta, gamma, log-normal) for use in PSA. We evaluated the impact of correlation between input parameters on CEA results using as an example the cost-effectiveness of herpes zoster (HZ) vaccination in persons aged ≥50 years. Vaccination was shown to be highly cost-effective under model assumptions and PSA, but so far, the existing CEA model did not address the correlation between input parameters.

Several studies reported that HZ incidence is higher in the older age groups [10-12]. A simulation of the incidence by age group without taking into account parameter correlation may lead to inconsistent values where the incidence in older age groups may be lower than the incidence in younger age groups, thereby leading to modeling results that are inconsistent with epidemiological data.

The objectives of this study are to raise awareness of the importance of parameter correlation in CEA when evaluating uncertainty and to provide access to an algorithm that can be efficiently incorporated into existing CEA models.

Methods

Correlation of parameters

In PSA, each input-parameter is considered as a random variable, associated with a probability distribution. For each MC run, parameters are randomly sampled from their assigned distributions and the model is executed to calculate a set of outcome values (i.e., Incremental Cost-Effectiveness Ratio [ICER]). The procedure is repeated for numerous iterations, usually 1,000 or more, and the outcomes presented on a probabilistic sensitivity analysis scatter plot showing the incremental costs per quality-adjusted life-year (QALY) gained distributed across the cost-utility plane [13].

Most often, parameters are sampled independently, but this is inadequate in the presence of parameter correlation. In the case of HZ incidence, epidemiological data show that the risk of HZ increases with increasing age. The HZ CEA model used here is stratified by age, and each age group is associated with its expected HZ incidence and distribution. When sampling HZ incidence values for each age group, the correlation between sample values should be maintained to reflect the fact that HZ incidence increases with age. Therefore, these parameters should be sampled together from a joint probability distribution of correlated parameters.

Generating a joint probability distribution from a set of uncorrelated (marginal) parameter distributions can be achieved by different approaches, such as the copula approach [8]. Here, we use a similar, well-known algorithm, i.e., the Cholesky decomposition matrix, which allows introducing pre-specified correlation among a set of uncorrelated variables.

Let us assume that the variables are correlated using the correlation matrix:

with an associated Cholesky decomposition matrix L, such that [14,15].

If we assume that vector Y represents the three independent variables , then the vector X = YL represents three correlated variables with correlation matrix C [14,15].

Sampling algorithm

The matrix is known as the Cholesky decomposition matrix of the correlation matrix . To simplify, we assumed that the correlation values between parameters are equal (=. We then simulated the correlated multivariate non-normal distribution following these steps:

We created a Cholesky decomposition matrix using the correlation value and the number of parameters correlated ();

We simulated uncorrelated normal random variables with a and a ;

We multiplied the uncorrelated normal variables by the Cholesky matrix to generate correlated normal random variables and calculated the probability associated with each value from the cumulative distribution function from a normal distribution;

We calculated the inverse probability of these values assuming a non-normal distribution (e.g., beta, gamma, log-normal).

The algorithm was implemented using visual basic for applications (VBA) in Microsoft Excel (2007) and further explanations are provided in the user guide (Additional file 1).

Implementation

The algorithm was implemented in a static multi-cohort Markov model developed to assess the cost-effectiveness of an adjuvanted recombinant zoster vaccine (RZV) vs no vaccination in adults aged ≥50 years [16]. Input data used in the model were retrieved from the RZV cost-effectiveness study recently published [17]. Cohorts were split into five age groups (i.e., 50–59, 60–64, 65–69, 70–79, and ≥80). The CEA was performed from a US societal perspective over a lifetime horizon. The primary outcome in the model was the ICER. The model results and ICERs including uncertainty were estimated by carrying out two PSAs consisting of 5,000 simulations each. In the first PSA, only HZ incidence was varied using multivariate distributions of correlated parameters generated by Cholesky decomposition. In the second PSA, besides HZ incidence most other input parameters were varied, except for all-cause mortality, vaccine price per dose, administration cost per dose, number of doses, the time between first and second dose, coverage of first dose vaccination, number of doses and discounting.

To assess the performance of the Cholesky decomposition algorithm in generating joint distributions of variables with a pre-specified correlation, we first simulated sets of parameters where the correlation was set to the desired values of 0.25, 0.5, 0.75, and 0.9. Next, we evaluated the impact of parameter correlation on PSA results. Two PSAs were carried out as described above: in the first PSA only correlation between HZ incidence parameters was considered, while in the second PSA, most model parameters were treated as correlated parameters. Furthermore, three different scenarios were evaluated, depending on the desired correlation values:

Zero correlation: no correlation, parameters are assumed to be independent (;

Moderate correlation (;

Strong correlation (.

Key parameters in the model and the distribution associated with each parameter in the PSA are presented in Additional file 2. Gamma distributions were assumed for all cost variables while beta distributions were applied for all other inputs [4,5]. The results of the PSA were presented separately on a cost-effectiveness plane and a cost-effectiveness acceptability curve.

Results

The Cholesky decomposition matrix has been embedded in a Microsoft Excel workbook via a VBA macro and is provided in Additional file1. We used a 5 × 5 matrix to match the number of age groups (i.e., 50–59, 60–64, 65–69, 70–79, and ≥80) used in our CEA example of HZ vaccination. The user can integrate this tool in any Excel-based PSA project to take into account the correlation between a set of correlated parameters. The correlation coefficient and type of distribution can be adapted by the user.

Adding the correlation between selected parameters using Cholesky decomposition added very little computational overhead. The actual runtime will depend on the CEA model and the computational platform.

Simulation results

Cholesky decomposition generated multivariate parameter distributions with correlations between age-dependent HZ incidence values close to the desired correlations, regardless of the underlying marginal distribution (beta or gamma distribution). The median correlation estimated through the simulation assuming a gamma distribution was 0.247, 0.496, 0.748, and 0.899 for the desired correlations of 0.25, 0.5, 0.75, and 0.9, respectively. Similarly, the median correlation estimated assuming a beta distribution was 0.248, 0.494, 0.746, and 0.895, respectively (Additional files 3 and 4).

Cost-effectiveness results

The base case results (e.g., the results of the CEA model) showed that RZV was associated with an ICER estimated at 11,863 USD per quality-adjusted life-year (QALY) gained compared to no vaccination.

For the PSAs varying incidence only, 90% of the simulations resulted in an ICER below 25,000 USD, 33,000 USD, and 38,000 USD using a correlation of 0.0, 0.5, and 0.9, respectively (Figure 1(a)). For the PSAs varying most parameters in the model, 90% of the simulations resulted in an ICER below 38,000 USD, 48,000 USD, and 58,000 USD using a correlation of 0.0, 0.5, and 0.9, respectively (Figure 1(b)).

Figure 1.

Cost-effectiveness acceptability curves varying (a) incidence parameters only; or (b) most parameters

Uncertainty around the incremental QALYs and incremental costs for the two scenarios is presented by scatterplots (Figure 2(a,b)). The QALYs gained ranged from approximately 1200 to 3600 while the costs varied from approximately 100 million in expenditure to savings of 145 million per 1 million individuals vaccinated. As expected, uncertainty was larger in scenario 2 in which most parameters were varied compared with scenario 1, in which only a subset of parameters related to HZ incidence was varied. Including correlation between the age-dependent HZ incidence also tended to increase uncertainty as illustrated by dots lying on the outside of the scatterplot cloud.

Figure 2.

ICER scatterplot of scenarios for PSA varying (a) incidence parameters only; or (b) most parameters

Discussion

PSAs provide a synthesis of parameter uncertainty and its impact on the cost-effectiveness of a given intervention; guidelines recommend the use of PSA to assess uncertainty in CEA and PSAs have been widely used in health economics evaluations [4]. However, a review of CEAs conducted for the National Institute for Health and Clinical Excellence (NICE) revealed that the majority of models ignore the potential correlation between input parameters [7]. The apparent lack of considering parameter correlation in CEAs might be due to the unfamiliarity of the analyst with the subject, programming issues, and the absence of explicit requirements and recommendations in applicable guidelines. The goal of this study was to overcome these hurdles by describing a known algorithm for handling parameter correlation and providing access to its implementation in Microsoft Excel, a program widely used for CEA modeling. Similar efforts have been undertaken by several other groups, using different methods to introduce correlation in PSA, with accompanying code provided for the statistical package R [8,18].

Our algorithm used a Cholesky decomposition to transform independently correlated normal random variables into correlated normal random variables. These variables were then transformed into non-normal variables resulting in correlated non-normal random variables. We simulated 5,000 values for each parameter included in the PSA; after 5,000 simulations, the results appeared to be stable and no additional simulations were needed [19]. The results showed that the desired correlations were met by the algorithm, as the median values were close to the desired correlations for both distributions. The scatter plot clouds produced by varying only the incidence parameters are more densely distributed than when most parameters are varied, which reflects the impact of the different parameters on the PSA analysis.

In our CEA model of HZ vaccination, including the (positive) correlation between age-dependent HZ incidence lead to higher uncertainty in CEA outcomes, i.e., incremental costs and incremental QALYs. As expected, the credibility interval changed dramatically when more parameters were introduced in the PSA. Other situations may occur where consideration of the correlation between input parameters may decrease the uncertainty in cost-effectiveness results [8]. In either case, the impact of parameter correlation provides a more realistic estimate of the uncertainty in CEA results. PSA may also highlight areas where additional information might be valuable to reduce uncertainty and thereby lending more credibility to model outcomes.

There are several limitations to the presented method. The results provided by this algorithm are imperfect in terms of the desired correlations and the correlations obtained based on simulations, e.g., for the desired correlation of 0.25 and assuming a beta distribution, the median was 0.247 with a range of 0.184 to 0.277 (Additional file 3). The skewed nature and ranges observed are similar to those presented by Iman and Conover 1982 [14].

Frequently, the true value of a correlation is unknown. Ideally, all available evidence, including expert opinion, should be used to estimate uncertainty both around input parameter values and correlation. While the presented method allows for propagating uncertainty in input parameters to assess uncertainty in CEA output variables, the correlation between input variables is considered as a fixed variable in the PSA. Other methods of simulating multivariate non-normal variables have been proposed which are more robust and allow for considering uncertainty in correlation itself, but these methods require iterative processes and are consequently less efficient from a computational perspective [6,18]. The described method is intended for situations where the exact correlation remains unknown, but where several plausible assumptions can be explored. Another consideration is the potential non-linearity of relationships, which may influence model calculations [8].

In conclusion, considering correlations in PSA may impact conclusions reached based on CEA results. Analysis of decisions taken by NICE on a series of CEA models suggests that the outcome depended on uncertainty and that high uncertainty was more likely leading to a negative decision [7]. Although the precise value of parameter correlation is often unknown, allowing for parameter correlation will help decision makers understand the potential impact on results and overall model uncertainty. Here, we developed an algorithm based on Cholesky decomposition, which simulates a correlated non-normal distribution taking any desired correlation value. The VBA code can be freely downloaded from the link attached to this paper.

Click here for additional data file.