Generalized Linear Models with Slab and Shrinkage Estimators

savvy_glm2 extends the classical glm2 function from the glm2 package by embedding a set of shrinkage-based methods within the iteratively reweighted least squares (IRLS) algorithm. These shrinkage methods (implemented via savvy_glm.fit2) are designed to improve convergence and estimation accuracy. The user can specify one or more methods through the model_class argument. When multiple methods are provided (default is c("St", "DSh", "SR", "GSR", "LW", "QIS")), the function can evaluate them in parallel (controlled by the use_parallel argument) and selects the final model based on the lowest AIC.

Usage

savvy_glm2(formula, family = gaussian, data, weights,
                 model_class = c("St", "DSh", "SR", "GSR", "LW", "QIS", "Sh"), subset,
                 na.action, start = NULL, etastart, mustart, offset,
                 control = list(...), model = TRUE,
                 method = "savvy_glm.fit2", x = FALSE, y = TRUE,
                 contrasts = NULL, use_parallel = FALSE, use_robust_start = FALSE, ...)

Arguments

formula: An object of class "formula": a symbolic description of the model to be fitted. As for glm2.
family: A description of the error distribution and link function to be used in the model. As for glm2.
data: An optional data frame, list or environment containing the variables in the model. As for glm2.
weights: An optional vector of weights to be used in the fitting process. As for glm2.
model_class: A character vector specifying the shrinkage method(s) to be used in the underlying fitter savvy_glm.fit2. Allowed values are "St", "DSh", "SR", "GSR", "LW", "QIS", and "Sh". If a single value is provided, only that method is executed. If multiple values are provided, the specified methods are evaluated (in parallel if use_parallel = TRUE) and the one with the lowest AIC is returned. By default, the value is c("St", "DSh", "SR", "GSR", "LW", "QIS"); note that the "Sh" method is considered only if explicitly included.
subset: An optional vector specifying a subset of observations to be used in the fitting process. As for glm2.
na.action: A function which indicates what should happen when the data contain NAs. As for glm2.
start: Starting values for the parameters in the linear predictor. As for glm2. If NULL and use_robust_start = TRUE, robust optimization-based starting values may be calculated automatically for fragile link functions (see Details).
etastart: Starting values for the linear predictor. As for glm2.
mustart: Starting values for the vector of means. As for glm2.
offset: An optional vector specifying a priori known component to be included in the linear predictor during fitting. As for glm2.
control: A list of control parameters to pass to the iterative fitting process. As for glm2.
model: A logical value indicating whether the model frame should be included as a component of the returned value. As for glm2.
method: The method to be used in fitting the model. The default is savvy_glm.fit2. This method uses IRLS with custom optimization methods to ensure better convergence by evaluating different fitting methods and selecting the best one based on AIC. As in glm2, the alternative method "model.frame" returns the model frame and does no fitting.
x: A logical value indicating whether the model matrix used in the fitting process should be returned as a component of the returned value. As for glm2.
y: A logical value indicating whether the response vector used in the fitting process should be returned as a component of the returned value. As for glm2.
contrasts: An optional list. See the contrasts.arg of model.matrix.default. As for glm2.
use_parallel: A logical value specifying whether to evaluate multiple shrinkage methods in parallel. Defaults to FALSE. Setting this to FALSE forces sequential evaluation.
use_robust_start: Logical. If TRUE, uses an optimization-based approach (via the CVXR package) to calculate robust starting values for fragile link functions (e.g., "sqrt"). Defaults to FALSE to save computational time, as standard initialization works well for most typical datasets.
...: Additional arguments to be passed to the low level regression fitting functions. As for glm2.

Value

The value returned by savvy_glm2 has exactly the same structure as that returned by glm2, except for:

method: the name of the fitter function used, which by default is savvy_glm.fit2.
chosen_fit: the name of the chosen fitting method based on AIC.

Details

savvy_glm2 improves upon the standard Generalized Linear Model (GLM) fitting process by incorporating shrinkage estimator functions (St_ost, DSh_ost, SR_ost, GSR_ost, LW_est, QIS_est, and optionally Sh_ost) within the IRLS algorithm. The function begins with initial parameter estimates and iteratively updates the coefficients using the specified shrinkage methods. When multiple methods are specified in model_class, they may be evaluated in parallel (if use_parallel = TRUE), and the final model is selected based on the lowest AIC. In cases where any candidate model returns an NA or non-finite AIC (such as when using quasi-likelihood families), the deviance is used uniformly as the selection criterion.

Robust Starting Values: In situations where the user does not provide start, etastart, or mustart, and the chosen family utilizes specific link functions (such as "log", "sqrt", "inverse", "1/mu^2", or "logit"), standard GLM initialization can sometimes lead to infinite deviance or divergence. If use_robust_start = TRUE, the underlying savvy_glm.fit2 automatically employs a data-driven optimization approach using the CVXR package to calculate stable starting coefficients.

References

Marschner, I. C. (2011). glm2: Fitting Generalized Linear Models with Convergence Problems. The R Journal, 3(2), 12–15. doi:10.32614/RJ-2011-012

Asimit, V., Avramescu, O., Chen, Z., Rivas, D., & Senatore, C. (2026). GLM Solutions via Shrinkage.

Asimit, V., Cidota, M. A., Chen, Z., & Asimit, J. (2025). Slab and Shrinkage Linear Regression Estimation. Retrieved from https://openaccess.city.ac.uk/id/eprint/35005/.

Ledoit, O. and Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis, 88(2):365–411. doi:10.1016/S0047-259X(03)00096-4

Ledoit, O. and Wolf, M. (2022). Quadratic shrinkage for large covariance matrices. Bernoulli, 28(3): 1519-1547. doi:10.3150/20-BEJ1315

Author

Ziwei Chen, Vali Asimit and Claudio Senatore
Maintainer: Ziwei Chen <Ziwei.Chen.3@citystgeorges.ac.uk>

Examples

set.seed(123)
n <- 200
p <- 5
X <- matrix(rnorm(n * p), n, p)
colnames(X) <- paste0("X", 1:p)

linear_predictor <- X %*% c(0.5, -0.5, 0.2, 0, 0)
prob <- 1 / (1 + exp(-linear_predictor))
y <- rbinom(n, 1, prob)

df <- data.frame(y = y, X)
fit1 <- savvy_glm2(y ~ ., data = df,
                   family = binomial(link = "logit"),
                   model_class = c("St", "DSh", "LW"),
                   use_parallel = FALSE)
#> Warning: In glm.fit2_St(x, y, control = control_list, family = family_type) :
#> 1: step size truncated due to increasing deviance
#> 2: inner loop 3; cannot correct step size
print(fit1$chosen_fit)
#> [1] "St"
print(coef(fit1))
#> (Intercept)          X1          X2          X3          X4          X5 
#>  0.11332198  0.31030998 -0.61146901  0.07782194  0.15388397 -0.10330695 

eta <- abs(X %*% c(0.5, -0.2, 0.1, 0, 0)) + 1
mu_quad <- eta^2
y_quad <- rpois(n, lambda = mu_quad)
df_quad <- data.frame(y = y_quad, X)

fit2 <- savvy_glm2(y ~ ., data = df_quad,
                   family = poisson(link = "sqrt"),
                   model_class = c("SR", "St", "QIS"),
                   use_robust_start = TRUE,
                   use_parallel = FALSE)
print(fit2$chosen_fit)
#> [1] "St"
print(coef(fit2))
#> (Intercept)          X1          X2          X3          X4          X5 
#>  1.44440092  0.04253551  0.08344906  0.03890617 -0.05256781  0.04984283