Accounting for parameter uncertainty in the definition of parametric distributions used to describe individual patient variation in health economic models

Background Parametric distributions based on individual patient data can be used to represent both stochastic and parameter uncertainty. Although general guidance is available on how parameter uncertainty should be accounted for in probabilistic sensitivity analysis, there is no comprehensive guidance on reflecting parameter uncertainty in the (correlated) parameters of distributions used to represent stochastic uncertainty in patient-level models. This study aims to provide this guidance by proposing appropriate methods and illustrating the impact of this uncertainty on modeling outcomes. Methods Two approaches, 1) using non-parametric bootstrapping and 2) using multivariate Normal distributions, were applied in a simulation and case study. The approaches were compared based on point-estimates and distributions of time-to-event and health economic outcomes. To assess sample size impact on the uncertainty in these outcomes, sample size was varied in the simulation study and subgroup analyses were performed for the case-study. Results Accounting for parameter uncertainty in distributions that reflect stochastic uncertainty substantially increased the uncertainty surrounding health economic outcomes, illustrated by larger confidence ellipses surrounding the cost-effectiveness point-estimates and different cost-effectiveness acceptability curves. Although both approaches performed similar for larger sample sizes (i.e. n = 500), the second approach was more sensitive to extreme values for small sample sizes (i.e. n = 25), yielding infeasible modeling outcomes. Conclusions Modelers should be aware that parameter uncertainty in distributions used to describe stochastic uncertainty needs to be reflected in probabilistic sensitivity analysis, as it could substantially impact the total amount of uncertainty surrounding health economic outcomes. If feasible, the bootstrap approach is recommended to account for this uncertainty. Electronic supplementary material The online version of this article (doi: 10.1186/s12874-017-0437-y) contains supplementary material, which is available to authorized users.


Additional file 1.1: Detailed description for the Bootstrap approach
Non-parametric bootstrapping is a statistical technique that can be used to construct an approximate sampling distribution for a statistic of interest, without the need for assumptions regarding the distribution of this statistic [1].
Several studies investigated the use of bootstrapping in health economics, e.g. for constructing confidence intervals for the incremental cost-effectiveness ratio and the incremental net benefit [2][3][4]. For a parameter of interest, we are interested in the value for the whole population β, which cannot be observed. Therefore, we try to find information about the value of this population parameter by drawing a random sample Y from this population and estimate ̂( Y) the parameter of interest based on this sample. The use of bootstrapping enables us to find information about the relation between the population parameter β and its estimate ̂( Y) by the relationship between an observed value for the parameter of interest ̂( ) and a value for the parameter of interest based on a bootstrap sample ̂( Y*) [1,3].
If non-parametric bootstrapping is applied, bootstrap sample Y* is constructed by resampling from the observed sample yobs with replacement [1,3].
The reasoning for applying non-parametric bootstrapping to reflect the uncertainty in distributions' parameter estimates is explained by the fact that estimates from a clinical trial are almost always obtained based on a part of the total population. Hence, it would be incorrect to assume that the distributions' parameters values estimated from the data are known with certainty, i.e. are correctly describing the entire population. There is a certain sampling error, as another clinical trial may yield different estimates for the distributions' parameters. By bootstrapping trial data, other trial results are simulated and other estimates of parameters will be found, which may generate different probabilistic sensitivity analysis outcomes when Monte Carlo simulation is applied. In fact, as the size n of the sample increases, the estimates obtained by bootstrapping converge to the population value [1,5].
In the Bootstrap approach the distributions' parameters are repeatedly, say r times, estimated based on different bootstrap samples * of the original data set , resulting in a set of {̂ * , ̂ * , …, ̂ * } parameter bootstrapestimates, where i equals the number of parameters required to define the distribution chosen to describe the stochastic uncertainty in the time-to-event data. These bootstrap samples * are obtained by resampling the original data set with replacement, such that the sizes of the original data set and the bootstrap sample are the same. These r sets of parameter bootstrap estimates ̂ * can then be used in the PSA, incorporating one set in each Monte Carlo sample.
In short, this approach can be translated into the following four steps: (1) Generate a feasible bootstrap sample of the original dataset, by resampling this dataset with replacement, such that the sample size of the bootstrap sample equals that of the original dataset.
(2) Fit the pre-specified distribution(s) to the bootstrap sample and record the estimated parameter values, e.g. of the shape and rate parameters of a Gamma distribution.
(3) Repeat (1) and (2) r times, where r equals the required number of PSA runs.
(4) Perform the PSA, using a different set of the r parameter values to define the distribution(s) for each PSA run.
Note that if multiple distributions are fitted in step (2) all of them need to be fitted on the same bootstrap sample to preserve correlation among distributions, which also applies to non-time-to-event distributions, such as Beta distributions to describe utilities, and other model parameters, such as probabilities.
The definition of a feasible bootstrap sample, as required in step (1), may vary between studies and can be difficult to decide on. When repeatedly fitting distributions to these random bootstrap samples, there might be samples for which it is not possible to fit a distribution. For example, a bootstrap sample might be drawn that contains only one time-to-event observation for a specific event, which makes it impossible to fit a Weibull distribution for the time to this event. This particular bootstrap sample can be considered infeasible for fitting a Weibull distribution, though it does contain information on the event of interest and, therefore, might be considered feasible. In such a scenario, it needs to be decided whether the bootstrap sample is to be excluded and a new sample will be drawn. Especially in case of small sample sizes or rare events these situations are likely to occur. Modelers need to be aware of this decision and should clearly communicate their choices in the corresponding publications.

Additional file 1.2: Detailed description for the MVNorm approach
The uncertainty in distributions' parameter estimates can be described by assuming these parameter estimates to be Normal distributed. This is appropriate for sufficiently large sample sizes according to the Central Limit Theorem, which states that for any population distribution of a parameter of interest, the distributions of the sample means ̅ will converge to Normal distributions as the sample size increases [4]. For a sufficient large sample size this indicates that the uncertainty surrounding an independent parameter can be defined according to its estimate ̂( ) and the Standard Error of the Estimate However, when estimating a multi-parameter distribution, these distribution's parameters are likely to be correlated and it is, therefore, incorrect to define separate Normal distributions using Equation 1 for each of the parameters.
However, multivariate Normal distributions can be used to draw correlated values for the parameters of interest.

Equation 2
Using multivariate Normal distributions, i.e. the MVNorm approach, the uncertainty surrounding the parameter estimates can be approximated using the following four steps: (1) Fit the pre-specified distribution to the original dataset and record the estimated parameter values, e.g.
of the shape and rate parameters of a Gamma distribution, and (calculate) the variance-covariance matrix.
(2) Define a multivariate Normal distribution from the parameters' estimates and their variance-covariance matrix according to (1).
(3) Draw r feasible sets of parameter values from the defined distribution (2), where r equals the required number of PSA runs.

Additional file 1.4: Description of the model used in the case study
The discrete event simulation (DES) model that was used in the case study was defined on patient-level using AnyLogic software and according to the ISPOR-SMDM Modeling Good Research Practice Task Force guidelines [6].
The model was defined to have the same health states as the model that was used for the original evaluation of the CAIRO3 study: post-induction, re-induction, salvage, and death (see figure below) [7].
Weibull distributions [8] were used to define all health state-specific time-to-event parameters and were estimated from the CAIRO3 trial data using the fitdist function of the fitdistrplus [9] package in R Statistical Software [10]. Events, i.e. transitions between health states, were based on patient-specific processing times, which were randomly drawn from the estimated Weibull distributions. Competing risks were handled by selecting the first event to occur based on the respective observed event probabilities and their corresponding state-specific time-to-event distributions [11]. For example, for a patient that is entering the re-induction state a random number was compared to the chance of progression to determine whether the patient would survive and progress to the salvage therapy state.
After the event was selected, the time to that event was randomly drawn from the corresponding Weibull distribution, i.e. time-to-progression or time-to-death.
In the DES model, 10,000 patients were simulated per treatment strategy. Patient-level outcomes were calculated based on the time patients had spent in each health state and were summarized to enable comparison of the two treatment strategies on population level. The DES model was validated according to good practices guidelines by assuring no unnecessary detail was present, structured "walk-throughs", comparing results with calculations by hand, extreme value analysis, trace analysis, sensitivity analysis, and cross validation with the model that was used for the original evaluation of the CAIRO3 study [12,13].

Post-Induction Re-Induction Salvage Death
Graphical representation of the model used in the case study