Fits a sparse Likelihood-based Additive Single-Effect Regression (LASER) model using iteratively blockwise coordinate ascent. The model represents the coefficient vector as a sum of sparse "single effects" and produces confidence sets for variable selection based on posterior model probabilities.
Arguments
- X
Numeric matrix of predictors (n × p). If a vector is provided, it will be converted to a single-column matrix.
- y
Response variable:
For GLMs: numeric vector of length n
For Cox models: numeric matrix with 2 columns (time, status) and n rows
- L
Integer specifying the number of components to fit (default: 10). Will be truncated to min(10, ncol(X)) if necessary.
- family
A family object specifying the model type (e.g.,
gaussian(),binomial(),poisson(),Gamma()) orlist(family="cox")for Cox regression.- coverage
Numeric in \([0,1]\) specifying the target coverage probability for confidence sets (default: 0.95).
- cor_threshold
Numeric in \([0,1]\) specifying the minimum absolute correlation required for variables to be grouped in the same confidence set (default: 0.5).
- standardize
Logical indicating whether to center and scale predictors before fitting (default: TRUE).
- decompose
Logical indicating whether to decompose theta in fitting (default: TRUE).
- shrinkage
Logical indicating whether to shrinkage parameters using pvals (default: TRUE).
- tol
Numeric specifying convergence tolerance for the expected log-likelihood between iterations (default: 5e-2).
- lambda
Numeric penalty weight for the truncated-L1 penalty (default: 0.0). If 0, no penalization is applied.
- tau
Numeric truncation parameter for the truncated-L1 penalty (default: 1e-5). Controls the transition from L1 to L0 regularization.
- ties
Character string specifying the method for handling tied event times in Cox regression: "efron" (default) or "breslow".
- max_iter
Integer specifying maximum number of coordinate ascent iterations (default: 100).
- seed
Integer seed for reproducibility (default: NULL).
Value
A list with class "glmsusie" containing:
- call
The matched call
- X
The model matrix
- y
The response vector/matrix
- family
The family object used
- theta
p × L matrix of estimated coefficients for each single effect
- intercept
Estimated intercept (NULL for Cox regression)
- pmp
p × L matrix of posterior model probabilities
- loglik
p × L matrix of log-likelihoods
- bic
p × L matrix of BIC values
- bic_diff
p × L matrix of BIC differences from null model
- evidence
p × L matrix of evidence
- bf
p × L matrix of Bayes factors
- marginal
Vector of marginal inclusion probabilities for each predictor
- kept
Logical vector indicating which effects were retained
- cs
List of confidence sets based on posterior probabilities
- niter
Number of iterations performed
- max_iter
Number of maximum iterations
- elapsed
Elapsed computation time in seconds
Details
The LASER model decomposes the coefficient vector into a sum of L sparse components. At each iteration, the algorithm cyclically updates one component while holding the others fixed. For each component, it fits a univariate model for each predictor, computes model probabilities (via BIC), and updates coefficients as probability-weighted averages. The approach extends traditional GLMs by providing Bayesian-inspired confidence sets for variable selection.
Supported model families include:
Gaussian linear regression
Binomial logistic regression
Poisson regression
Gamma regression
Other GLM family regression
Cox proportional hazards regression
See also
summary.glmsusie for summarizing model results,
coef.glmsusie for extracting coefficients,
plot.glmsusie for plotting results
Examples
if (FALSE) { # \dontrun{
# Gaussian linear regression example
set.seed(42)
n <- 100
p <- 50
X <- matrix(rnorm(n*p), n, p)
colnames(X) <- paste0("X", 1:p)
true_beta <- c(rep(1, 5), rep(0, p-5))
y <- X %*% true_beta + rnorm(n)
# Fit model with 3 components
fit <- glmsusie(X, y, L = 3, family = gaussian())
# Examine results
summary(fit)
plot(fit)
# Extract coefficients
coef(fit)
coef(fit, intercept = TRUE)
# Cox regression example
X <- matrix(rnorm(100*10), 100, 10)
colnames(X) <- paste0("X", 1:10)
times <- rexp(100, rate = exp(0.5 * X[,1] + 0.5 * X[,2]))
status <- rbinom(100, 1, 0.7)
y_cox <- cbind(times, status)
fit_cox <- glmsusie(X, y_cox, L = 2, family = list(family = "cox"))
summary(fit_cox)
} # }