Penalized Regression on RSS (Summary Statistics) Objective

Generalizes lassosum_rss() to support LASSO, MCP, SCAD, L0, L0L1, and L0L2 penalties. Uses coordinate descent on the objective $$\beta^T R \beta - 2 \beta^T z + \mathrm{penalty}(\beta)$$ where $R$ is a (possibly pre-shrunk) LD matrix and $z = \hat\beta / \sqrt{n}$.

Usage

penalized_rss(
  bhat,
  LD,
  n,
  penalty = c("lasso", "MCP", "SCAD", "L0", "L0L1", "L0L2"),
  lambda = exp(seq(log(1e-04), log(0.1), length.out = 20)),
  gamma = NULL,
  alpha = 1,
  lambda0 = 0,
  lambda2 = 0,
  thr = 1e-04,
  maxiter = 10000,
  max_swaps = 100
)

Arguments

bhat: Numeric vector of marginal effect estimates (length p).
LD: A list of LD correlation matrices (one per block), as in lassosum_rss().
n: GWAS sample size (positive scalar).
penalty: Penalty type: "lasso", "MCP", "SCAD", "L0", "L0L1", or "L0L2".
lambda: Numeric vector of regularization parameter values along which to trace a solution path (warm-started, largest-first). For LASSO/MCP/SCAD this is the primary penalty strength; for L0 variants it controls the L1 component.
gamma: Concavity parameter for MCP (default 3) or SCAD (default 3.7). Ignored for LASSO and L0 variants.
alpha: Elastic-net mixing for MCP/SCAD: $l_1 = \lambda \alpha$, $l_2 = \lambda (1-\alpha)$. Default 1 (pure L1, no ridge).
lambda0: L0 penalty weight (number of non-zeros). Required for L0 variants; ignored otherwise. Default 0.
lambda2: L2 penalty weight for L0L2 variant. Default 0.
thr: Convergence threshold. Default 1e-4.
maxiter: Maximum coordinate descent iterations per lambda. Default 10000.
max_swaps: Maximum swap rounds for L0 variants. Default 100. Set to 0 to disable swaps.

Value

A list with components:

beta: p x length(lambda) matrix of coefficient estimates.
lambda: The lambda values used.
conv: Convergence indicators (1 = converged).
loss: Quadratic loss at each lambda.
fbeta: Full penalized objective at each lambda.
nparams: Number of non-zero coefficients at each lambda.
beta_est: Coefficient vector at the lambda minimizing fbeta.

Examples

if (FALSE) { # \dontrun{
set.seed(42)
p <- 10; n <- 100
bhat <- rnorm(p, sd = 0.1)
R <- diag(p)
# MCP
penalized_rss(bhat, list(blk1 = R), n, penalty = "MCP")
# SCAD
penalized_rss(bhat, list(blk1 = R), n, penalty = "SCAD")
# L0
penalized_rss(bhat, list(blk1 = R), n, penalty = "L0", lambda0 = 0.01,
              lambda = c(0))
} # }