Authors: Giuliano Netto Flores Cruz, Keegan Korthauer
Categories: Research Article, Bayesian, R package, clinical decision‐making, clinical prediction models, decision curve analysis, diagnostic tests
Source: Statistics in Medicine
Doi: 10.1002/sim.10277
Authors: Giuliano Netto Flores Cruz, Keegan Korthauer
Clinical decisions are often guided by clinical prediction models or diagnostic tests. Decision curve analysis (DCA) combines classical assessment of predictive performance with the consequences of using these strategies for clinical decision‐making. In DCA, the best decision strategy is the one that maximizes the net the net number of true positives (or negatives) provided by a given strategy. Here, we employ Bayesian approaches to DCA, addressing four fundamental concerns when evaluating clinical decision (i) which strategies are clinically useful, (ii) what is the best available decision strategy, (iii) which of two competing strategies is better, and (iv) what is the expected net benefit loss associated with the current level of uncertainty. While often consistent with frequentist point estimates, fully Bayesian DCA allows for an intuitive probabilistic interpretation framework and the incorporation of prior evidence. We evaluate the methods using simulation and provide a comprehensive case study. Software implementation is available in the bayesDCA R package. Ultimately, the Bayesian DCA workflow may help clinicians and health policymakers adopt better‐informed decisions.
In Decision Curve Analysis (DCA), we are typically interested in estimating the net benefit of adopting a given clinical decision strategy [1]. A decision strategy may be simply treating all patients under suspicion of a given condition because the condition is so deadly that any costs of a potentially unnecessary intervention are outweighed by devastating costs of neglecting necessary treatment, even if the risk is small—e.g., treating potentially aggressive cancer. In this “Treat all” strategy, intervention happens in all patients regardless of the true underlying disease status—e.g., whether the patient's cancer is aggressive or not. Conversely, an intervention may be high risk while the underlying condition is of favourable prognosis—e.g., surgical removal of a stable noncancerous brain tumour. In this case, a reasonable decision strategy is to not treat any patient—the “Treat none” strategy.
In general, the relative risk conferred by the disease and the treatment does not always clearly side with either the “Treat all” or “Treat none” strategies. In addition, there is often uncertainty around a patient's disease status or prognosis. In this case, a decision strategy could be based on a prediction model that estimates, e.g., a patient's likelihood of having aggressive disease. If the patient's likelihood is above a decision threshold t, then we intervene. Beyond the probability of having a disease right now (diagnostic setting), the threshold t could also be the probability of a future event like death, hospitalization, or disease progression (prognostic setting). The same idea serves for binary tests, in which case we intervene if the test is positive.
In the context of DCA, the Net Benefit (NB) at the decision threshold t can be written as [1]:
(1) NBt=(TPt−FPt·wt)/n
where TPt and FPt are corresponding true and false positive counts, n is the total sample size, and wt=t/(1−t). Given a decision threshold t, this definition (1) fixes the weight of each true positive at 1, which mathematically implies a relative weight of wt for each false positive. This allows decision analysis without the need to specify the absolute costs of each potential outcome (true and false positive/negative). Instead, we rely on a clinically‐motivated decision threshold t, which properly weights true and false positives/negatives based on the clinical context [2]. To consider a range of relative weights (e.g., due to disagreement between clinicians or even patients' preferences), the decision curve is then constructed by plotting NBt for a reasonable range of decision thresholds t. At each decision threshold, the best decision strategy is the one that maximizes the net benefit [1]. Definition (1) can also be further refined by deriving strategy‐specific costs (known as “test harm” terms).
When interpreting DCA results, there may be special considerations for making a decision that changes a well‐accepted practice. Beyond the one‐time cost of the implementation process itself, there is always the risk of implementing the “wrong” strategy whose apparent optimality was an artifact of chance. As with any other estimate, the observed net benefit is subject to random variation in the data. Addressing this uncertainty, however, depends on how we understand risk. Under risk neutrality, we only care about expected gains (or costs), and uncertainty quantification does not change which strategy should be used, although it can motivate further research [3, 4, 5]. If the strategy with the highest observed net benefit indeed has the highest true net benefit, then foregoing implementation due to uncertainty incurs harm as we will be sticking with a worse strategy until more data are available. On the other hand, prematurely replacing the current Standard of Care (SoC) with a new decision strategy poses potentially irrecoverable costs to individual patients, healthcare institutions, and even healthcare systems [5, 6, 7]. This setting motivates risk aversion, which may require a more careful assessment of uncertainty to prevent premature implementation.
The work presented here regards the Bayesian estimation of net benefit and is a priori indifferent as to whether end users are risk‐neutral or not—a debate we do not aim to settle. In effect, end users can decide how to interpret uncertainty in the net benefit estimates. To describe the full potential of the method, we include a case study that assumes risk aversion when assessing a new decision strategy to potentially replace a hypothetical SoC. In this case, if we are too uncertain about the net benefit gain from adopting the new decision strategy, then more data may be desirable before changing current practice. Given the superiority of a decision strategy in terms of observed net benefit, different levels of uncertainty may be compatible with context‐specific costs to ultimately justify implementation. Still, the risk‐neutral reader is free to interpret uncertainty as a motivation for further research only, without changing the conclusions based on the observed estimates, if desired.
Therefore, uncertainty quantification around estimates of net benefit and their differences allows two potential interpretations, depending on the reader's risk profile. Under risk neutrality, uncertainty informs whether more research is needed but does not change our decision based on the currently‐available point estimates. Under risk aversion, uncertainty may inform whether the available data provides enough evidence to change well‐established clinical practice. Both interpretations of uncertainty in DCA motivate further research, but only the risk‐averse wait for more data before changing clinical practice currently in place. In what follows, we employ Bayesian approaches to DCA, addressing four fundamental concerns when evaluating clinical decision Which strategies are clinically useful?What is the best available decision strategy?Which of two competing strategies is better?What is the expected net benefit loss associated with the current level of uncertainty?
Uncertainty around these concerns is natural to address in the Bayesian context; however, the approach remains agnostic to its interpretation. While often consistent with the frequentist approach in terms of point estimates (e.g., under vague priors), the fully Bayesian estimation of net benefit allows for an intuitive probabilistic interpretation of DCA results as well as for the principled incorporation of prior evidence. Our proposal builds on the work from Wynants et al. [8] which proposed Bayesian DCA for evidence synthesis in meta‐analysis of data from multiple settings (e.g., multiple hospitals). We adapt their binary outcome model to allow for the incorporation of prior information in the single‐setting case and propose an alternative formulation for survival outcomes. We then compare the methodology with Frequentist alternatives using simulation and provide a case study with openly available data. The proposed approaches are implemented in the freely available bayesDCA R package.
Suppose we have access to validation data to assess one or more clinical decision strategies, be it a pre‐specified clinical prediction model or a binary diagnostic/prognostic test. We will use DCA to decide whether any of the strategies under investigation are clinically useful and which of them is the best. We will also examine the difference between strategies, which may be useful to evaluate alternatives under scenarios where one or more strategies are unavailable. Finally, we would like to have a sense of whether the current study is precise enough so that new studies with the same population are not necessary. In what follows, we first describe the Bayesian estimation of decision curves for binary outcomes and then extend it to survival outcomes.
The Net Benefit formulation in (1) can be rewritten in terms of the outcome prevalence (p) and the threshold‐specific Sensitivity (Set) and Specificity (Spt). For a decision threshold t:
(2) NBt=Set·p−(1−Spt)·(1−p)·wt
During (external) validation of predictive models or binary tests, we can estimate the parameters above using a conjugate Beta‐Bernoulli joint model for the indicator variables of positive predictions and disease status (see Section 9.2 for full model specification). Given conjugacy, the full posterior distribution of the parameter vector p,Set,Spt is known in closed
(3) Beta(p|D+α0,ND+β0)×Beta(Set|TPt+α1,FNt+β1)×Beta(Spt|TNt+α2,FPt+β2)
where TPt,FPt,TNt,FNt,D,andND represent the total number of true and false positives, true and false negatives, and individuals with and without the disease, respectively. The (α·,β·) terms are parameters of the independent Beta prior distributions. Within bayesDCA, we set α·=β·=1 as a default, representing uniform priors on the (0,1) interval, though the user may choose different priors as well. We suggest a more informative prior based on the expected relationship between the decision threshold and sensitivity/specificity in the Supporting information (Section Informative priors in Bayesian DCA for binary outcomes and Figures S1 and S2). In addition to conjugacy, the factorization in (3) is due to parameter orthogonality in the likelihood function and implies posterior independence. This means that Markov‐Chain Monte Carlo (MCMC) is not we can combine samples from the marginal posteriors to easily generate valid samples from the joint posterior, making estimation particularly fast—typically a fraction of a second for an entire DCA.
The model (3) can be seen as a single‐setting version of the model from Wynants et al. [8], which models the parameters of interest on the logit scale using a multivariate normal (MVN) distribution. Given the meta‐analysis context, there are multiple sensitivities, specificities, and prevalences to be considered at each threshold; their MVN formulation accounts for correlation across these parameters in a random‐effects fashion. Here, however, there is only one triple of sensitivity, specificity, and prevalence at each threshold, so there is no correlation estimate to be modeled. In fact, jointly modeling positive predictions and disease outcomes makes the parameters orthogonal by construction, and, hence, their estimators are also independent (see Section 9.2). On the other hand, the model from Wynants et al. [8] benefits from partial pooling, being more appropriate for evidence synthesis in the multiple‐setting context (e.g., multiple validation samples from different hospitals or a meta‐analysis of multiple studies). The random effects formulation is a natural choice to handle data from multiple settings due to the expected heterogeneity in the net benefit, which arises from variation in setting‐specific prevalences, sensitivities, and specificities [8]. Finally, unlike sampling from (3), Wynants et al. [8] method requires MCMC and is, therefore, expected to be considerably slower.
We use data from the GUSTO‐I trial [9], a large randomized study involving thrombolytic treatments for Acute Myocardial Infarction, as an example. Figure 1 shows the resulting Bayesian DCA with the Frequentist counterpart superimposed. While closed‐form asymptotic versions also exist [10, 11, 12], throughout we will employ the bootstrap approach as the Frequentist version of DCA as it is not restricted to the large sample size setting. Although point estimates mostly coincide, notice that the bootstrap‐based DCA collapses to zero as the threshold increases—in particular, as the threshold approaches the maximum observed risk prediction. The same does not happen with the Bayesian approach, which continues to naturally propagate uncertainty through the net benefit equation even in the absence of events, though with an increasing influence of the prior distribution. This behavior can be useful, e.g., in settings with small effective sample sizes, such as when most risk predictions are concentrated on one side of the decision threshold. While a decision threshold higher than all risk predictions in the population does imply zero net benefit, in small samples this may happen by chance alone—in which case a positive net benefit would require more informative priors. Nonetheless, the bayesDCA R package warns the user if no events were observed above some decision threshold. Notice that, at very high thresholds, high sensitivity and specificity are required to yield a positive net benefit because of the increasing weight of false positives (i.e., large wt). Thus, unless under strong prior knowledge or in the face of substantial evidence of high sensitivity and specificity, the posterior distribution will tend toward net clinical harm in this region. The Bayesian DCA, therefore, indicates that harm is likely at higher thresholds for this example using vague priors.
![FIGURE 1: Bayesian DCA captures uncertainty across the entire decision curve. Bayesian DCA was computed using the bayesDCA R package, while the Frequentist alternative used the bootstrap‐based rmda package. An example model was built using data from the GUSTO‐I trial [9], while DCA was constructed with held‐out data (N=500, 36 events). The bootstrap intervals and point estimates collapse to zero as the threshold approaches the maximum risk prediction. The Treat all curve is nearly identical for both approaches, and the Bayesian version is shown. Intervals correspond to 95% confidence and credible intervals for Frequentist and Bayesian methods, respectively. Notice that overlapping uncertainty intervals do not mean there is no difference between decision curves due to correlation in the sampling distribution of net benefit estimators.](SIM-43-6042-g007.jpg)
Beyond uncertainty propagation, three key advantages come naturally with the Bayesian approach for binary outcomes implemented in bayesDCA. First, one can easily use full external information to estimate the prevalence parameter. In case‐control studies ^1^ , DCA needs to be adjusted for the population prevalence, which is usually done by plugging in a point estimate and completely ignoring uncertainty [10]. Here, we may simply sample p from the posterior (3) using data from an external cross‐sectional study or the source population cohort. Equivalently, this could be seen as constructing an informative prior distribution from external prevalence information and sampling from the prior. Regardless of the source, the prevalence posterior distribution is then used to compute the net benefit. This takes into account the uncertainty in the prevalence parameter because we are using the raw prevalence data instead of plugging in a point estimate. The second advantage that our Bayesian approach brings is the option to use informative priors to improve estimation—in terms of both uncertainty and point estimates. For instance, we know that Se is higher for small thresholds and lower for large thresholds—and that the reverse is true for Sp. The informative prior suggested in the Supporting information takes advantage of this reasoning (Section Informative priors in Bayesian DCA for binary outcomes and Figures S1 and S2). Additionally, prior elicitation is straightforward for binary tests since their sensitivity and specificity are fixed across thresholds [14].
The third and most notable advantage of our Bayesian approach for DCA is the ability to interrogate posterior decision curves with an intuitive probabilistic interpretation. Since we have access to the full posterior distribution of the net benefit across all thresholds of interest, we can compute arbitrary functions to help us interpret the DCA output probabilistically. Since this advantage applies to Bayesian DCA in general and not just to binary outcomes, we first propose a method of Bayesian DCA for survival outcomes and then describe the proposed probabilistic interpretation framework based on the interrogation of the posterior decision curves.
Many decision strategies address prognostic problems. For such survival outcomes, we must rewrite the net benefit formula to be able to account for censoring as
(4) NBtτ=[1−Sτ|r^τ>t]·Pr^τ>t−Sτ|r^τ>t·Pr^τ>t·wt
where r^τ is the predicted risk of the event of interest at time τ (a probability in the case of a prognostic model and 0 or 1 in the case of a prognostic test). At time τ and threshold t, the probability of a positive prediction is given by Pr^τ>t, and Sτ|r^τ>t is the survival probability given a positive prediction.
To estimate (4), we jointly model survival times T (with censoring indicators C) and the indicator of positive predictions Z with Weibull and Bernoulli likelihoods, respectively (see Section 9.2 for full model specification). Under independent priors, the posterior distribution factorizes due to parameter orthogonality
(5) π(p,θ1|Data)∝π(p|𝒟0)×π(θ1|𝒟+)
where 𝒟0={Zi}i=1n is the set of positive prediction indicators, 𝒟+={Ti,Ci}i∈[n]:zi=1 is the survival dataset for patients with positive predictions. Here, p=Pr^τ>t while θ1=(α1,σ1) represents Weibull shape and scale parameters, respectively. Although the resulting posterior distribution does not have a closed form, it does benefit from orthogonality between the Weibull and the Bernoulli components, allowing us to estimate them separately. Hence, we put a Beta prior on p (uniform by default, as before) to take advantage of conjugacy, while θ1 is estimated with MCMC using Stan [15]. In bayesDCA, the default priors for the Weibull parameters are
(6) α1∼Half‐Student‐t(5,0,1.5)σ1∼Half‐Student‐t(30,0,100)
These priors put a nearly equal prior probability on increasing and decreasing hazards and are largely vague with respect to scale. bayesDCA also allows user specification of the prior parameters above and provides a Gamma prior option instead of Half‐Student‐t. Once sampling is done, we then combine the draws from the posterior distributions of θ1 and p to compute the net benefit given in (4).
We can now interrogate the posterior decision curves for all decision strategies under investigation. The entire interpretation framework proposed in the next section is immediately available for both survival and binary outcomes, highlighting once again the advantages of the proposed Bayesian approach.
The main advantage of Bayesian DCA is the ability to arbitrarily interrogate posterior distributions of decision curves. This allows for a probabilistic interpretation framework that helps understand the degree of uncertainty imposed by the currently available data on the observed decision curves. We may ask, for what is the posterior probability that the model under investigation is useful at a given threshold? Following the definition in Wynants et al. [8], that
(7) P(useful)=P(NBmodel>max{NBtreat all,NBtreat none})
where NBtreat none is always zero. However, if we have two or more competing models, another natural question may what is the probability that my model is the best decision strategy available? For instance, for a given “model 1”:
(8) P(best)=P(NBmodel1>max{NBtreatall,NBtreatnone,NBmodel2,NBmodel3,⋯})
Notice that Wynants et al. [8] employ the term P(useful) for both Equations (7) and (8), whereas here we define P(useful) and P(best) separately. This is because the best decision strategy might not be available everywhere, so the usefulness of remainder strategies is still relevant in that case. The above definitions can be extended to an arbitrary number of models or tests and be computed across all decision thresholds. Within bayesDCA, one may also compute pairwise comparisons between, say, two models with very similar net benefits. For instance, compute the probability that a given strategy beats another by at least c net benefit units (i.e., net true positives):
(9) P(NBmodel1−NBmodel2>c)c≥0
Although we would take c=0 in general, c could be positive if considering, e.g., different test harms between model 1 and model 2. Upon full posterior interrogation, we may reach a better understanding of the net benefit profiles of the decision strategies under if there is large uncertainty around the decision curves, the above probabilities will be inconclusive (i.e., far away from 0% and 100%), so we may opt to collect more data before making a decision.
One way to directly quantify the expected consequences of the current level of uncertainty is to compute the Expected Value of Perfect Information (EVPI) for model validation [16]:
(10) EVPI=Emax−maxEEmax=E[max{NBtreat all,NBtreat none,NBmodel1,NBmodel2,⋯}] maxE=max{E[NBtreat all],E[NBtreat none],E[NBmodel1],E[NBmodel2],⋯}
where the maximum in Emax is computed for each draw of the joint posterior distribution, whereas maxE is a maximum of posterior means. The EVPI may be seen as the expected net benefit loss due to the current level of uncertainty in the decision curves. For instance, if the EVPI is 0.1, picking the observed best decision strategy is associated with an expected loss of 0.1 net true positives as compared to picking the actual best decision strategy.
It is important to notice another advantage of the parametric Bayesian approach to computing the EVPI. As a sample statistic, the EVPI is expected to decrease monotonically with the sample size [16]. In a simulation study using the GUSTO‐I trial data, however, Sadatsafavi et al. [16] showed that small effective sample sizes may cause the observed EVPI to escape this monotonic behaviour, especially at very low or very high decision thresholds. One way to avoid this issue is to employ informative prior distributions in Bayesian DCA. Reproducing the simulation code from Sadatsafavi et al. [16], we show how our Bayesian approach can recover the expected monotonic behaviour of the EVPI in the Supporting information (Section Informative priors preserve EVPI monotonic behaviour and Figure S3).
In the next section, we provide a simulation study of the empirical performance of the suggested Bayesian approaches. In the following section, we present a case study to highlight the bayesDCA workflow, where we fully interrogate the posterior decision curves to quantify uncertainty around the answers to the four fundamental questions mentioned above. Here we focus on binary outcomes, though the same workflow for survival outcomes is also implemented and easily accessible through the bayesDCA R package.
To test the approach for binary outcomes, we simulate a population with an underlying logistic regression model as the data‐generating process and select an example model to be evaluated using DCA. To represent a range of scenarios, we vary the outcome prevalence and the maximum achievable discrimination, measured by AUC, representing the setting's signal‐to‐noise ratio. To resemble a common scenario of overfitting, the example model is its predictions are overly extreme (i.e., too close to zero or to one). Each simulation run emulates a different external validation or test dataset selected at random from the setting's population, with which we perform DCA. The sample size for each simulated dataset is set so that the expected number of events is 100. While potentially oversimplified, these simulation scenarios are designed to cover a sufficiently broad range of net benefit values. See Section 9.3 for a full description of the simulations. Average performance statistics from the fixed models for each setting are provided in Table S1.
For each setting, we ran 1,000 simulations performing both Bayesian DCA and Frequentist DCA for comparison. Figure 2 shows the resulting distributions of estimation errors. As expected, the distributions of point estimate errors are nearly identical for the Bayesian and the Frequentist approaches in almost all simulation settings and decision thresholds. This is in general, with 100 expected events, any influence of the default vague priors is unnoticeable. One exception is the setting with an AUC of 0.65 and a prevalence of 1%, in particular at a threshold of 75%. Given the low prevalence and discrimination, very high thresholds such as 75% often yield no positive prediction (i.e., no risk prediction above the decision threshold), so the effective sample size for that threshold is low and the prior matters more. This is the setting in which bootstrap‐based point estimates and intervals are expected to collapse to zero, while the Bayesian estimates are regularized by the prior distribution. Still, the discrepancy between the two approaches is negligible in absolute terms—see Figure S4 for a plot of the point estimates in the original scale with the true net benefit overlayed.

We also assessed the empirical coverage of 95% uncertainty intervals for both Bayesian and bootstrap‐based Frequentist methods. Although not required for a valid Bayesian analysis, users of the bayesDCA R package may feel more comfortable having Frequentist calibration of credible intervals (Cr.I.), at least under vague priors. As shown in Figure 3, the empirical coverage of both the Bayesian and Frequentist decision curves closely matches the nominal value of 95% for the entire range of decision thresholds in most simulation settings. The only exception is, once again, the setting with an AUC of 0.65 and a prevalence of 1%. At thresholds of 10% or above, the bootstrap intervals show increasing undercoverage, reaching an empirical coverage of nearly 0% for the 75% threshold. This happens because although the true net benefit is not exactly zero in this setting due to the presence of positive predictions in the population, the bootstrap intervals tend to collapse to zero due to the lack of positive predictions in the observed samples. In contrast, under the same problematic scenario, Bayesian intervals show only a slight overcoverage. Despite its vagueness, the default prior distribution yields valid, credible intervals with reasonable Frequentist calibration. Moreover, the Bayesian credible intervals are not significantly wider than the Frequentist confidence intervals, except when the Frequentist method fails. Overall, Figure S5 shows that the Bayesian intervals have reasonable width, even when the Frequentist counterpart collapses to zero.

Finally, as shown in Figure S6, the Bayesian DCA for binary outcomes is orders of magnitude faster than its bootstrap‐based Frequentist counterpart. Here, we are sampling 4000 draws from the posterior distribution (3) and using 500 bootstrap samples—the default in bayesDCA and rmda, respectively. Moreover, computation time significantly increases with the overall sample size for the bootstrap case, but not in the Bayesian case. Given the 100 expected events fixed for each simulation run, simulation settings with prevalences of 30%, 5%, and 1% imply overall sample sizes around 333, 2,000, and 10,000, respectively. While this is expected to impact bootstrap speed, the computation time for the proposed Bayesian DCA is virtually unaffected because (3) only depends on simple summary statistics. In our experience, running Bayesian DCA for binary outcomes with bayesDCA takes no more than a second in most cases (using a standard laptop with 12 GB of RAM and no parallelization).
We follow the same simulation strategy as above. The underlying populations follow Weibull distributions with covariates satisfying the proportional hazards (PH) assumption, while censoring times are uniformly distributed. The decision strategies being investigated are based on fixed PH Cox models with exaggerated coefficients—i.e., miscalibrated due to overfitting. For each setting, we simulate 1000 datasets, with which we perform both Bayesian and Frequentist DCA using a 12‐month prediction horizon. We vary the underlying C‐statistics and 1‐year survival rate. For brevity, we show the results for C‐statistics 0.6 and 0.9 with a 1‐year survival of 10%. A full description of the simulations is provided in the Section 9.3, and further simulation settings are reported in the Figures S7–S9. Average performance statistics from the fixed models for each setting are provided in Table S2.
As shown in Figure 4, the Bayesian point estimates and uncertainty intervals generally behave similarly to the ones from the Frequentist approach. At decision thresholds of 50% or above, the effective sample size can be small, which causes a slight bias for both approaches depending on the simulation setting—here seen at the 75% threshold for the Frequentist approach under C‐statistic 0.9, see also Figure S7. The Bayesian uncertainty intervals show reasonable empirical coverage overall, though we observed undercoverage for very high or very low thresholds depending on the simulation setting. While the setting with C‐statistic 0.9 shown in Figure 4 represents the worst‐case scenario, most coverage probabilities remained above 90% across all settings—Figure 83. Mean absolute percentage errors of point estimates were comparable between Bayesian and Frequentist approaches across all simulation settings—Figure S9. The Bayesian approach was significantly faster than the bootstrap‐based alternative—Figure S10. However, because our bootstrap implementation was not optimized for speed, this should be seen as a worst‐case scenario for the Frequentist case. Finally, the dcurves R package may fail to produce net benefit estimates in some simulation settings for higher thresholds—Figure S11. This is because the default Kaplan‐Meier method implemented in the package (likely intentionally) does not produce survival estimates at time points beyond observed data.

In summary, our Bayesian approach offers an accurate alternative for estimating decision curves for both binary and survival outcomes. It also enables uncertainty quantification for survival outcomes and provides much faster quantification of uncertainty for binary outcomes. These methods are implemented and easily accessible in the bayesDCA R package. While the simulations shown in this section apply the default weakly informative priors, different priors may be specified by the user if desired to further improve estimation performance. The bayesDCA R package allows sampling from the prior only so that users can easily visualize the consequence of the specified priors on the net benefit scale—see Figures S1 and S2.
Clinical prediction models are commonly employed to predict cancer diagnosis. For example, the ADNEX model predicts an individual's risk of having ovarian cancer with high discrimination (AUC>0.9) and adequate calibration [17]. The model employs clinical and ultrasound features from patients under suspicion of ovarian cancer due to the known or suspected presence of adnexal masses. Previously, Wynants et al. (2019) warned about the importance of utility‐based decision thresholds and used the ADNEX model as a motivating example [2]. While in practice thresholds from 10% to 40% are used for ADNEX [18], the authors suggest that a reasonable decision threshold may be as low as 6% due to the high cost of false negatives, which can cause late detection and treatment of aggressive cancer. Here, we will expand on their example using Bayesian DCA in a hypothetical scenario.
Suppose we perform an external validation study to assess the predictive performance of the ADNEX model. We wish to know whether the model should replace the Standard of Care (SoC) currently in a hypothetical diagnostic test with 81% sensitivity and 88% specificity. Using the publicly‐available data from Wynants et al. [2] (N=2403, 980 events), the Bayesian DCA results are shown in Figure 5.
![FIGURE 5: Illustration of Bayesian DCA for the ADNEX model and hypothetical Standard of Care diagnostic test. Bayesian DCA was computed using the bayesDCA R package and publicly available data from Wynants et al. [2] (N=2403, 980 events). A hypothetical Standard of Care diagnostic test was simulated to have 81% sensitivity and 88% specificity. Lines are posterior means, and uncertainty intervals are 2.5% and 97.5% posterior percentiles (i.e., 95% credible intervals). Notice that overlapping uncertainty intervals do not mean there is no difference between decision curves due to correlation in the sampling distribution of net benefit estimators.](SIM-43-6042-g001.jpg)
In terms of the point estimate of net benefit, the original ADNEX model is superior for most decision under no additional costs to implement or use, ADNEX‐based decisions would be the best strategy to follow. However, there is substantial uncertainty around the clinically motivated threshold of 6%, where it is not immediately clear if the ADNEX superiority is simply due to chance. We may require a low degree of uncertainty before replacing the current Standard of Care (or default strategies) to prevent having to undo the implementation of a worse strategy. However, judging based on the overlap of uncertainty intervals is particularly misleading in DCA due to high correlations across net benefit estimates. Thus, inspecting the DCA alone, under the threshold of 6%, it is unclear whether more evidence confirming ADNEX superiority is required before implementing it into clinical practice.
Additionally, uncertainty intervals from the ADNEX model and the hypothetical SoC test start overlapping at higher decision thresholds as well. The SoC test becomes superior in terms of estimated net benefit at very high thresholds. Here, we consider only decision thresholds below 50% due to the assumption that, in the case of ovarian cancer diagnosis, the cost of false negatives is generally larger than the cost of false positives.
From its decision curve alone, it is clear that the ADNEX model is clinically useful for higher thresholds, but there is considerable uncertainty at thresholds below 10%. Besides visual inspection of the decision curves, how can we quantify our uncertainty about which decision strategies are clinically useful? We can answer this question by interrogating the posterior at the clinically motivated threshold of 6%, there is over 99.9% posterior probability that the ADNEX model is clinically useful—Figure 6A. As the threshold increases, this posterior probability is consistently close to 100%, and the SoC test becomes useful as well. Importantly, we are only able to speak of P(useful) at all and inspect posterior distributions for any decision strategy because we are adopting a Bayesian approach [16]. While Sadatsafavi et al. [16] employ Bayesian bootstrapping to estimate P(useful), no natural counterpart exists under the Frequentist interpretation.

The results above may be counterintuitive due to the overlapping uncertainty intervals in Figure 5. However, notice that P(useful) at 6% for the ADNEX model is determined by the difference between its net benefit and the one from the Treat all strategy—Figure 6B, also notice the inset element within the plot for a zoomed view at the 6% threshold. Their posterior distributions are highly correlated (R = 0.98) due to their shared prevalence parameter, so the estimated net benefit difference is very 0.014 (95% Cr.I. 0.009–0.018). However, to justify the implementation of the ADNEX model at the 6% threshold, we must be confident that this gain in net benefit also overcomes any additional cost of using the model in daily practice, if this isn't already factored in. Since the strategy cost may depend on the local context, we simply report the estimated net benefit difference and its uncertainty. The consumer of the DCA results can then reason if they are confident enough that this estimated difference overcomes their context‐specific costs. If not, more data may still be required prior to implementation.
Although usefulness is important (i.e., being better than Treating all and Treat none), ideally we would like to use the best decision strategy available. Among the strategies under consideration, we can once again interrogate the posterior distribution to quantify our uncertainty about what is the best decision strategy at each decision threshold—Figures 6C and 6D. We do this by comparing each strategy against the best of the remaining strategies (i.e., the “best competitor”). For a given strategy at a given threshold, if the difference in net benefit against its best competitor is positive, then this strategy is better than all remaining strategies.
For instance, at the 6% threshold, the best competitor against the ADNEX model is the Treat all strategy—see Figure 5. In this case, therefore, P(best) and P(useful) over 99.9%. As expected, the probability that the ADNEX model is the best available strategy is virtually 100% for most thresholds. For higher thresholds, however, there is increasingly more overlap between the uncertainty intervals from the ADNEX model and the SoC test. This translates into an increasingly higher P(best) for the SoC test and a progressively lower P(best) for the ADNEX model. At the 50% threshold, there is 84% posterior probability that the SoC test is the best decision strategy available.
From Figures 5 and 6, using the ADNEX model or treating all the patients is likely the best we can do under low decision thresholds. However, the decision‐maker may be interested in higher thresholds. For instance, we may classify patients as high‐risk if their predicted probability of cancer is higher than the prevalence in the present dataset (41%) [2]. The motivation for this is to maintain the proportion of high‐risk patients, according to the model predictions, close to the actual disease prevalence [19]. From Figure 6, we see that P(best) for the ADNEX model at the 41% threshold is far less around 74%. Moreover, the difference in net benefit between the ADNEX model and its best competitor (the SoC test in this case) is almost unnoticeable. Since we have full posterior distributions, we can directly compute a pairwise comparison between the ADNEX model and the SoC test, which is shown in Figure 7. The superiority of the ADNEX model over the SoC test is clear for small thresholds. However, for high thresholds, the difference is smaller. At the 41% threshold, the estimated difference is 0.003 (95% Cr.I. −0.014–0.019). At this threshold, the SoC test is the best competitor against the ADNEX model, so the probability that the ADNEX model is better than the SoC test is again 74%—matches the corresponding P(best) for ADNEX. If no additional strategy costs are considered, maximizing observed net benefit would lead us to favour the implementation of the ADNEX model. However, there is still a 100%–74% = 26% posterior probability that the current Standard of Care is better than the model for this decision threshold. In the face of such large uncertainty, risk aversion may lead us to oppose model implementation unless more data is made available.

Finally, one might wonder if more data are needed to fully characterize the best decision strategy across all thresholds in the target population. It might be that the consequences of uncertainty at, e.g., very high thresholds are too costly in terms of net benefit. With more data, maybe we could confirm if the ADNEX model is indeed superior at the 41% threshold, for instance. To directly quantify the consequences of the current level of uncertainty, Figure 8 shows the validation EVPI for the present case study. As a summary of the posterior distributions, the EVPI shows small, noisy variations. The highest EVPI across all the decision thresholds of interest is around 0.003. Whether this value is relevant depends, again, on the specific context [16]. For instance, if 1,000 clinical decisions regarding the presence of ovarian cancer are made every year for a given population, then the current level of uncertainty is associated with missing at most three net true positive cases of ovarian cancer per year. If 10,000 decisions are made yearly, we could be missing out on up to 30 net true positives every year.

In summary, using the ADNEX model is the best decision strategy across small thresholds assuming no additional strategy costs according to the data. If there exists additional cost of using the model and we accept the decision threshold of 6%, treating all patients or gathering additional data may be preferred. For the 41% threshold, even though the net benefit point estimate for the ADNEX model is slightly higher than for the SoC test, the uncertainty may be too high to justify replacing the Standard of Care currently in place if we are not risk‐neutral. For decision thresholds of 42% and higher, we are increasingly confident that the current Standard of Care is the best decision strategy. Finally, collecting more data from the same population is expected to yield substantial net benefit gain if thousands of clinical decisions are made every year.
Bayesian decision curve analysis was first proposed in the context of meta‐analysis [8]. An immediate advantage was the possibility of calculating P(useful), the probability that a decision strategy is clinically useful, which can also be defined as the probability of being the best decision strategy available [8]. Here, we attempt to clarify terminology by defining the latter as P(best) instead, which differs from P(useful) if we are evaluating more than one predictive model or test at the same time.
The key advantage of Bayesian DCA is the ability to fully interrogate posterior decision curves with an intuitive probabilistic interpretation. We may compute multiple quantities of interest to quantify uncertainty around the answers to fundamental questions such as which decision strategies are useful, what is the best decision strategy, which of two competing strategies is better, and what is the expected consequence of the current level of uncertainty in terms of expected loss in net benefit. These may help the interpretation of DCA if the uncertainty is too large, then more data may be needed to properly choose strategies.
On the other hand, the bootstrap‐based approaches for DCA currently implemented are limited to estimating confidence intervals around net benefit [20]. Moreover, the coverage of these intervals is limited by the distribution of the under a low effective sample size, the bootstrap fails. In external validation studies, overall sample size calculations often result in a few hundred cases required, and studies with around 100 cases are common [21, 22, 23]. Depending on the distribution of the risk predictions, an external validation study with an arguably sufficient overall sample size may still yield low effective sample sizes within the range of decision thresholds of interest. This may pose a problem when the observed risk predictions are concentrated on one side of a given decision threshold (typically very high or very low thresholds), ultimately hiding the possibility of net harm. A low effective sample size may also cause the EVPI to misbehave, varying non‐monotonically with the overall sample size—see figure 5 in Sadatsafavi et al. [16]. While Sadatsafavi et al. [16] proposed a Bayesian bootstrap approach for estimating net benefit, our work proposes specific parametric Bayesian models that do not suffer from the above challenges. Sadatsafavi et al. [16] also mention the need for fully Bayesian methods, which allow the specification of informative priors on specific parameters of net benefit (e.g., prevalence, sensitivity, and specificity). Our approach offers such a possibility.
From a Frequentist perspective, the main limitations of Bayesian DCA are the choice of priors and the interpretation of the uncertainty intervals. We chose weakly informative priors as default in bayesDCA and provided a simulation study to address this valid concern. Under these priors, Bayesian intervals for the binary outcomes case showed near‐perfect empirical coverage, with similar width to the bootstrap intervals. When stronger priors exist (e.g., from past studies), this knowledge may be used to obtain more accurate estimates. For the survival outcomes case, though empirical coverage of the Bayesian approach was reasonable overall, we did observe slight undercoverage in some simulation settings. This might be due to known biases in the estimation of the shape parameter of the Weibull distribution [24]. Future research is warranted to further improve this method, potentially involving reparametrization of the Weibull likelihood to further leverage parameter orthogonality [25]. Informative prior elicitation for Bayesian DCA, in both binary and survival cases, may also be further developed in future work.
Another possible limitation of the Bayesian approach is the computation time due to sampling. In the binary outcomes case, we leveraged parameter orthogonality and posterior independence to allow particularly fast estimation. The bayesDCA implementation is multiple times faster than its bootstrap counterpart and typically takes less than a second, even for large sample sizes. In the survival outcomes case, we do need to use MCMC to estimate the Weibull parameters, and, hence, the Bayesian approach is slower than in the binary outcomes case. From our experience, a typical DCA may take from three to five minutes in this case (using a standard laptop with 12 GB of RAM and 8 cores and running four MCMC chains in parallel, taking 4,000 draws each).
There is debate regarding the value of uncertainty quantification in DCA [26]. From the risk‐neutral point of view, uncertainty informs whether more data is warranted to confidently compare decision strategies, but it wouldn't alter our choice of strategy given the data that is available at the moment of the clinical decision. When forced to choose between two strategies, any net benefit gain (however small) fully justifies choosing one strategy over the other. Yet, having sufficient statistical power to detect a difference in net benefit slightly above zero can be simply unfeasible. If our policy is to forego the implementation of decision strategies with positive (estimated) net benefit gain because of uncertainty, then, on average, we will be losing net benefit—that is, incurring harm. Under risk neutrality, thus, uncertainty motivates future research, but it cannot impact the strategy choice given the available data [5]. One example that arguably satisfies risk neutrality is when choosing between two widely used prognostic models, such as the QRISK2 score and the Framingham equations for predicting the 10‐year risk of cardiovascular disease [27].
When it comes to new decision strategies, however, the decision to replace the well‐accepted Standard of Care may warrant a more conservative approach [5, 7, 28]. Implementing a new clinical prediction model or diagnostic test that turns out, later, to be worse than current practice poses potentially irrecoverable costs. Challenges such as infrastructure development, training, physician adoption, and regulatory approvals cannot be easily reversed. Importantly, patients harmed by new technologies due to premature implementation may face serious unwanted consequences. From the risk‐averse point of view, it may be preferable to forego uncertain net benefit gain to ensure no further loss.
For example, Zhao et al. (2024) externally validated a series of risk prediction models (EsVan models) for bloodstream infection in febrile children with cancer, suggesting that antibiotic treatment could be avoided in patients with a predicted risk below 10% [29]. In routine practice, many participating study sites administered empiric antibiotics to all patients from the target population. If implemented, model‐assisted medical decisions would withhold antibiotics for children predicted to be at low risk. In this setting, model adoption may face risk aversion from prescribing physicians and patient caregivers, who may still prefer to ignore the predictions until they are sufficiently confident of the model's clinical utility. In fact, such aversion to implementation risk was apparent in related work assessing outcomes of patients whose care was guided by the EsVan models [29]. Among 834 patients classified as low‐risk by the model, 176 (21%) still received empirical intravenous antibiotics at presentation due to judgment of the treating physician and other protocol‐defined safety triggers. While model adoption is expected to be partial (e.g., due to clinical factors assumed to be unrelated to model predictions), the establishment of criteria in the study protocol to ignore low‐risk predictions is consistent with a risk‐averse mindset toward implementation.
From a decision theoretical perspective, the possibility of postponing implementation due to uncertainty arises once we consider the implementation of a new model that replaces a given Standard of Care as a decision in itself, which admits its own utility function. The gain in net true positives (or net true negatives) in future model‐based clinical decisions is the reward from an appropriate implementation decision. There is no guarantee that the utility function for the implementation decision is linear in the net benefit scale. This is analogous to having a non‐linear utility for financial reward in a betting decision. More specifically, risk neutrality assumes that the utility loss due to premature implementation of a harmful new model is equal to the utility loss due to postponing the adoption of a useful new model. A risk‐averse decision maker, however, perceives the former as greater than the latter [30, 31]. Per reviewer request, a more formal treatment of these issues is presented in the Supporting information (Section Implementation Decision Rules).
Thus, risk aversion implies that substantial uncertainty around an estimated net benefit gain may prevent the immediate implementation of a newly developed clinical prediction model or diagnostic test. This is not to advocate for questionable inferential procedures such as threshold‐specific p‐values but to point out that caution may be warranted against the implementation of a new strategy that alters well‐established clinical practice when, for instance, its P(best) varies just above 50% for a range of decision thresholds of interest. Judging whether implementation is appropriate given the observed level of uncertainty depends on multiple factors, including implementation costs, clinical context, feasibility of acquiring more data, and, importantly, the decision maker's subjective aversion to the potential risk of replacing the Standard of Care with a worse decision strategy.
The role of risk neutrality goes beyond the field of clinical prediction models and diagnostic tests. When assessing new health technologies, more research may be needed before a final regulatory approval due to uncertainty around an apparent net benefit gain [32]. In this scenario, a relatively risk‐neutral policy grants “approval with research” (AWR), a preliminary approval that allows the technology to be used while further research is conducted. Conversely, a more risk‐averse policy approves “only in research” (OIR), in which case the technology is restricted to research settings only until the required research is finalized [33, 34]. While approving research is often feasible, only research approvals may be necessary if there are significant irrecoverable costs or when approval itself may discourage the necessary research [35]. Thus, while uncertainty initially motivates the need for further research, it may ultimately delay the approval of health technologies until additional, confirmatory data is made available. For the specific case of clinical prediction models and diagnostic tests, it follows that an apparent but uncertain net benefit gain might not be enough to justify implementation, even if future research is planned, as AWR‐like policies are not always appropriate.
Although the above debate warrants further discussion, the present work proposes a Bayesian approach for the estimation of net benefit that is indifferent to whether the end user operates under risk neutrality and how they interpret uncertainty in DCA. The bayesDCA R package allows for an easy and comprehensive characterization of uncertainty in DCA, including its expected consequences through EVPI calculation. While omitting uncertainty may hinder decisions from risk‐averse stakeholders, full reporting of uncertainty allows transparent assessment of DCA results from both the risk‐neutral and the risk‐averse perspectives.
In sum, we propose a method for Bayesian decision curve analysis and provide a freely available implementation in the bayesDCA R package. We hope our contribution will be relevant for studies that involve the external validation of clinical prediction models as well as studies evaluating diagnostic and prognostic tests. Ultimately, the Bayesian DCA workflow may help clinicians and health policymakers adopt informed decisions when choosing and implementing clinical decision strategies.
All analyses were performed with R version 4.2.3 within a fixed Docker image [31, 36]. All code and data to reproduce the results are available on GitHub (https://github.com/giulianonetto/paper‐bayesdca). For the survival case, the bayesDCA R package employs Markov Chain Monte Carlo based on the No U‐Turn Sampler implemented within Stan and accessed via the rstan R package [15, 37]. Processing of posterior samples and data visualization employed the tidyverse meta‐package and the patchwork package [38, 39]. Parallel processing for the simulations was implemented using the furrr package [40]. Pipeline management was implemented using the targets package [41]. The bayesDCA R package is freely available at https://github.com/giulianonetto/bayesdca.
An illustrative model was built using the GUSTO‐I trial dataset [9]. From the full dataset (N=40,830), we held out a randomly selected validation set (N=500, 36 events) and trained a simple logistic regression model on the remaining data based on age, systolic blood pressure, pulse, and Killip class (I–IV). We ran both Bayesian and Frequentist DCA of the fitted model on the validation data to provide an initial illustration under large uncertainty.
We now describe the model used to estimate the net benefit for each decision strategy at each decision threshold. For each patient i=1,2,…,N, suppose we observe the pair (Di,Zi) where Di=1{ithpatient has disease} and Zi=1{ithpatient has positive prediction}. We then
(11) D∼Bernoulli(θ0)Z|D=1∼Bernoulli(θ1)Z|D=0∼Bernoulli(θ2)
where θ0 is the prevalence, θ1 is the sensitivity, and θ2 is 1 minus the specificity. The likelihood function for the observed data 𝒟={(Di,Zi)}i=1N given the parameter vector θ=(θ0θ1θ2)
(12) L(𝒟|θ)=∏i=1Np(di,zi)=∏i=1Np(di)p(zi|di)=∏i=1NBer(di|θ0)×Ber(zi|θ1)di×Ber(zi|θ2)1−di
The parameters θ0, θ1, and θ2 are orthogonal because the likelihood function factorizes
(13) L(𝒟|θ0,θ1,θ2)=∏i=1NBer(di|θ0)×∏i=1NBer(zi|θ1)di×∏i=1NBer(zi|θ2)1−di
Notice that each product term depends on only one parameter. From a Frequentist perspective, this implies that the maximum likelihood estimators for the parameters of interest are asymptotically independent [42]. In the Bayesian approach, under independent priors, the factorization above implies posterior independence. In our case, let θj∼Beta(αj,βj) for j=0,1,2 be our independent priors. The posterior distribution π(θ|𝒟) is then proportional
[Beta(α0,β0)∏i=1NBer(di|θ0)]×[Beta(α1,β1)∏i=1NBer(zi|θ1)di]×[Beta(α2,β2)∏i=1NBer(zi|θ2)1−di]
which is the numerator of the Bayes' Theorem formula and simplifies to
(14) [θ0D+α0−1(1−θ0)ND+β0−1]×[θ1TP+α1−1(1−θ1)FN+β1−1]×[θ2FP+α2−1(1−θ2)TN+β2−1]
where D=∑idi and ND=N−D are the numbers of patients with and without the disease, respectively, TP=∑idizi represents the total number of true positives, FN=∑idi(1−zi) of false negatives, FP=∑i(1−di)zi of false positives, and TN=∑i(1−di)(1−zi) of true negatives. The expression above can already be recognized as the joint density of three independent Beta random variables. Now, let ℬ(a,b)=∫01ta−1(1−t)b−1dt be the Beta function,
(15) p(𝒟)=∫01∫01∫01π(θ0)π(θ1)π(θ2)×L(𝒟|θ)dθ0dθ1dθ2 =∏j=02ℬ(αj,βj)−1 ×∫01[θ0D+α0−1(1−θ0)ND+β0−1]dθ0 ×∫01[θ1TP+α1−1(1−θ1)FN+β1−1]dθ1(from(14)) ×∫01[θ2FP+α2−1(1−θ2)TN+β2−1]dθ2 =∏j=02ℬ(αj,βj)−1×ℬ(D+α0,ND+β0) ×ℬ(TP+α1,FN+β1)×ℬ(FP+α2,TN+β2)
Putting (14) and (15) together and noticing that the normalizing constant that multiplies expression (14) for the posterior distribution is [∏j=02ℬ(αj,βj)−1]/p(𝒟), we have by Bayes' Theorem:
(16) π(θ|𝒟)=Beta(θ0|D+α0,ND+β0)×Beta(θ1|TP+α1,FN+β1)×Beta(θ2|FP+α2,TN+β2)
Notice that the joint posterior distribution is the product of the marginal posterior distributions from each parameter—i.e., posterior independence. Hence, we can simply draw from the marginal posteriors and combine the marginal samples to form a draw from the joint posterior. This result holds for any sample size and not only asymptotically.
Given samples from the joint posterior, we then compute the posterior net benefit as NB|𝒟=θ1·θ0−wt·θ2·(1−θ0), where wt=t/(1−t). The posterior net benefit for the treat all strategy is given by NBall|𝒟=θ0−wt·(1−θ0). Also, within bayesDCA, we sample specificity θ3=1−θ2 directly instead of θ2 to make communication easier with end users—in agreement with Equation (3)—and to make prior specification potentially more intuitive. Finally, notice that θ0 is a common parameter shared by all thresholds and all decision strategies, whereas θ1 and θ2 are threshold‐ and strategy‐specific.
For each patient i=1,2,…,N, suppose we observe (Ti,Ci,Zi) where Ti is the observed survival time for the ith patient, Ci is the censoring indicator, and Zi=1{ithpatient haspositive prediction}, i.e., Zi is 1 if the predicted risk of an event at time horizon τ is above the decision threshold. We then
(17) Z∼Bernoulli(p)(T,C)|Z=1∼Weibull‐Censored(α1,σ1)(T,C)|Z=0∼Weibull‐Censored(α2,σ2)
The likelihood function for the observed data 𝒟={(Ti,Ci,Zi)}i=1N given the parameter vector θ=(p,α1,σ1,α2,σ2)
(18) L(𝒟|θ)=∏i=1NW‐Cens(ti,ci|α1,σ1)zi×W‐Cens(ti,ci|α2,σ2)1−zi×Bern(zi|p)=∏i=1NW‐Cens(ti,ci|α1,σ1)zi×∏i=1NW‐Cens(ti,ci|α2,σ2)1−zi×∏i=1NBern(zi|p)
(19) :=L1(𝒟+|θ1)×L2(𝒟−|θ2)×L3(𝒟0|p)
where we represent the data as 𝒟+={Ti,Ci}i∈[n]:zi=1 (survival data for patients with positive predictions), 𝒟−={Ti,Ci}i∈[n]:zi=0 (survival under negative predictions), 𝒟0={Zi}i=1n (positive prediction indicators), and Weibull parameters as θ1=(α1,σ1) and θ2=(α2,σ2). Hence, under proper independent
(20) π(θ|𝒟)∝π(θ1)L1(𝒟+|θ1)×π(θ2)L2(𝒟−|θ2)×π(p)L1(𝒟0|p)
which then implies posterior independence, and we can sample each parameter independently from their marginal posteriors. In particular, we place a conjugate Beta(a, b) prior on the Bernoulli parameter p so that we have a closed‐form solution for its marginal posterior — default being a=1 and b=1. As there is no such closed form for the Weibull likelihood under right‐censoring, we need to employ MCMC to sample from the marginal posterior of θ1. Within bayesDCA, this is implemented using Stan to sample θ1 as a parameter and p as a generated quantity [15]. The default priors for θ1 are as in (6), though Gamma priors are also allowed. Due to parameter orthogonality and posterior independence, we don't need to sample the nuisance parameter θ2. We can then compute the posterior conditional survival at the time horizon τ as S=exp(τ/σ1)α1 and the posterior net benefit as NB|𝒟=(1−S)·p−wt·S·p. Notice that here both S and p are threshold‐ and strategy‐specific. For the treat‐all strategy, p is a fixed number at 1.
For each simulation setting, we simulate a large population dataset (N=2×106) from which we randomly generate reasonably‐sized samples to perform DCA. We set an expected value of 100 events for the simulated samples and varied the sample size according to each setting's outcome prevalence or incidence. Bayesian DCA was compared with the Frequentist counterpart using available open‐source software, including the packages rmda and dcurves [20, 43].
For the binary outcome simulations, the underlying data‐generating process is as
yi∼Bernoulli(pi)fori=1,2,…,Npi=logit−1(ziTβ),p^i=logit−1(ziTβ^)ziT=1xi1xi2,xi1,xi2∼iidExp(1)
where β is a vector of true coefficients used to generate the data, and β^ is a vector of “estimated” coefficients from the model under investigation—the model we are validating with DCA. The true underlying disease probabilities are represented by pi, and the estimated probabilities are p^i. We choose values of β and β^ to yield settings with a range of values for outcome prevalence, signal‐to‐noise ratio, and model discrimination. The signal‐to‐noise ratio is represented by the maximum possible AUC any model could achieve in a given setting (i.e., the AUC computed using the true disease probabilities pi). The discrimination of the model under investigation is the true AUC computed using p^i in the entire population data.
Here, we choose β and β^ to fix the true prevalence at 1%, 5%, or 30%, and the maximum AUC at either 0.65 or 0.85—these represent low and high signal‐to‐noise ratio settings, respectively, with varying prevalence. We choose β^ so that the example model approximates the maximum AUC in a given setting as well as the true underlying prevalence, but with poor calibration slope (i.e., the coefficients for xi1 and xi2 are exaggerated). This resembles a common scenario of overfitting in risk prediction. The fixed parameters for each simulation setting are provided in Table 1.
For the survival outcome simulations, we employed the simsurv package [44]. Briefly, we sample survival times from a Weibull distribution with specified shape and scale parameters as well as two standard normal covariates with fixed coefficients (under proportional hazards assumption). Then, we sample censoring times from a uniform distribution with support between zero and 24 months, representing a maximum follow‐up time of two years. For each patient in the population, the observed time is the minimum between survival and censoring times. Table 2 shows the fixed parameters for each simulation setting. Notice that simsurv employs a different Weibull parameterization than the one used by Stan and model (18); while the shape parameter (γ in simsurv) is the same, the scale is denoted λ=σ−α. Table 2 uses simsurv notation for easy reproduction.
Because we are not aware of any Frequentist implementation of DCA for survival outcomes that provides uncertainty intervals, we employed custom R code to get bootstrap standard error estimates for net benefit. Using 500 bootstrap samples (same as the default in the rmda package), we passed each bootstrap sample to the dcurves R package. Standard errors were computed as the standard deviation of the resulting bootstrap distribution for the net benefit at each decision threshold.