Probabilistic threshold analysis by pairwise stochastic approximation for decision-making under uncertainty

The concept of probabilistic parameter threshold analysis has recently been introduced as a way of probabilistic sensitivity analysis for decision-making under uncertainty, in particular, for health economic evaluations which compare two or more alternative treatments with consideration of uncertainty on outcomes and costs. In this paper we formulate the probabilistic threshold analysis as a root-finding problem involving the conditional expectations, and propose a pairwise stochastic approximation algorithm to search for the threshold value below and above which the choice of conditionally optimal decision options changes. Numerical experiments for both a simple synthetic testcase and a chemotherapy Markov model illustrate the effectiveness of our proposed algorithm, without any need for accurate estimation or approximation of conditional expectations which the existing approaches rely upon. Moreover we introduce a new measure called decision switching probability for probabilistic sensitivity analysis in this paper.


EVPPI versus decision switching probability
Here we show an example of the two-alternative decision models for which the larger EVPPI does not necessarily mean the larger decision switching probability. Let D = {d 1 , d 2 } and θ = (θ 1 , θ 2 ) with θ 1 and θ 2 independently following the standard normal distribution. The functions f d1 and f d2 are given by respectively, with some constant 0 < α < 1.
Since Let us consider the partial perfect information on θ 1 . As we have the conditionally optimal option becomes d opt (θ 1 ) = d 1 for − √ 1 + α < θ 1 < √ 1 + α, d 2 for θ 1 where d opt (θ 1 ) is not unique for θ 1 = ± √ 1 + α. Now the decision switching probability for θ 1 is given by where erf denotes the error function, while the EVPPI for θ 1 is given by A similar computation leads to the decision switching probability for θ 2 : There results are compared in Supplementary Figure 1. The left panel shows the behaviors of the decision switching probability for θ 1 (in blue) and θ 2 (in orange) as functions of α ∈ (0, 1), respectively, whereas the right panel shows those of the EVPPI. The decision switching probability monotonically decreases both for θ 1 and θ 2 as α increases. It is important that the decision switching probability for θ 2 is much larger than that for θ 1 when α is small. This means that knowing the value of θ 2 exactly has a higher chance of changing the optimal option from d 1 to d 2 . On the contrary, knowing the exact value of θ 1 has a high chance of keeping the optimal option unchanged, making the partial prefect information useless in terms of decision making. The two curves intersect at α ≈ 0.958362, beyond which the decision switching probability for θ 1 becomes larger.
The EVPPI decreases both for θ 1 and θ 2 as α increases. Such a behavior is consistent with the decision switching probability. However, for all α ∈ (0, 1), the EVPPI for θ 1 is larger than that for θ 2 , which is opposite to the decision switching probability when α < 0.958362. From these results, we see that knowing the exact value of θ 1 has a low chance of changing the optimal option from d 1 to d 2 but an increment of the expected cost-effectiveness gained from such decision switching is more significant than knowing the exact value of θ 2 . Through this example, we can conclude that the larger EVPPI does not necessarily mean the larger decision switching probability and vice versa.

Problem setting
Here we test our proposed approach for an additional example introduced by Heath and Baio [2]. The model compares two different chemotherapy treatments: a standard of care (d 1 ) and a novel treatment (d 2 ). The cost and the quality-adjusted life year (QALY) for each treatment are modelled by a homogeneous, discrete-time Markov process with four states (recorvery, home care, hospital care, and death) depending on whether patients experience side effects or not. The difference between two treatments is in the probability of experiencing side effects. Here π d1 denotes the probability of experiencing side effects for the standard of care, and ρ denotes the relative percentage in the number of side effects for the new treatment from that for the standard of care. That is, the probability of experiencing side effects for the new treatment is given by the product ρπ d1 . The input vector θ consists of 12 random variables π d1 , ρ, Γ 1 , Γ 2 , λ 1 , λ 2 , C Death , C Hosp , C Home , Q Recov , Q Hosp , Q Home , which are denoted by θ 1 , . . . , θ 12 , respectively, in this order. The model contains three constants C d1 , C d2 , λ, while two parameters γ 1 and with SE d1 ∼ Bin(10 3 , π d1 ) and SE d2 ∼ Bin(10 3 , ρπ d1 ) being the conditionally independent random variables given π d1 , where Bin(n, p) denotes the Binomial distribution with n independent trials and success probability p for each trial. Moreover, the matrix M ∈ R 4×4 ≥0 denotes the transition matrix of the time-homogeneous Markov process given by It is clear that the functions f d1 and f d2 are non-linear with respect to some of input variables, π d1 , ρ, Γ 1 , Γ 2 , λ 1 and λ 2 , whereas they remain linear with respect to the other variables. Even through the input variables are assumed mutually independent, such non-linearity of the functions makes it quite hard to compute the probabilistic parameter threshold exactly. Nevertheless, since only two treatments are compared in this model, the probabilistic parameter threshold can be estimated by solving a single root-finding problem without any need for hypothesis testing, that is, the fourth item in Algorithm 2 (from the main text [1]) can be skipped.

Results and discussion
We consider estimating the probabilistic parameter thresholds for 4 input variables θ 1 , θ 2 , θ 5 and θ 6 , respectively. We use Algorithm 2 (from the main text [1]) without the fourth step and with M = 1, T = 10 4 to estimate the threshold K j . In the first item of Algorithm 2, we generate θ 1 j from the marginal distribution of θ j . We set the sequence of step sizes to where σ(θ j ) denotes the standard deviation for the marginal distribution of the variable θ j . As the variables θ 1 , θ 5 and θ 6 are defined over the interval [0, 1], if the recursion (2) in the main text [1] yields θ t+1 j < 0 or θ t+1 j > 1, we put θ t+1 j = 0 or θ t+1 j = 1, respectively. Moreover, as the product θ 1 θ 2 (= ρπ d1 ) is used as a parameter to generate a binomial random variable SE d2 , we must have 0 ≤ θ 1 θ 2 ≤ 1. In our experiments, instead of restricting the support of θ 2 , we simply consider SE d2 ≡ 0 if θ 1 θ 2 < 0 and SE d2 ≡ 10 3 if θ 1 θ 2 > 0. The averaged outputs Θ t j are considered as a sequence of our threshold estimates. We carry out 20 independent runs for each considered variable.
Supplementary Figure 2 shows the convergence behaviors of the estimates Θ t j as functions of t, obtained from the second item of Algorithm 2. For all 20 runs, the estimates converge to a similar value, indicating that the threshold K j consists of only one element. The standard error decays with the convergence rate of approximately t −1/2 , or even better for the variables θ 2 and θ 6 . This way our proposed approach works quite well not only for the simple synthetic testcase but also for such a complicated Markov model.
The second column of Supplementary Table 2 summarizes the final estimates for the threshold K j obtained by our proposed approach after T = 10 4 iteration steps. Using these results, the intervals of θ j where d opt (θ j ) is equal to d 1 and d 2 are identified as shown in the third and fourth columns, respectively. For this model we have d opt (∅) = d 1 , so that the decision switching probability for a variable θ j is given by  which can be computed by the marginal cumulative distribution function. The resulting values are shown in the fifth column. The variable θ 6 yields the minimum decision switching probability, whereas the variables θ 2 and θ 5 yield the approximately equal probabilities about 0.2. As a reference, we estimate the EVPPI for each variable by using the nested Monte Carlo estimator with 2 18 outer samples for θ j and 2 10 inner samples for θ −j . Note that the EVPI is estimated to be 36.81 by using the standard Monte Carlo estimator with 2 24 random samples for θ. The estimated EVPPI values are shown in the sixth column. The variables θ 1 and θ 6 for which the decision switching probabilities are small yield the relatively small EVPPI values. Although the decision switching probabilities for θ 2 and θ 5 are almost the same, the EVPPI for θ 5 is twice larger than that for θ 2 . It is interesting to see that a variable with larger EVPPI does not necessarily have a larger decision switching probability not only for an artificial model considered in Supplementary Information 1 but also for a more realistic model for health economic evaluations.