Create a design-aware nested K-fold cross-validation plan

design_nested_cv() creates a lazy design-aware resampling object for nested K-fold cross-validation. It computes the retained row universe from a design_spec() object and creates outer and inner K-fold 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 outer split. The materialized outer fold contains analysis, assessment, and an inner design-aware CV object. Use design_fold() again on outer$inner to materialize inner folds for model selection within the outer analysis rows.

Usage

design_nested_cv(
  spec,
  data,
  n_outer_folds = 5L,
  n_inner_folds = 5L,
  n_iters = 1L,
  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 folds. Objects that inherit from data.frame, including tibbles, are accepted.
n_outer_folds: (numeric scalar: 5L; whole number in [2, unit_n])
Number of outer folds. unit_n is the number of retained rows when group is NULL and the number of retained groups otherwise. Whole-number numeric values such as 5 and 5.0 are accepted and stored as integers.
n_inner_folds: (numeric scalar: 5L; whole number in [2, unit_n - ceiling(unit_n / n_outer_folds)])
Number of inner folds. The upper bound ensures that every outer-analysis set has enough sampling units for inner cross-validation. Whole-number numeric values such as 5 and 5.0 are accepted and stored as integers.
n_iters: (numeric scalar: 1L; whole number in [1, .Machine$integer.max])
Number of independent outer K-fold partitions. Use values greater than one for repeated nested cross-validation.
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 resampling. Use a one-sided formula such as ~ patient_id to keep rows with the same group value together.
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_nested_cv_resamples", "design_resamples") with elements:

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

Each outer split contains analysis, assessment, iter, outer_fold, seed, and inner. Each inner split contains analysis, assessment, iter, outer_fold, inner_fold, and seed. Use design_fold() to convert outer or inner splits into model-ready inputs.

Details

Workflow

A typical nested validation workflow starts with design_spec() to declare the design rules. Call design_nested_cv() after setting the random-number seed to create reproducible outer and inner folds. For each outer split, use the inner folds to choose a candidate model or tuning parameter. Then fit the selected candidate on the full outer analysis input and evaluate it on the outer assessment input.

Keeping split construction separate from fold materialization prevents assessment rows from determining factor levels, dropped categorical terms, column order, or other fitted design-schema details. Outer assessment rows are not used during inner model selection.

Retained rows

The retained row universe is computed once before folds are assigned. 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 outer and inner split stores row identities from the retained universe. When a fold is materialized, row_index values in analysis and assessment refer to row numbers in the original data. Assessment row identities are preserved after split construction so prediction errors can be aligned back to the original data. Assessment predictors and responses are still encoded with the fold-specific analysis schema.

Splits

n_outer_folds controls the number of outer folds in each outer K-fold partition. n_inner_folds controls the number of inner folds within each outer analysis set. n_iters controls the number of independent outer K-fold partitions. The total number of outer splits is n_outer_folds * n_iters. Fold assignment is over retained rows when group is NULL and over retained groups when group is supplied.

Within each iteration, every retained row appears in exactly one outer assessment set. Within each outer split, every outer-analysis row appears in exactly one inner assessment set. Across different iterations, the same row can appear in different outer and inner folds. Call base::set.seed() immediately before design_nested_cv() for reproducible split assignment.

Each outer split carries a seed for its final refit, and each inner split carries a seed for its candidate sweep. design_fold() copies the relevant seed onto the materialized outer and inner folds. See resample_nested_cv() for the seeding recipe and its rationale.

Operational variables

group and strata are optional operational controls for assigning folds. 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.

Response stratification

Use strata = NULL for unstratified folds. Use a one-sided formula such as strata = ~ treatment to balance folds by observed covariate levels. Use strata = ~ 1 to stratify by the raw response.

Response stratification evaluates the response before any response_type conversion. For ordinal-count workflows with rare levels, strata = ~ 1 can reduce the chance that an analysis fold has an all-zero declared-level response column. See resample_nested_cv() for additional stratification and grouping details.

`response_type`

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

The response type is fixed when the nested CV object is created. Every materialized outer and inner fold uses the same response representation for analysis and assessment inputs. Create a separate design nested-CV object when a different response representation is needed.

Examples

#----------------------------------------------------------------------------
# Estimate prediction error after inner-fold model selection
#----------------------------------------------------------------------------
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"
)

candidate_models <- c("age_only", "full")

candidate_columns <- function(inputs, candidate) {
  if (candidate == "age_only") {
    return(c("(Intercept)", "age"))
  }
  colnames(inputs$x)
}

candidate_rmse <- function(fold, candidate) {
  cols <- candidate_columns(fold$analysis, candidate)

  fit <- stats::lm.wfit(
    x = fold$analysis$x[, cols, drop = FALSE],
    y = fold$analysis$y,
    w = fold$analysis$weights,
    offset = fold$analysis$offset
  )

  predicted <- drop(
    fold$assessment$x[, cols, drop = FALSE] %*% fit$coefficients
  )
  predicted <- predicted + fold$assessment$offset
  sqrt(mean((fold$assessment$y - predicted)^2))
}

set.seed(42L)
ncv <- design_nested_cv(
  spec,
  analysis,
  n_outer_folds = 3L,
  n_inner_folds = 2L,
  strata = ~ risk_group
)

nested_results <- lapply(seq_along(ncv$splits), function(i) {
  outer <- design_fold(ncv, i)

  inner_rmse <- vapply(candidate_models, function(candidate) {
    fold_rmse <- vapply(seq_along(outer$inner$splits), function(j) {
      inner <- design_fold(outer$inner, j)
      candidate_rmse(inner, candidate)
    }, numeric(1))
    mean(fold_rmse)
  }, numeric(1))

  selected <- names(which.min(inner_rmse))

  data.frame(
    outer_split = outer$id,
    selected_model = selected,
    rmse = candidate_rmse(outer, selected)
  )
})

nested_results <- do.call(rbind, nested_results)
nested_results
#>   outer_split selected_model      rmse
#> 1      split1           full 0.8564407
#> 2      split2           full 0.8901924
#> 3      split3           full 0.6704498
mean(nested_results$rmse)
#> [1] 0.8056943

#----------------------------------------------------------------------------
# Keep repeated patient observations together
#----------------------------------------------------------------------------
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_ncv <- design_nested_cv(
  repeated_spec,
  repeated,
  n_outer_folds = 3L,
  n_inner_folds = 2L,
  group = ~ patient_id,
  strata = ~ treatment
)

outer_fold <- design_fold(grouped_ncv, 1L)
inner_fold <- design_fold(outer_fold$inner, 1L)
intersect(inner_fold$analysis$case_id, inner_fold$assessment$case_id)
#> character(0)