Scientific Papers

# Stochastic cost-effectiveness analysis on population benefits | Cost Effectiveness and Resource Allocation In order to further clarify the properties of the method proposed in this paper and demonstrate its association and difference with classical CEA tools, we perform numerical simulations to calculate the results, differences, and effectiveness of various types of uncertain CEA tools. Further, we take an empirical example to show that the mixed structure of interventions could get a higher population benefit, and the tool in this paper helps to find it.

### Stochastic system settings

We choose the mean–variance system to describe the stochasticity of the simulation system to simply and clarify the simulation processes and results. It is assumed that the joint distribution of the costs and QALY outcomes of individuals receiving various interventions follows an overall joint multinormal distribution. Therefore, characterizing all the stochasticity can be accomplished by setting the means, standard deviations, and the matrix of the correlation coefficient. For each intervention, it is assumed that the costs and QALY outcomes of each patient taking this intervention are identically distributed but not independent. And, there is a certain correlation between the costs and QALY outcomes of individual patients. The correlation between patient costs of different interventions is only related to the intervention they take, but without distinguishing individuals. Same assumption applies to QALY outcomes. In addition, it is assumed that the cost of any patient is independent of the QALY outcome of other patients, which means that the cost of patient A does not correlate with the QALY outcome of any patient B.

Correlation matrix can express the structure of correlation much clearer. Specifically, for each group of patients taking the same intervention, the costs and QALY outcomes of each patient have unique mean and variance values. That is, for any i and mi, we have

\begin{aligned}& {{\mathbb{E}}}[{\text{QALY}}_{m^i }^i ] = \mu_Q^i ,\;{\text{ Var}}[{\text{QALY}}_{m^i }^i ] = (\sigma_Q^i )^2 {{\ominus }}, \\ &{{\mathbb{E}}}[{\text{Cost}}_{m^i }^i ] = \mu_C^i ,\quad\; {\text{ Var}}[{\text{Cost}}_{m^i }^i ] = (\sigma_C^i )^2 . \hfill \\ \end{aligned}

(19)

And for the correlation coefficient of each cost, there is

$$\forall i,j,\;{\text{Corr}}[{\text{Cost}}_{m_1 }^i ,{\text{Cost}}_{m_2 }^j ] = \left\{ {\begin{array}{ll} {1,} &\quad {i = j,\;m_1 = m_2 } \\ {\rho_C^{ii} ,} &\quad {i = j,\;m_1 \ne m_2 } \\ {\rho_C^{ij} ,} &\quad {i \ne j,} \\ \end{array} } \right.$$

(20)

and for each QALY there is a similar equation

$$\forall i,j,\;{\text{Corr}}[{\text{QALY}}_{m_1 }^i ,{\text{QALY}}_{m_2 }^j ] = \left\{ {\begin{array}{ll} {1,} &\quad {i = j,\;m_1 = m_2 } \\ {\rho_Q^{ii} ,} &\quad {i = j,\;m_1 \ne m_2 } \\ {\rho_Q^{ij} ,} &\quad {i \ne j.} \\ \end{array} } \right. \,$$

(21)

In addition, the cross-sectional relationship between cost and QALY output is assumed as

$$\forall i,j,\;{\text{Corr}}[{\text{Cost}}_{m_1 }^i ,{\text{QALY}}_{m_2 }^j ] = \left\{ {\begin{array}{ll} {\rho_{QC}^i ,} &\quad {i = j,\;m_1 = m_2 } \\ {\rho_{QC}^{ii} = 0,} &\quad {i = j,\;m_1 \ne m_2 } \\ {\rho_{QC}^{ij} = 0,} &\quad {i \ne j.} \\ \end{array} } \right.$$

(22)

The structure of each correlation is schematically shown in the following Fig. 1.

### Simulation methods and settings

Given a selected intervention structure, the numerical simulation uses Monte Carlo sampling to form a column vector that combines the costs and QALY outputs of each patient under the assumed parameters. A sufficient number of samples are generated and the statistical features of the overall results (such as population benefits and total costs) are calculated. Then, the results of the objective function under the given intervention structure are obtained. Following this, different stochastic CEA methods are applied to this problem to get their own judgements on optimal interventions. Performance of population benefits are compared with these optimal interventions. Moreover, simulations are performed again for each method’s recommended intervention structure modification, showing the effect of recommendations for different stochastic CEA methods’ judgements. Some key intermediate results are also recorded and presented.

#### Simulation process

The purpose of simulation mainly lies in three aspects. Firstly, illustrating through numerical examples that, in terms of population benefits, the mixture of different interventions could have a better return-to-risk result than an absolutely unique intervention. Secondly, under different intervention structures, calculating the risk-adjusted expected cost and risk-adjusted ICER proposed in this paper, and verifying the process of inducing the optimal risk-return ratio of population benefits in a marginal sense. Thirdly, comparing the decisions generated by different stochastic CEA methods for the same problem, testing the performance of different methods on population benefits, and checking their suggested direction of optimization.

Specifically, we first select two interventions from all optional ones, separately calculate their risk-return levels of population benefits, and then mixes them in different proportions to form intervention structures. Secondly, under different intervention structures, the risk-adjusted expected cost and risk-adjusted ICER of the two interventions are calculated, along with the corresponding decisions that the criterion makes for optimizing the target function. Then, to show the effectiveness of each method, a specific intervention structure is selected, and different stochastic CEAs are used for making determinations between any two interventions. The changes in the performance of population benefits are measured after adjusting the intervention structure according to the corresponding determination. These changes can show the effectiveness of different methods.

In the simulation process, it is necessary to simulate the QALY outcomes and costs under the selected intervention and evaluate the population benefits quantitatively. In detail, under each selected intervention structure M, the evaluation of the return-to-risk ratio of the population benefits level consists of the following steps in the simulation process:

1. 1)

Divide the M simulated patients into N groups, where the ith group has Mi patients according to the structure of M.

2. 2)

Combine the random variables of cost and QALY output of each patient into a random vector and calculate the mean and covariance matrix of the random vector.

3. 3)

Generate 10,000 samples of the group using the calculated mean and covariance matrix as the parameters of a high-dimensional multinormal distribution.

4. 4)

Measure the population benefits of each sample of the group.

5. 5)

Summarize the population benefits of each sample and measure their statistical characteristics and return-to-risk ratios.

Since the simulation can only generate a limited number of samples, selecting extremely tailed risk measures may cause the simulation results to be unstable. Therefore, we repeat the evaluating steps mentioned above 100 times independently and use the average value as the final result. Finally, we use empirical probability of correctly recommendations to estimate the accuracy of different CEAs.

#### Parameter settings and reference methods selection

To balance the adequacy of alternative interventions and the interpretability of results, we select a total of N = 5 interventions and assume that the patient population consists of M = 100 individuals. The number of samplings of the outcome results of any patient subgroup is chosen to be 10,000, and the number of repeated experiments for the return-to-risk evaluation of population benefits is chosen to be 100 times.

The assumed distribution parameters within intervention subgroups are shown in Table 1:

In addition, the correlation coefficients between different individual outcomes among different groups are set to $$\forall i \ne j,\;\rho_C^{ij} = 0.005,\;\rho_Q^{ij} = 0.002$$, and the monetary constant, $$\lambda$$, is set to $$\lambda = 2000\;\ /{\text{QALY}}$$ unless specifically mentioned.

Besides the method proposed in this paper, we select four different stochastic CEA methods as reference methods to participate in the evaluation together. The first is the classical ICER method, where the criterion measure is

$${\text{ICER}}_{i,j} ({{\mathbf{M}}}) = \frac{{{{\mathbb{E}}}[{\text{Cost}}^j ] – {{\mathbb{E}}}[{\text{Cost}}^i ]}}{{{{\mathbb{E}}}[{\text{QALY}}^j ] – {{\mathbb{E}}}[{\text{QALY}}^i ]}}.$$

(23)

The second is the Sortino ratio method (see Elbasha ) based on NMB, which compares the Sortino ratios (a kind of return-to-risk ratio using the expected benefit dividend by the downside variance) of different interventions to determine their priority. The third and fourth methods are the probabilistic methods based on CEAC, which determine the superiority of the interventions by judging whether the α percentile of the ICER distribution, $$Q_\alpha (\frac{{{\text{Cost}}^j – {\text{Cost}}^i }}{{{\text{QALY}}^j – {\text{QALY}}^i }})$$, satisfies the requirement of the willing-to-pay level ($$\lambda$$). The third method is relatively aggressive, using α = 5%, while the fourth method is relatively conservative, using α = 95%. Comparisons are made between the method proposed in this paper and these four reference methods.

### Simulation results

The simulation results are mainly presented in three parts. The first part shows the superiority of mixed intervention compared to unique intervention from the perspective of population benefits. The second part displays the relationship between changes in the risk-adjusted ICER and the intervention structures, which demonstrates the marginal optimality under different intervention structures. The third part shows the preference decisions of different interventions by different stochastic CEA methods under specific intervention structures, which demonstrates the validity of the decision by the tool of this paper.

#### Population benefits with mixed intervention

In this subsection, we evaluate the population benefits under different intervention structures through simulation, and the incremental population benefits introduced by intervention mixing. It is verified that mixed interventions can lead to better return-to-risk ratios. Specifically, intervention 4 and intervention 5 are chosen as the controlled subgroups of this part, while the first three interventions have proportions of 0 in the structure, i.e., M1, M2, M3 = 0 to avoid extra interference. Then, the level of M4/M gradually varies from 0 to 1, and the corresponding return-to-risk ratios of the population benefits are calculated. The results are shown in Fig. 2.

From Fig. 2, it can be seen that as the proportion of intervention 4 gradually increases from 0 to 1, the return-of-risk ratios of population benefits shows a concave pattern, which firstly increasing and then decreasing. The target function reaches its optimum at around M4/M = 0.4. The optimum level, 0.0227, is significantly higher than the return-of-risk ratio of 0.0173 at M4/M = 0 and 0.0127 at M4/M = 1. This suggests that if considering the medical decisions with intervention 4 and intervention 5, using a 4:6 structure to mix both interventions is more efficient in a population view than using either intervention alone. This validates that it can be superior to mix interventions for population benefits.

Following the previous subsection, we calculate risk-adjusted expected cost of both interventions and the risk-adjusted ICER for switching from intervention 4 to intervention 5 marginally as we gradually adjust the proportion of M4/M in the [0, 1] range. Through the simulation results, it is shown that the risk-adjusted ICER changes can lead the intervention structure to achieve the optimum. The simulation results are shown in Fig. 3.

In Fig. 3, the red curve represents the risk-adjusted expected cost of intervention 4, and the blue curve represents that of intervention 5. It can be seen that as the proportion of intervention 4 increases, the risk-adjusted expected costs’ difference between intervention 5 and intervention 4 changes from extremely positive to extremely negative. This is caused by changes in the diversification structure. Specifically, when the proportion of intervention 4 is small, there are many cases using intervention 5 in the population. At this time, increasing one unit of intervention 5 cases is equivalent to adding a highly positive correlated risk to the total risk pool, resulting in a relatively high level of risk allocated to the newly added one. This makes the risk-adjusted expected cost of the new unit of intervention 5 high at this situation. In contrast, the risk-adjusted expected cost of intervention 4 is relatively low because it is relatively independent with the risk pool. This pattern validates that the risk-adjusted expected cost does take the marginal impact of intervention restructuring on the risk pool into account in the risk adjustment.

In addition, the black curve in Fig. 3 represents the evaluation results of the risk-adjusted ICER from intervention 4 to 5. It can be seen that this curve intersects with the $$\lambda = 2000$$ threshold at around M4/M = 0.4, the optimal structure point. To the left of the intersection point, the risk-adjusted ICER is higher than the monetized constant, indicating that switching one unit of patients from intervention 4 to intervention 5 is not cost-effective and has a negative impact on the population benefits. Therefore, the proportion of intervention 5 needs to be reduced. On the contrary, to the right of the intersection point, the opposite is true. It is worth noting that the intersection of the black curve in Fig. 3 with $$\lambda = 2000$$ coincides perfectly with the position where the return-of-risk ratio of population benefits gets its optimum in Fig. 2. This also confirms that the risk-adjusted ICER can induce changes in the intervention structure to achieve the optimum.

#### Comparison of methods

In this subsection, we selected a specific intervention structure (equal-weighted) as the starting point to examine the decisions of different stochastic CEAs for adjusting interventions and quantified the corresponding impact on the distribution of the population benefits. These results are shown in Fig. 4.

“Traditional ICER” represents the results of the traditional ICER tool, “Risk-adjusted ICER” represents our proposed method, “Sortino Ratio” represents the decision made using the Sortino ratio tool, “CEAC at 5%” represents the aggressive CEAC method, and “CEAC at 95%” represents the conservative CEAC method.

In Fig. 4, each sub-plot presents the information in a 5 × 5 matrix, where the grid in the ith row and jth column of any sub-plot represents the decision made by the corresponding stochastic CEA tools or its correctness, for changing an unit of patients from intervention i to intervention j. The decisions are shown in the sub-plots of the first row, where + 1 represents a decision to accept the change, − 1 represents a decision to accept an inverse change, and 0 represents no modification. The correctness of the decision results is displayed in the sub-plots of the second row, where the number in each grid represents the percentage of scenarios where population benefits are improved after the changes are taken based on the corresponding decision, within the 100 times of simulated scenarios. Additionally, the cells on the diagonal of each sub-plot do not contain any information.

From this comparison, it can be seen that the decisions of the traditional ICER and CEAC methods differ from our proposed method and the Sortino method to some extent. Specifically, the traditional ICER does not consider the stochasticity and only requires the ICER to be below the monetized constant in an expected view for each accepted intervention change. And, the CEAC method, similar to the traditional ICER, does not use the information contained in the current intervention structure. The Sortino Ratio method obtained similar decisions with respect to our proposed method, which can be attributed to the similarity in risk-adjusted measures between the two methods. However, the Sortino Ratio method still does not consider the current intervention structure, so there are still differences. These differences are not very significant due to the low level of correlation assumed in the simulation.

Also, it can be found that only our proposed method can almost perfectly ensure that the decisions to achieve an accuracy over 50%, which also reaches the highest level among all methods in expectation (The absolute value of the accuracy is related to the number of simulated samples, so the relative level of the accuracy of each method is mainly compared.). In addition, except for the judgementFootnote 1 between intervention 2 and 1, our proposed method achieved the highest accuracy rate among all methods for any pair of comparisons. This also validates the effectiveness of our proposed method.

### An empirical example

In this subsection, we discuss an empirical example using the information extracted from recent research on comparing direct oral anticoagulants compared to low-molecular-weight-heparins for treatment of cancer associated venous thromboembolism, see Muñoz et al. . It is worth noting that we just use some of the parameters from the research, but do not try to answer or challenge any problem or conclusion about this health technology assessing problem in real world. For short, we just borrow some parameters from this work to show that the population can benefit from mixed interventions, and we can find the optimal mix structure with our tools.

There are four kinds of interventions introduced in this article, the low molecular weight heparin (LMWH), the Apixaban, the Rivaroxaban and the Edoxaban. All these four interventions could be used for treatment of cancer associated venous thromboembolism. The paper shows the expectation of costs and QALYs (12 months) of all these four interventions, and we can roughly extract the coefficient of variation (CV, standard error divided by mean) form the confidence intervals and simulation figures. With some more assumptions, like $$\lambda = 30000$$, the parameters in our example are shown in Table 2.

To show the benefit from mixed interventions and the optimality of our proposed structure, we compare six situations. Four of them are using these four interventions alone, the fifth is an equal weighting structure and the final is the structure generated from our tool. The mean, variance and the return-to-risk ratio are calculated for all these six situations, and the results are shown in Table 3.

We can find in Table 3 that none of the four interventions alone could make the return-to-risk ratio of population benefit more than 0.100, but just a simple mix with equal weight can make it 0.112, which clearly shows the benefit from mixing interventions. Moreover, simple mix does benefit, but only the optimal structure generated by our tool reaches the highest ratio of 0.148. This is around 0.7 times higher than solely use the intervention suggested by classical CEA tools.

The optimal structure is {0.00%, 32.62%, 33.22%, 34.16%}, which gives zero weight on LMWH and nearly equal weights on Apixaban, Rivaroxaban and Edoxaban. This consistent with the judgment of the original article that the intervention of LMWH is dominated. This also shows the reliably of our tool. 