Skip to contents

Repeated Cross Validated Elastic Net Regression

rcv_glmnet.Rd

Repeatedly runs cross-validated elastic-net models. Coefficient estimates and out-of-sample performance are aggregated for the optimal hyperparameter of each run. For (generalized) linear and multinomial/Cox models it wraps glmnet::cv.glmnet(). For ordinal outcomes it wraps ordinalNet::ordinalNetTune().

Usage

rcv_glmnet(
  data,
  formula,
  engine = c("auto", "glmnet", "ordinalNet"),
  repeats = 50L,
  n_folds = 5L,
  glmnet_s = c("lambda.min", "lambda.1se"),
  ordinal_metric = c("loglik", "misclass", "brier", "devPct"),
  standardize = TRUE,
  symmetry = FALSE,
  weights = NULL,
  offset = NULL,
  subset = NULL,
  na_action = stats::na.omit,
  drop_unused_levels = TRUE,
  xlev = NULL,
  sparse = FALSE,
  ...
)

Arguments

data

(data.frame)
The rectangular dataset containing the response and predictors referenced in formula.

formula

(formula)
Use y ~ . for an additive model using all column in data. For ordinalNet, the response must be an ordered factor (or a factor if you intentionally want nonparallel forms).

engine

(Scalar character: c("auto","glmnet","ordinalNet"))
Selects which backend to use. "auto" chooses "ordinalNet" when the response is ordered, otherwise "glmnet".

repeats

(Scalar integer: [1, Inf))
Number of repeated cross-validated fits to run and aggregate.

n_folds

(Scalar integer: [2, nrow(data)-1])
Number of folds for the outer CV. Mapped to nfolds for glmnet and to nFolds for ordinalNetTune().

glmnet_s

(Scalar character: c("lambda.min","lambda.1se"))
Selection rule used to extract coefficients from cv.glmnet() per repeat.

ordinal_metric

(Scalar character: c("loglik","misclass","brier","devPct"))
Metric to optimize across \(\lambda\) for ordinalNetTune(): maximize "loglik" or "devPct", minimize "misclass" or "brier".

standardize

(Scalar logical)
If TRUE, standardize predictors within folds for both engines.

symmetry

(Scalar logical: FALSE)
For glmnet engine only. If TRUE, categorical predictors are one-hot encoded (k indicators) and the model intercept is dropped in glmnet models. If FALSE, categorical predictors are treatment encoded (k-1 indicators) and the model intercept is kept in glmnet models.

Treatment encoding (k-1 dummies) breaks symmetry among levels of a factor predictor. The omitted level's effect is absorbed into the intercept, which is not penalized, while the other levels' coefficients are shrunk by the elastic-net penalty. This can arbitrarily bias the model results depending on which level is chosen as the reference. Using one-hot encoding ensures all levels are treated equally by the regularization process.

If your goal is estimation of prediction performance, treatment encoding is preferred. If your goal is variable selection, one-hot encoding with dropped intercept is preferred.

As recommended by the ordinalNet authors, categorical variables are always one-hot encoded (k indicators) for this model type.

weights

(Scalar character: NULL)
The column name for the observation weights.

offset

(Scalar character: NULL)
The column name for the offset variable.

subset

(expression or integer index)
Filter rows prior to NA handling and model fit.

na_action

(function)
Used by stats::model.frame() to handle missing data.

drop_unused_levels

(Scalar logical)
If TRUE, drop unused factor levels when constructing the model frame and matrix.

xlev

(named list)
A named list of character vectors giving the full set of levels to be assumed for each factor.

sparse

(Scalar logical)
If TRUE, build a sparse model matrix via Matrix::sparse.model.matrix().

...

Additional engine-specific arguments passed through.
For glmnet, forwarded to cv.glmnet() (e.g., family, alpha, type.measure, foldid, etc.). For ordinalNet, forwarded to ordinalNetTune() (e.g., alpha, family, link, parallelTerms, nonparallelTerms, penaltyFactors, etc.).

Value

A list with components:

  • engine: the backend used ("glmnet" or "ordinalNet").

  • coef_summary: data.frame with columns term, estimate_mean, estimate_sd, and coef_var summarizing \(\hat\beta\) across repeats.

  • coef_matrix: numeric matrix of size p × repeats with per‑repeat coefficients (columns aligned to coef_summary$term).

  • performance: numeric vector of length repeats with the out-of-sample measure per repeat; attribute "measure" specifies the metric.

  • lambda_selected: numeric vector (length repeats) of the selected \(\lambda\) per repeat.

  • call: the matched call for reproducibility.

Details

Model engines

  • If engine = "auto" (default), an ordinal response (an ordered factor) triggers the ordinalNet engine; otherwise the glmnet engine is used, unless overridden.

  • For glmnet, cross-validation uses cv.glmnet() and selects \(s \in \{\lambda_{min}, \lambda_{1se}\}\); the mean CV loss at the selected \(\lambda\) is reported per repeat.

  • For ordinalNet, cross-validation uses ordinalNetTune() (single outer CV across a grid of \(\lambda\)). The best \(\lambda\) is chosen by maximizing average out-of-sample log-likelihood or deviance explained, or minimizing average misclassification or Brier score, consistent with ordinalNetTune() output.

Design matrices from formulas

  • Predictors are created from the formula using model.matrix().

  • For ordinalNet, categorical predictors are one-hot encoded (k indicator columns for k levels) as recommended to avoid reference level dependence of the fit.

  • For glmnet, categorical predictors are treatment encoded (k-1 indicators) if symmetry = FALSE. Categorical predictors are one-hot encoded (k indicators) and the intercept is dropped if symmetry = TRUE.

Aggregation of estimates

Let \(\hat\beta^{(r)} \in \mathbb{R}^p\) denote the coefficient vector (including intercepts) chosen by cross-validation at repeat \(r=1,\dots,R\). The function returns the across-repeat mean \(\bar\beta = \frac{1}{R}\sum_{r=1}^R \hat\beta^{(r)}\), standard deviation \(s_j = \sqrt{\frac{1}{R-1}\sum_{r=1}^R(\hat\beta^{(r)}_j-\bar\beta_j)^2}\), and coefficient of variation \(\mathrm{CV}_j = s_j / |\bar\beta_j|\) for each parameter \(j\).

Performance

  • For glmnet, the out-of-sample metric recorded per repeat is cv.fit$cvm[index(s)], i.e., the mean cross-validated error at the selected \(\lambda\). Its definition depends on type.measure (MSE, deviance, class error, AUC, MAE, C for Cox); see ?cv.glmnet.

  • For ordinalNet, the performance per repeat is the fold-averaged value of the chosen summary (log-likelihood, misclassification, Brier, or deviance explained) at the selected \(\lambda\).

Numerical notes

  • For ordinal models we scale within folds when standardize = TRUE, matching ordinalNet’s CV behavior.

  • For glmnet, cv.glmnet() computes a master \(\lambda\)-sequence and aligns per‑fold fits; selection rules lambda.min and lambda.1se follow the package defaults.

References

Wurm MJ, Rathouz PJ, Hanlon BM (2021). “Regularized Ordinal Regression and the ordinalNet R Package.” Journal of Statistical Software, 99(6), 1–42. doi:10.18637/jss.v099.i06 .

Simon N, Friedman J, Hastie T, Tibshirani R (2011). “Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent.” Journal of Statistical Software, 39(5), 1–13. doi:10.18637/jss.v039.i05 .

Tay JK, Narasimhan B, Hastie T (2023). “Elastic Net Regularization Paths for All Generalized Linear Models.” Journal of Statistical Software, 106(1), 1–31. doi:10.18637/jss.v106.i01 .

Examples

# Gaussian (glmnet)
set.seed(1)
n <- 120
p <- 8
d <- data.frame(
  y = rnorm(n),
  x = matrix(rnorm(n*p), n, p)
)
fit_g <- rcv_glmnet(
  data = d,
  formula = y ~ .,
  engine = "glmnet",
  repeats = 3,
  n_folds = 5,
  glmnet_s = "lambda.min",
  alpha = 0.5,
  family = "gaussian"
)
fit_g$coef_summary
#>          term estimate_mean estimate_sd coef_var
#> 1 (Intercept)     0.1097109           0        0
#> 2         x.1     0.0000000           0       NA
#> 3         x.2     0.0000000           0       NA
#> 4         x.3     0.0000000           0       NA
#> 5         x.4     0.0000000           0       NA
#> 6         x.5     0.0000000           0       NA
#> 7         x.6     0.0000000           0       NA
#> 8         x.7     0.0000000           0       NA
#> 9         x.8     0.0000000           0       NA
mean(fit_g$performance)
#> [1] 0.7957022
attr(fit_g$performance, "measure")
#>                  mse 
#> "Mean-Squared Error" 

# Binomial (glmnet)
set.seed(2)
d2 <- transform(d, y = rbinom(n, 1, plogis(x.1 - 0.5*x.2)))
fit_b <- rcv_glmnet(
  data = d2,
  formula = y ~ .,
  engine = "glmnet",
  repeats = 3,
  n_folds = 5,
  alpha = 0.8,
  family = "binomial",
  type.measure = "class"
)
fit_b$coef_summary
#>          term estimate_mean estimate_sd   coef_var
#> 1 (Intercept)   -0.34170447  0.02001299 0.05856814
#> 2         x.1    0.95750081  0.22185926 0.23170660
#> 3         x.2   -0.26002228  0.15598852 0.59990444
#> 4         x.3   -0.02861764  0.02529668 0.88395407
#> 5         x.4   -0.12100001  0.10488302 0.86680170
#> 6         x.5   -0.11645270  0.10091680 0.86659045
#> 7         x.6    0.00000000  0.00000000         NA
#> 8         x.7    0.18669857  0.16179421 0.86660656
#> 9         x.8    0.22298404  0.19111526 0.85708042

# Ordinal (ordinalNet)
set.seed(3)
n <- 150
p <- 5
X <- matrix(rnorm(n*p), n, p)
eta <- 0.8*X[,1] - 0.6*X[,2]
cut <- c(-0.5, 0.7)
z <- cut(eta + rlogis(n), breaks = c(-Inf, cut, Inf), labels = FALSE)
d3 <- data.frame(y = ordered(z), X)
fit_o <- rcv_glmnet(
  data = d3,
  y ~ .,
  engine = "ordinalNet",
  repeats = 3,
  n_folds = 5,
  alpha = 0.5,
  family = "cumulative",
  link = "logit",
  ordinal_metric = "loglik",
  parallelTerms = TRUE,
  nonparallelTerms = FALSE
)
head(fit_o$coef_summary)
#>            term estimate_mean estimate_sd    coef_var
#> 1 (Intercept):1    -0.4490301 0.014843840 0.033057559
#> 2 (Intercept):2     0.5916671 0.001990385 0.003364028
#> 3            X1    -1.0821751 0.059692778 0.055160002
#> 4            X2     0.5954707 0.039756607 0.066765007
#> 5            X3    -0.1028114 0.019800259 0.192588145
#> 6            X4     0.0000000 0.000000000          NA
mean(fit_o$performance)
#> [1] -26.95848
attr(fit_o$performance, "measure")
#> [1] "mean out-of-fold loglik (higher=better)"