Create a design-aware bootstrap plan

design_bootstrap() creates a lazy design-aware resampling object for bootstrap validation. It computes the retained row universe from a design_spec() object and creates bootstrap split indices on those rows. It does not fit a statistical model and does not build a full-data design matrix.

Use design_fold() to materialize one bootstrap split. Each materialized split fits its design schema on the bootstrap analysis sample only and then encodes the full retained-row assessment set and out-of-bag rows with that analysis schema.

Usage

design_bootstrap(
  spec,
  data,
  n_boot = 1000L,
  strata = NULL,
  group = NULL,
  response_type = c("response", "ordinal_counts")
)

Arguments

spec: (design_spec)
Unfitted design specification returned by design_spec().
data: (data.frame)
Analysis data used to create bootstrap samples. Objects that inherit from data.frame, including tibbles, are accepted.
n_boot: (numeric scalar: 1000L; whole number in [1, .Machine$integer.max])
Number of bootstrap replicates. Whole-number numeric values such as 100 and 100.0 are accepted and stored as integers.
strata: (formula: NULL)
Stratification specification. Use NULL for no stratification. Use ~ 1 to stratify by the raw response. Use a one-sided formula such as ~ treatment or ~ site + treatment to stratify by observed columns.
group: (formula: NULL)
Grouping specification. Use NULL for row-level bootstrap sampling. Use a one-sided formula such as ~ patient_id to draw whole groups with replacement.
response_type: (Scalar character: "response")
Response output mode used by design_fold(). Must be one of "response" or "ordinal_counts".

Value

A list of class c("design_bootstrap_resamples", "design_resamples") with elements:

spec: the original design_spec() object.
data: the original analysis data.
row_index: retained row identities used for bootstrap construction.
response_type: response representation used by materialized splits.
response_levels: declared response levels for ordinal-count output, or NULL.
resamples: the underlying index-only bootstrap_resamples object.
splits: a named list of bootstrap split indices.

Each split contains analysis, assessment, oob, and seed. Use design_fold() to convert a split into model-ready inputs.

Details

Workflow

A typical bootstrap validation workflow starts with design_spec() to declare the design rules. Call design_bootstrap() after setting the random-number seed to create reproducible bootstrap samples. Then iterate over seq_along(bs$splits), materialize each split with design_fold(), fit the model on fold$analysis, and evaluate predictions on fold$assessment, fold$oob, or both.

The full retained-row assessment input is useful for bootstrap optimism correction. The oob input is useful for out-of-bag error estimation.

Keeping split construction separate from fold materialization prevents assessment and out-of-bag rows from determining factor levels, dropped categorical terms, column order, or other fitted design-schema details.

Retained rows

The retained row universe is computed once before bootstrap samples are drawn. The na_action value stored in spec controls whether rows with missing predictor or response values are kept, removed, or rejected. Weights and offsets are checked on retained rows using the same rules as design_fit().

Each split stores row identities from the retained universe. When a split is materialized, row_index values in analysis, assessment, and oob refer to row numbers in the original data. Assessment and OOB row identities are preserved after split construction so prediction errors can be aligned back to the original data. Assessment and OOB predictors and responses are still encoded with the bootstrap-analysis schema.

Splits

n_boot controls the number of bootstrap replicates. Each split contains:

analysis: retained row identities drawn with replacement. Duplicate row identities are expected and their draw order is preserved.
assessment: the full retained-row universe in retained-row order.
oob: retained row identities not present in unique(analysis). The OOB set can be empty for a bootstrap split.
seed: single positive integer for reproducible model fits on this split.

Bootstrap assignment is over retained rows when group is NULL and over retained groups when group is supplied. When group is supplied, all rows from a drawn group enter analysis together. If retained group sizes differ, the length of analysis can vary across bootstrap splits. Each bootstrap split advances the current random-number generator state. Call base::set.seed() immediately before design_bootstrap() for reproducible split assignment.

Each split carries a seed, a single positive integer for reproducible model fits on that split. design_fold() copies this seed onto the materialized fold. See resample_bootstrap() for the seed recipe and its rationale.

Operational variables

group and strata are optional operational controls for drawing bootstrap samples. They are evaluated on the retained rows after the design specification computes the resampling row universe.

If the design formula uses ., variables used only for group or strata must not enter the predictor set through . expansion. Put those variables in case_id, weights, offset, or exclude, or write an explicit predictor formula. Operational variables may also be named explicitly as predictors when that is the intended statistical design.

Stratification

Use strata = NULL for unstratified bootstrap sampling. Use a one-sided formula such as strata = ~ treatment to draw samples within observed covariate levels. Use strata = ~ 1 to stratify by the raw response.

Response stratification evaluates the response before any response_type conversion. Each bootstrap stratum must contain at least two sampling units. See resample_bootstrap() for additional stratification and grouping details.

`response_type`

response_type = "response" returns the ordinary evaluated response in each materialized split. response_type = "ordinal_counts" returns weighted ordinal-count response matrices in analysis$y, assessment$y, and oob$y. Ordinal-count splits require a two-sided bare-symbol factor or ordered-factor response with non-missing retained responses. They also require weights when splits are materialized because the response matrix stores weighted counts.

The response type is fixed when the bootstrap object is created. Every materialized split uses the same response representation for analysis, assessment, and OOB inputs. Create a separate design bootstrap object when a different response representation is needed.

Examples

#----------------------------------------------------------------------------
# Estimate optimism-corrected prediction error
#----------------------------------------------------------------------------
library(bkmodel)

set.seed(1L)
analysis <- data.frame(
  patient_id = sprintf("P%02d", 1:36),
  age = round(rnorm(36, mean = 60, sd = 8)),
  treatment = factor(rep(c("control", "active"), 18)),
  stage = factor(rep(c("I", "II", "III"), each = 12)),
  exposure = runif(36, min = 0.5, max = 2),
  inverse_prob_weight = runif(36, min = 0.8, max = 1.2),
  chart_score = rnorm(36)
)
analysis$risk_group <- interaction(
  analysis$treatment,
  analysis$stage,
  drop = TRUE
)

analysis$outcome <- with(
  analysis,
  5 + 0.25 * age + 1.5 * (treatment == "active") +
    0.8 * (stage == "III") + log(exposure) + rnorm(36)
)

spec <- design_spec(
  outcome ~ .,
  weights = ~ inverse_prob_weight,
  offset = ~ log(exposure),
  case_id = ~ patient_id,
  exclude = ~ chart_score + risk_group,
  encoding = "dummy",
  intercept = TRUE,
  na_action = "na.omit",
  novel_levels = "fail"
)

weighted_rmse <- function(observed, predicted, weights) {
  sqrt(stats::weighted.mean((observed - predicted)^2, weights))
}

fit_lm <- function(inputs) {
  stats::lm.wfit(
    x = inputs$x,
    y = inputs$y,
    w = inputs$weights,
    offset = inputs$offset
  )
}

predict_lm <- function(fit, inputs) {
  predicted <- drop(inputs$x %*% fit$coefficients)
  predicted + inputs$offset
}

full_fit <- design_fit(spec, analysis)
full_inputs <- design_data(full_fit)
full_model <- fit_lm(full_inputs)
apparent_rmse <- weighted_rmse(
  full_inputs$y,
  predict_lm(full_model, full_inputs),
  full_inputs$weights
)

set.seed(42L)
bs <- design_bootstrap(
  spec,
  analysis,
  n_boot = 25L,
  strata = ~ risk_group
)

optimism <- vapply(seq_along(bs$splits), function(i) {
  fold <- design_fold(bs, i)
  model <- fit_lm(fold$analysis)

  bootstrap_error <- weighted_rmse(
    fold$analysis$y,
    predict_lm(model, fold$analysis),
    fold$analysis$weights
  )

  original_error <- weighted_rmse(
    fold$assessment$y,
    predict_lm(model, fold$assessment),
    fold$assessment$weights
  )

  original_error - bootstrap_error
}, numeric(1))

apparent_rmse
#> [1] 0.7509329
apparent_rmse + mean(optimism)
#> [1] 0.8419219

#----------------------------------------------------------------------------
# Estimate out-of-bag prediction error
#----------------------------------------------------------------------------
oob_rmse <- vapply(seq_along(bs$splits), function(i) {
  fold <- design_fold(bs, i)
  model <- fit_lm(fold$analysis)

  weighted_rmse(
    fold$oob$y,
    predict_lm(model, fold$oob),
    fold$oob$weights
  )
}, numeric(1))

mean(oob_rmse)
#> [1] 0.8664824

#----------------------------------------------------------------------------
# Draw repeated patient observations as groups
#----------------------------------------------------------------------------
repeated <- data.frame(
  patient_id = rep(sprintf("P%02d", 1:18), each = 2),
  visit = rep(c("baseline", "follow_up"), times = 18),
  treatment = factor(rep(rep(c("control", "active"), each = 9), each = 2)),
  age = rep(round(rnorm(18, mean = 60, sd = 8)), each = 2),
  outcome = rnorm(36)
)

repeated_spec <- design_spec(
  outcome ~ age + treatment + visit,
  case_id = ~ patient_id,
  encoding = "dummy",
  intercept = TRUE
)

set.seed(7L)
grouped_bs <- design_bootstrap(
  repeated_spec,
  repeated,
  n_boot = 10L,
  group = ~ patient_id,
  strata = ~ treatment
)

grouped_fold <- design_fold(grouped_bs, 1L)
grouped_fold$analysis$case_id[duplicated(grouped_fold$analysis$case_id)]
#>  [1] "P12" "P16" "P11" "P15" "P17" "P17" "P17" "P12" "P12" "P17" "P17" "P17"
#> [13] "P17" "P06" "P08" "P03" "P04" "P08" "P08" "P02" "P06" "P06" "P06" "P06"
#> [25] "P05"