savvy_glm.fit2: Generalized Linear Models Fitting with Slab and Shrinkage...

In savvyGLM: Generalized Linear Models with Slab and Shrinkage Estimators

Generalized Linear Models Fitting with Slab and Shrinkage Estimators

Description

savvy_glm.fit2 is a modified version of glm.fit that enhances the standard iteratively reweighted least squares (IRLS) algorithm by incorporating a set of shrinkage estimator functions. These functions (namely St_ost, DSh_ost, SR_ost, GSR_ost, LW_est, QIS_est, and optionally Sh_ost) provide alternative coefficient updates based on shrinkage principles. The function allows the user to select one or more of these methods via the model_class argument. Additionally, it optionally implements a robust optimization-based method to find starting values for fragile link functions.

Usage

savvy_glm.fit2(x, y, weights = rep(1, nobs),
                        model_class = c("St", "DSh", "SR", "GSR", "LW", "QIS", "Sh"),
                        start = NULL, etastart = NULL, mustart = NULL,
                        offset = rep(0, nobs), family = gaussian(),
                        control = list(), intercept = TRUE,
                        use_parallel = FALSE, use_robust_start = FALSE)

Arguments

x	A numeric matrix of predictors. As for glm.fit.
y	A numeric vector of responses. As for glm.fit.
weights	An optional vector of weights to be used in the fitting process. As for glm.fit.
model_class	A character vector specifying the shrinkage model(s) to be used. Allowed values are "St", "DSh", "SR", "GSR", "LW", "QIS", and "Sh". If a single value is provided, only that method is run. If multiple values are provided, the function runs the specified methods in parallel (if use_parallel = TRUE) and returns the best one based on Akaike Information Criterion (AIC). When the user does not explicitly supply a value for model_class, the default is c("St", "DSh", "SR", "GSR", "LW", "QIS") (i.e. "Sh" is not considered).
start	Starting values for the parameters. As for glm.fit. If NULL and use_robust_start = TRUE, robust starting values may be calculated automatically (see Details).
etastart	Starting values for the linear predictor. As for glm.fit.
mustart	Starting values for the mean. As for glm.fit.
offset	An optional offset to be included in the model. As for glm.fit.
family	A description of the error distribution and link function to be used in the model. As for glm.fit.
control	A list of parameters for controlling the fitting process. As for glm.fit.
intercept	A logical value indicating whether an intercept should be included in the model, default is TRUE. As for glm.fit.
use_parallel	Logical. If TRUE, enables parallel execution of the fitting process. Defaults to TRUE. Set to FALSE for serial execution.
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.

Details

savvy_glm.fit2 extends the classical Generalized Linear Model (GLM) fitting procedure by evaluating a collection of shrinkage-based updates during each iteration of the IRLS algorithm. When multiple shrinkage methods are specified in model_class (the default is c("St", "DSh", "SR", "GSR", "LW", "QIS"), thereby excluding "Sh" unless explicitly provided), if the user explicitly includes "Sh" (for example, model_class = c("St", "DSh", "SR", "GSR", "LW", "QIS", "Sh") or model_class = c("St", "Sh")), then the method Sh_ost is also evaluated.

The function can run these candidate methods in parallel (if use_parallel = TRUE). The final model is selected based on the lowest AIC. In cases where any candidate model returns an NA or non-finite AIC (which may occur when using quasi-likelihood families), the deviance is used uniformly as the model-selection criterion instead.

Robust Starting Values: In situations where the user does not provide start, etastart, or mustart, and the chosen family utilizes a fragile link function (specifically "log", "sqrt", "inverse", "1/mu^2", or bounded links like "logit"), standard GLM initialization can lead to infinite deviance or immediate divergence. If use_robust_start = TRUE, savvy_glm.fit2 employs a novel, data-driven optimization approach using the CVXR package to calculate stable starting coefficients (start). This process minimizes the squared difference between the linear predictor and a safely transformed target mean, enforcing positivity constraints strictly when required by the link function.

In essence, the function starts with initial estimates (either robustly generated or user-supplied) and iteratively updates the coefficients using the chosen shrinkage method(s). The incorporation of these shrinkage estimators aims to improve convergence properties and to yield more stable or parsimonious models under various data scenarios.

Shrinkage Estimators Used in the IRLS Algorithm:

Stein Estimator (St)

Computes a single multiplicative shrinkage factor applied uniformly to all coefficients. The factor is chosen to minimize the overall mean squared error (MSE) of the OLS estimator. It is best suited when the predictors are similarly scaled so that uniform shrinkage is appropriate.

Diagonal Shrinkage (DSh)

Extends the Stein estimator by calculating an individual shrinkage factor for each coefficient. Each factor is derived from the magnitude and variance of the corresponding coefficient, providing enhanced flexibility when predictor scales or contributions differ.

Simple Slab Regression (SR)

Penalizes the projection of the OLS estimator along a fixed, predefined direction. It uses a closed-form solution under a slab (fixed-penalty) constraint, making it suitable when prior information suggests shrinking the coefficients along a known direction.

Generalized Slab Regression (GSR)

Generalizes SR by incorporating multiple, data-adaptive shrinkage directions. It typically employs the leading eigenvectors of the design matrix as the directions, providing adaptive regularization that effectively handles multicollinearity.

Ledoit-Wolf Linear Shrinkage (LW)

Implements linear shrinkage towards a one-parameter matrix where all variances are equal and all covariances are zero. This provides a well-conditioned estimator for large-dimensional covariance matrices, based on the methodology of Ledoit and Wolf (2004b).

Quadratic-Inverse Shrinkage (QIS)

A nonlinear shrinkage estimator derived under Frobenius loss, Inverse Stein's loss, and Minimum Variance loss. It is highly effective for large covariance matrices and is based on the quadratic shrinkage approach of Ledoit and Wolf (2022).

Shrinkage Estimator (Sh)

Computes a non-diagonal shrinkage matrix by solving a Sylvester equation. It transforms the OLS estimator to achieve a lower mean squared error. This estimator is evaluated only when "Sh" is explicitly included in the model_class argument.

Value

The value returned by savvy_glm.fit2 has the same structure as the value returned by glm.fit. It includes the following components:

coefficients	the estimated coefficients.
residuals	the working residuals, that is the residuals in the final iteration of the IWLS fit.
fitted.values	the fitted mean values, obtained by transforming the linear predictors by the inverse of the link function.
R	the upper-triangular factor of the QR decomposition of the weighted model matrix.
rank	the numeric rank of the fitted linear model.
qr	the QR decomposition of the weighted model matrix.
family	the family object used.
linear.predictors	the final linear predictors.
deviance	the deviance of the final model.
aic	the AIC of the final model.
null.deviance	the deviance for the null model.
iter	the number of iterations used.
weights	the final weights used in the fitting process.
prior.weights	the weights initially supplied.
df.residual	the residual degrees of freedom.
df.null	the residual degrees of freedom for the null model.
y	the response vector used.
converged	a logical value indicating whether the IRLS iterations converged.
boundary	a logical value indicating whether the algorithm stopped at a boundary value.
time	the time taken for the fitting process.
chosen_fit	the name of the chosen fitting method based on AIC.

Author(s)

Ziwei Chen, Vali Asimit and Claudio Senatore
Maintainer: Ziwei Chen <[e-mail]>

References

Marschner, I. C. (2011). glm2: Fitting Generalized Linear Models with Convergence Problems. The R Journal, 3(2), 12–15. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3150/20-BEJ1315")}

See Also

glm.fit, glm, glm.fit2

Examples

# Example 1: Standard Poisson regression with log link
set.seed(123)
n <- 100
p <- 5
X <- cbind(1, matrix(rnorm(n * p), n, p))
true_beta <- c(0.5, -0.2, 0.1, 0, 0, 0)

mu <- exp(X %*% true_beta)
y <- rpois(n, lambda = mu)

# Fit using savvy_glm.fit2 with the Stein (St) estimator
fit1 <- savvy_glm.fit2(x = X, y = y,
                       family = poisson(link = "log"),
                       model_class = "St")

print(fit1$coefficients)

# Example 2: Demonstrating Robust Starting Values for a Quadratic Link
eta <- abs(X %*% true_beta) + 1
mu_quad <- eta^2

# Generate Poisson data using the quadratic mean
y_quad <- rpois(n, lambda = mu_quad)
fit2 <- savvy_glm.fit2(x = X, y = y_quad,
                       family = poisson(link = "sqrt"),
                       model_class = c("SR", "St", "LW"),
                       use_parallel = TRUE,
                       use_robust_start = TRUE)

print(fit2$chosen_fit)
print(fit2$coefficients)

savvyGLM documentation built on May 8, 2026, 9:06 a.m.