Bayespmtools Package Tutorial

The ISARIC Model

As a case study, we present sample size calculations for the ISARIC 4C model, a risk prediction model for predicting deterioration in patients hospitalized from COVID-19 Gupta et al. (2021). The investigators used electronic hospital records from all 9 health regions of the UK to develop and validate this model. London was left out for external validation while the other 8 regions were used in creating the model.

We assume we plan to validate this model in the London region.

Specifying Evidence (Uncertainty Quantification)

The investigators used leave-one-region out cross validation to estimate out-pf-sample performance metrics of the model. We use this analysis as the source of evidence on model performance for the London region. Based on the exchangeability assumption, the performance of the model in the London region will be a random draw from the predictive distribution of performance metrics observed across the other 8 regions of the UK.

The Bayesian sample size calculation approach requires specifying prior distributions on 1) outcome prevalence (prev), 2) c-statistic (cstat), 3) calibration slope (cal_slp), 4) one of calibration intercept (cal_int), observed-to-expected (O/E) ratio (cal_oe), or mean calibration (cal_mean, difference between average observed and predicted risks - which is the one we use for this example).

These distributions are derived from a meta-analysis of internal-external validation results from the original study. Details of this meta-analysis are provided elsewhere in Sadatsafavi et al. (2026). Here, we directly use the predictive distributions from this meta-analysis.

Evidence can be specified in different ways. Perhaps the most straightforward is a list of R formulas, as below:

evidence <- list(
    prev~beta(mean=0.428, sd=0.030), 
    cstat~beta(mean=0.761, cih=0.773), 
    cal_mean~norm(-0.009, 0.125),  #mean and SD
    cal_slp~norm(0.995, 0.024)     #mean and SD
  )

Note the flexibility in specifying the parameters for distributions. For prevalence, we are parameterizing the Beta distribution by its mean and SD; for c-statistic, we are using mean and the upper bound of 95% confidence intervals (i.e., 97.5th percentile). These values are internally mapped to the parameters of their respective distribution. If parameters are not named, they are taken as the natural parameters for the distribution as defined in R, as in the last two components.

Currently recognized distributions are norm (Normal), logitnorm (Logit-normal), probitnorm (Probit-normal), and beta (Beta). Note the flexibility in specifying the distributions.

Sample size targets and rules

Targets of sample size / precision calculations are generally divided into two groups: statistical metrics of model performance, and the net benefit (NB) as a measure of clinical utility.

For statistical metrics, the user can choose among c-statistic (cstat), mean calibration error (cal_mean), O/E ratio (cal_oe), calibration intercept (cal_int), and calibration slope (cal_slp). For our example, we use the same targets as chosen in Riley et al. (2021): The target 95% confidence interval widths are as follows: C-Statistic: 0.10, O/E ratio: 0.22, Calibration Slope: 0.3.

For these metrics, two types of rules are definable: those that target the expected CI width (ECIW), and those that target a given quantile of these metrics (QCIW). The former can be seen as the Bayesian equivalent of the frequentist sample size calculation in Riley et al. (2021). The latter is unique to our Bayesian framework and can be used to devise ‘assurance-type’ rules: determining the sample size that provides a given level of assurance that the CI width will be below a desired target.

For our case study, we calculate sample sizes corresponding to the expected value of these metrics, as well as sample sizes that will provide 90% assurance that the future CI width will be at most as wide as the specified targets above. The expected confidence intervals widths (eciw) only require a target confidence interval width. Meanwhile, the quantile confidence interval widths (qciw) require the target confidence interval width to be specified along with the desired assurance (quantile). In our example we use a quantile of 0.9 for 90% assurance.

For net benefit (NB) as a decision-theoretic metric, this framework suggests two sample size rules. The first is the Optimality Assurance. An Optimality Assurance of 90% indicates that we would like to be at least 90% confident that the strategy that we declare as the most optimal (highest NB) in the sample is truly the most optimal strategy. The assurance.nb parameter is therefore set at 0.9. The second sample size rule is the voi.nb, and is defined as the ratio of the Expected Value of Sample Information (EVSI) at a given sample size to Expected Value of Perfect Information (EVPI). A ratio of 0.8 means that the planned study is expected to reduce NB loss due to uncertainty by 80%. In our case study, we do not impose a VoI criteria. Instead, we will later calculate VoI metrics at sample size components for other targets. This will help us inspect the efficiency of sample sizes in terms of expected gain in NB.

targets_samp <- list(eciw.cstat=0.1,
                eciw.cal_oe=0.22,
                eciw.cal_slp=0.30,
                qciw.cstat=c(0.9, 0.1),
                qciw.cal_oe=c(0.9, 0.22),
                qciw.cal_slp=c(0.9, 0.30),
                oa.nb=0.90)

Sample size calculator (bpm_valsamp())

This is the main function call that computes the sample size requirements based on our evidence and targets. It returns the required sample size for external validation for each parameter requested in the targets.

samp <- bpm_valsamp(evidence=evidence, #Evidence as a list
                   targets=targets_samp, #Targets (as specified above)
                   n_sim=10000, #Number of Monte Carlo simulations
                   threshold=0.2) #Risk threshold for NB calculations

Let’s take a look at the sample size components.

#Print results
print(samp$results)
#>   eciw.cstat  eciw.cal_oe eciw.cal_slp   qciw.cstat  qciw.cal_oe qciw.cal_slp 
#>          347          440         1022          396          517         1172 
#>        oa.nb 
#>          306

As can be seen above, different components require different sample size estimates. The largest sample size is 1172, dictated by the assurance desired for the calibration slope criterion. The smallest component, at 306 pertains to Optimality Assurance.

Computing precision / VoI for a given sample size (bpm_valprec())

This function returns the parameter values based on the given sample sizes and targets. The utility of this function is in fixed designs where we are interested in evaluating the sufficiency of the available sample size for various targets.

For this example, we will be using the outputs from bpm_valsamp to verify that the 95% CI widths and VoI values for sample size components derived from bpm_valprec are indeed close to the desired ones. As such, we use the same evidence element defined earlier, and create a new target element target_prec, that contains the same parameters target_samp but defined as T for present eciw targets, or 0.9 for present qciw targets.

N_prec <- samp$results  #Sample size components from the main function call.

targets_prec = list(eciw.cstat = T, eciw.cal_oe = T, eciw.cal_slp = T, qciw.cstat = 0.9, qciw.cal_oe = 0.9, qciw.cal_slp = 0.9, oa.nb=T)

Note that for ECIW (as well as Optimality Assurance), we just express our desire to have them calculated. For QCIW values, the quantile value needs to be supplied.

prec <- bpm_valprec(
  N = N_prec, #Sample sizes
  evidence = evidence, #Evidence as a list
  dist_type="logitnorm", #Distribution type for calibrated risks
  method="sample", #Sample based or two-level ("2s") method
  targets = targets_prec, #Targets specified above
  n_sim = 10000, #Number of Monte Carlo simulations
  threshold=0.2) #Risk threshold for NB VoI calculations

The main output of bpm_valprec() is a matrix named ‘results’ that contains all requested precision / VoI metrics at all requested sample sizes. The rows of this matrix are named as the requested sample size (as character), while the columns are named after sample size rules. This allows us to extract the relevant cells from this matrix using the short code below:

res <- prec$results[cbind(samp$results, names(samp$results))]
names(res) <- names(samp$results)

print(res)
#>   eciw.cstat  eciw.cal_oe eciw.cal_slp   qciw.cstat  qciw.cal_oe qciw.cal_slp 
#>   0.10081432   0.21673020   0.30495915   0.09995627   0.21972258   0.30001731 
#>        oa.nb 
#>   0.88560000

These values should be close to the requested precision / assurance / VoI targets. If this is not the case, one should potentially run the analysis with larger values of n_sim, as well as checking if evidence is specified reasonably.