Getting Started with srlars

Introduction

srlars fits the Fast and Scalable Cellwise-Robust Ensemble (FSCRE) algorithm: a competitive ensemble of n_models sparse sub-models, built on a cellwise-robust foundation (Detect Deviating Cells imputation and wrapping-based robust correlations). “Cellwise” robustness matters because contamination in practice often corrupts individual cells of a data matrix rather than whole observations – a handful of bad measurements scattered across otherwise-good rows – which classical observation-level robust methods are not designed to handle.

This vignette builds one small simulated example and reuses it throughout: we simulate a contaminated dataset, fit srlars() with its default settings, and then see what changes when we adjust its two ensemble-shape controls – how much sub-models are allowed to share variables (max_share) and how small a sub-model is allowed to end up (n_min). It closes with a short, code-free note on cv.srlars(), which chooses max_share automatically instead of by hand as we do below.

library(srlars)
library(mvnfast)
library(cellWise)

Building a test case: contaminated, high-dimensional data

The scenario we simulate: 500 candidate predictors, only a fraction of which are truly related to the response, and a training sample smaller than the number of predictors – a setting where the number of unknowns outnumbers the number of observations. The truly active predictors sit in a few correlated blocks (so that some predictors are legitimately more informative than others), and the true predictor-response relationship is otherwise sparse. We then contaminate a fraction of the cells of the training predictors, leaving the test set clean.

We start with the pieces that describe the true relationship: a block-correlation structure among the active predictors, and a sparse coefficient vector that is nonzero only for those active predictors.

set.seed(100)

n <- 50              # training observations
m <- 2000            # test observations
p <- 500             # candidate predictors
p.active <- 75       # truly active predictors, in blocks below
group.size <- 15     # active predictors per correlated block
n_models <- 10       # ensemble size (K)

# Active predictors sit in correlated blocks; everything else is independent noise.
sigma.mat <- matrix(0, p, p)
sigma.mat[1:p.active, 1:p.active] <- 0.1 # weak correlation across blocks
for (g in 0:(p.active / group.size - 1)) {
  idx <- (g * group.size + 1):(g * group.size + group.size)
  sigma.mat[idx, idx] <- 0.7 # stronger correlation within a block
}
diag(sigma.mat) <- 1

# A sparse, moderate-signal true coefficient vector
true.beta <- c(runif(p.active, 0, 5) * (-1) ^ rbinom(p.active, 1, 0.7),
              rep(0, p - p.active))
sigma <- as.numeric(sqrt(t(true.beta) %*% sigma.mat %*% true.beta)) # signal-to-noise = 1

With the true relationship fixed, generating the actual training and test sets is just sampling predictors and adding noise to the response – the test set stays clean throughout, as a genuine holdout should:

x_train <- mvnfast::rmvn(n, mu = rep(0, p), sigma = sigma.mat)
y_train <- as.numeric(x_train %*% true.beta + rnorm(n, 0, sigma))
colnames(x_train) <- paste0("V", 1:p)

x_test <- mvnfast::rmvn(m, mu = rep(0, p), sigma = sigma.mat)
y_test <- as.numeric(x_test %*% true.beta + rnorm(m, 0, sigma))
colnames(x_test) <- colnames(x_train)

Finally, we contaminate 15% of the cells of the training predictors only. Rather than replacing values with arbitrary noise, each contaminated row’s affected cells are set to a correlation outlier: a combination of values that looks unremarkable one variable at a time, but distorts the multivariate dependence structure DDC and wrapping are specifically designed to catch. The exact linear algebra behind that (contam_correlation() below) isn’t essential reading – what matters is that it plants exactly this kind of cellwise, dependence-breaking contamination into x_train:

contam_correlation <- function(X, prop, sigma_mat, gamma = 3) {
  n <- nrow(X); p <- ncol(X)
  idx <- sample.int(n * p, size = round(n * p * prop))
  rows <- ((idx - 1) %% n) + 1
  cols <- ((idx - 1) %/% n) + 1

  for (i in 1:n) {
    J <- cols[rows == i]
    if (length(J) == 0) next
    if (length(J) == 1) { X[i, J] <- gamma * 3; next }
    SigmaJ <- sigma_mat[J, J, drop = FALSE]
    vmin <- eigen(SigmaJ, symmetric = TRUE)$vectors[, length(J)]
    denom <- mahalanobis(t(vmin), center = rep(0, length(J)), cov = SigmaJ)
    X[i, J] <- gamma * sqrt(length(J)) * (vmin / sqrt(denom))
  }
  X
}

x_train <- contam_correlation(x_train, prop = 0.15, sigma_mat = sigma.mat)

Lastly, one small helper we’ll reuse for every fit below: precision and recall of the selected variables against the known active set, and out-of-sample MSPE (scaled by the noise variance, so that 1 is roughly what a correctly-specified model would achieve):

get_metrics <- function(fit) {
  coefs <- as.numeric(coef(fit))[-1]
  sel <- which(coefs != 0)
  truth <- which(true.beta != 0)

  preds <- as.numeric(predict(fit, x_test))

  c(Precision = length(intersect(sel, truth)) / max(length(sel), 1),
   Recall = length(intersect(sel, truth)) / length(truth),
   MSPE = mean((y_test - preds)^2) / sigma^2,
   `Mean sub-model size` = mean(vapply(fit$active.sets, length, integer(1))))
}

Fitting srlars()

With the data in hand, fitting the default ensemble is a single call. By default, max_share = 1: the n_models sub-models are fully disjoint, so no variable can be selected by more than one of them.

fit_default <- srlars(x_train, y_train,
                      n_models = n_models,
                      tolerance = 1e-4,
                      x_preprocess = "ddc",
                      y_preprocess = "wrap",
                      cor_estimator = "wrap",
                      cv_preprocess = "global",
                      cv_fit = "huber",
                      cv_loss = "huber",
                      cv_folds = 5,
                      compute_coef = TRUE)

Each sub-model gets its own, disjoint set of variables:

knitr::kable(
  data.frame(`Sub-model` = seq_len(n_models),
            `Variables selected` = vapply(fit_default$active.sets, length, integer(1))),
  align = "c"
)

metrics_default <- get_metrics(fit_default)
knitr::kable(t(round(metrics_default, 3)), caption = "srlars() at the default max_share = 1")

Out of 75 truly active predictors, this run recovers a recall of 0.64 at a precision of 0.68 – exactly how well any particular run does will vary with the random contamination draw, but the shape of the result (some but not all of the true signal recovered, most of what’s selected genuinely active) is typical of this kind of high-dimensional, heavily contaminated setting.

coef() and predict() both average across the ensemble’s sub-models, so they work exactly like the corresponding methods for an ordinary fitted regression:

knitr::kable(t(round(coef(fit_default)[1:6], 3)),
            col.names = c("Intercept", paste0("V", 1:5)))

knitr::kable(t(round(predict(fit_default, x_test[1:5, ]), 2)),
            col.names = paste("Test row", 1:5))

Controlling how much sub-models can share: max_share

max_share (from 1 to n_models) relaxes the disjointness above: a variable may now be selected by up to max_share sub-models instead of just one. At the other extreme, max_share = n_models removes the restriction entirely – sub-models are then completely free to converge on the same variables. Refitting on the exact same data, changing only max_share:

fit_shared <- srlars(x_train, y_train,
                     n_models = n_models,
                     max_share = n_models,
                     tolerance = 1e-4,
                     x_preprocess = "ddc",
                     y_preprocess = "wrap",
                     cor_estimator = "wrap",
                     cv_preprocess = "global",
                     cv_fit = "huber",
                     cv_loss = "huber",
                     cv_folds = 5,
                     compute_coef = TRUE)

metrics_shared <- get_metrics(fit_shared)
knitr::kable(
  rbind(`max_share = 1 (default)` = round(metrics_default, 3),
       `max_share = n_models`    = round(metrics_shared, 3))
)

Here, allowing unrestricted sharing moves recall from 0.64 to 0.15: with fewer independent sub-models exploring different variables, the ensemble as a whole tends to cover less of a broad, many-active-variable truth like this one. That trade is not universal, though – letting sub-models refit the same, well-supported variables can help on datasets that don’t need broad exploration in the first place. cv.srlars() (see the closing section) picks max_share for you by cross-validation, rather than requiring this kind of manual comparison.

Guaranteeing a minimum sub-model size: n_min

The selection loop’s stopping rule is evaluated once per round, across the entire ensemble at once: it halts as soon as no sub-model’s next candidate variable shows a sufficient cross-validated improvement. On some datasets that can leave sub-models quite small. n_min sets a floor under that: sub-models below it keep receiving their best available variable even when it doesn’t clear the usual improvement bar (though it can never violate the max_share restrictions above – only the improvement requirement is relaxed).

fit_floor <- srlars(x_train, y_train,
                    n_models = n_models,
                    n_min = 10,
                    tolerance = 1e-4,
                    x_preprocess = "ddc",
                    y_preprocess = "wrap",
                    cor_estimator = "wrap",
                    cv_preprocess = "global",
                    cv_fit = "huber",
                    cv_loss = "huber",
                    cv_folds = 5,
                    compute_coef = TRUE)

metrics_floor <- get_metrics(fit_floor)
knitr::kable(
  rbind(`n_min = NULL (default)` = round(metrics_default, 3),
       `n_min = 10`               = round(metrics_floor, 3))
)

Forcing a floor of 10 variables per sub-model raises the mean sub-model size from 7.1 to 10.8, which typically trades some precision (a few of the forced-in variables are not genuinely active) for higher recall – useful when the default stopping rule is cutting sub-models off before they’ve captured much real signal, but not something to reach for by default.

Choosing max_share automatically: cv.srlars()

Rather than comparing max_share values by hand as above, cv.srlars() chooses it by an outer cross-validation loop scored on held-out ensemble prediction error, and returns the refit at the cross-validated optimum. The expensive cellwise-robust preprocessing stage is computed once per outer fold and reused across every candidate value, so this is efficient relative to a naive grid search. The object it returns is classed so that coef() and predict() work on it exactly as they do above – no new methods to learn. See ?cv.srlars for details.