Bootstrap stability of near-optimal symptom combinations

Repeats the entire exhaustive search on bootstrap resamples of the respondents, and reports how often each symptom combination stays on the plateau and how often it is the outright winner. This answers the question a single search cannot: is the best-performing combination genuinely best, or did it win by a margin that would not survive a different sample of patients?

Usage

bootstrap_stability(
  fit,
  data,
  delta = 1,
  n_boot = 2000,
  seed = 123,
  chunk_size = 2000,
  show_progress = TRUE
)

Arguments

fit: A fitted exhaustive search from score_all_combinations.
data: The data frame that was passed to score_all_combinations, with the 20 PCL-5 item columns symptom_1 through symptom_20.
delta: Numeric. Plateau width in percentage points of balanced accuracy (default 1). A single value.
n_boot: Integer. Number of bootstrap replicates (default 2000).
seed: Integer. Random seed. The function restores the calling session's random number state when it exits.
chunk_size: Integer. Number of combinations processed per block (default 2000). Affects speed and memory only, never the result.
show_progress: Logical. If TRUE (default), report a runtime estimate and display a progress bar.

Value

An object of class ptsdiag_stability: a list of three data frames.

per_set: one row per combination, with combination_id, rank and balanced_accuracy carried over from fit, pi (proportion of replicates in which the combination lay within delta of that replicate's best), p_argmax (proportion of replicates in which it was the best, ties shared equally), n_argmax_outright and n_argmax_tied.
per_replicate: one row per replicate, with replicate, ba_best, plateau_size, plateau_proportion, n_tied_best, n_cases and n_noncases.
summary: one row, with delta, n_boot, n_boot_used, n_candidates, and the mean and 2.5th / 97.5th percentiles of the plateau size and proportion across replicates.

Details

Individuals are resampled with replacement, never combinations. In each replicate the balanced accuracy of every combination is recomputed and the reference point BA_best(b) is taken as the maximum within that replicate. This keeps every comparison on the scale of its own replicate: a resampled balanced accuracy is measured against a maximum drawn from the same resample, not against the original sample's maximum, which resamples would fall short of for reasons having nothing to do with the combination.

Read pi with the selection in mind. For a combination fixed in advance it estimates how often that combination is near-optimal in samples like this one. For the combination that happened to win in this sample it is optimistic, because bootstrap replicates share most of their observations with the sample that singled it out; treat the winner's pi as an upper bound rather than an out-of-sample rate. The plateau size, which does not depend on which combination was selected, is not affected by this.

The computation is arranged as matrix products rather than a loop over combinations. Because the reference diagnosis and every decision rule depend on the data only through the binarized responses (item score >= 2), respondents are first collapsed to unique response patterns. Resampling individuals is then equivalent to drawing pattern counts from a multinomial distribution, and the counts of true and false positives for every combination in every replicate follow from two matrix multiplications. A full 2000-replicate run over all 38,760 six-item combinations takes minutes rather than days.

Before resampling begins, the function recomputes balanced accuracy for a spread of combinations from data and checks the result against fit. This catches a mismatched data argument as well as any disagreement between the matrix path and the package's per-row rule.

Examples

# \donttest{
# A compact 4-symptom search and 50 replicates keep the example fast;
# a published analysis would use the full search and n_boot = 2000
ptsd_data <- rename_ptsd_columns(simulated_ptsd[1:120, ],
                                 id_col = c("patient_id", "age", "sex"))
fit <- score_all_combinations(ptsd_data, n_symptoms = 4, n_required = 3,
                              show_progress = FALSE)

stability <- bootstrap_stability(fit, ptsd_data, n_boot = 50,
                                 show_progress = FALSE)
stability$summary
#>   delta n_boot n_boot_used n_candidates plateau_size_mean plateau_size_lower
#> 1     1     50          50         4845             57.58                  9
#>   plateau_size_upper plateau_proportion_mean plateau_proportion_lower
#> 1            240.775              0.01188442              0.001857585
#>   plateau_proportion_upper
#> 1               0.04969556
head(stability$per_set)
#>   combination_id rank balanced_accuracy   pi   p_argmax n_argmax_outright
#> 1      6_7_11_12    1         1.0000000 1.00 0.32877549                 9
#> 2       4_6_7_12    2         0.9954955 0.88 0.05387592                 0
#> 3       4_6_7_11    3         0.9909910 0.68 0.01133175                 0
#> 4      6_7_11_13    4         0.9909910 0.80 0.01219151                 0
#> 5      6_7_11_17    5         0.9909910 0.72 0.01018441                 0
#> 6      6_7_11_19    6         0.9909910 0.64 0.01520290                 0
#>   n_argmax_tied
#> 1            41
#> 2            15
#> 3             5
#> 4             8
#> 5             5
#> 6             4
# }