grasps: Groupwise Regularized Adaptive Sparse Precision Solution

In grasps: Groupwise Regularized Adaptive Sparse Precision Solution

Groupwise Regularized Adaptive Sparse Precision Solution

Description

Provide a collection of statistical methods that incorporate both element-wise and group-wise penalties to estimate a precision matrix.

Usage

grasps(
  X,
  n = nrow(X),
  membership,
  penalty,
  diag.ind = TRUE,
  diag.grp = TRUE,
  diag.include = FALSE,
  lambda = NULL,
  alpha = NULL,
  gamma = NULL,
  nlambda = 10,
  lambda.min.ratio = 0.01,
  growiter.lambda = 30,
  tol.lambda = 0.001,
  maxiter.lambda = 50,
  rho = 2,
  tau.incr = 2,
  tau.decr = 2,
  nu = 10,
  tol.abs = 1e-04,
  tol.rel = 1e-04,
  maxiter = 10000,
  crit = "BIC",
  kfold = 5,
  ebic.tuning = 0.5
)

Arguments

X	An n \times p data matrix with sample size n and dimension p. A p \times p sample covariance matrix with dimension p.
n	An integer (default = nrow(X)) specifying the sample size. This is only required when the input matrix X is a p \times p sample covariance matrix with dimension p.
membership	An integer vector specifying the group membership. The length of membership must be consistent with the dimension p.
penalty	A character string specifying the penalty for estimating precision matrix. Available options include: "lasso": Least absolute shrinkage and selection operator \insertCitetibshirani1996regression,friedman2008sparsegrasps. "adapt": Adaptive lasso \insertCitezou2006adaptive,fan2009networkgrasps. "atan": Arctangent type penalty \insertCitewang2016variablegrasps. "exp": Exponential type penalty \insertCitewang2018variablegrasps. "lq": Lq penalty \insertCitefrank1993statistical,fu1998penalized,fan2001variablegrasps. "lsp": Log-sum penalty \insertCitecandes2008enhancinggrasps. "mcp": Minimax concave penalty \insertCitezhang2010nearlygrasps. "scad": Smoothly clipped absolute deviation \insertCitefan2001variable,fan2009networkgrasps.
diag.ind	A logical value (default = TRUE) specifying whether to penalize the diagonal elements.
diag.grp	A logical value (default = TRUE) specifying whether to penalize the within-group blocks.
diag.include	A logical value (default = FALSE) specifying whether to include the diagonal entries in the penalty for within-group blocks when diag.grp = TRUE.
lambda	A positive numeric vector specifying the grid for the regularization parameter. The default is NULL, which generates its own lambda sequence based on nlambda and lambda.min.ratio.
alpha	A numeric vector in [0,1] specifying the grid for the mixing parameter balancing the element-wise individual L1 penalty and the block-wise group L2 penalty. An alpha of 1 corresponds to the individual penalty only; an alpha of 0 corresponds to the group penalty only. The default value is a sequence from 0.1 to 0.9 with increments of 0.1. When both alpha and lambda are supplied, all Cartesian combinations of their values are evaluated.
gamma	A numeric value specifying the additional parameter fo the chosen penalty. The default value depends on the penalty: "adapt": 0.5 "atan": 0.005 "exp": 0.01 "lq": 0.5 "lsp": 0.1 "mcp": 3 "scad": 3.7
nlambda	An integer (default = 10) specifying the number of lambda values to generate when lambda = NULL.
lambda.min.ratio	A numeric value > 0 (default = 0.01) specifying the fraction of the maximum lambda value \lambda_{\max} to generate the minimum lambda \lambda_{\min}. If lambda = NULL, a lambda grid of length nlambda is automatically generated on a log scale, ranging from \lambda_{\max} down to \lambda_{\min}.
growiter.lambda	An integer (default = 30) specifying the maximum number of exponential growth steps during the initial search for an admissible upper bound \lambda_{\max}.
tol.lambda	A numeric value > 0 (default = 1e-03) specifying the relative tolerance for the bisection stopping rule on the interval width.
maxiter.lambda	An integer (default = 50) specifying the maximum number of bisection iterations in the line search for \lambda_{\max}.
rho	A numeric value > 0 (default = 2) specifying the ADMM augmented-Lagrangian penalty parameter (often called the ADMM step size). Larger values typically put more weight on enforcing the consensus constraints at each iteration; smaller values yield more conservative updates.
tau.incr	A numeric value > 1 (default = 2) specifying the multiplicative factor used to increase rho when the primal residual dominates the dual residual in ADMM.
tau.decr	A numeric value > 1 (default = 2) specifying the multiplicative factor used to decrease rho when the dual residual dominates the primal residual in ADMM.
nu	A numeric value > 1 (default = 10) controlling how aggressively rho is rescaled in the adaptive-rho scheme (residual balancing).
tol.abs	A numeric value > 0 (default = 1e-04) specifying the absolute tolerance for ADMM stopping (applied to primal/dual residual norms).
tol.rel	A numeric value > 0 (default = 1e-04) specifying the relative tolerance for ADMM stopping (applied to primal/dual residual norms).
maxiter	An integer (default = 1e+04) specifying the maximum number of ADMM iterations.
crit	A character string (default = "BIC") specifying the parameter selection criterion to use. Available options include: "AIC": Akaike information criterion \insertCiteakaike1973informationgrasps. "BIC": Bayesian information criterion \insertCiteschwarz1978estimatinggrasps. "EBIC": extended Bayesian information criterion \insertCitechen2008extended,foygel2010extendedgrasps. "HBIC": high dimensional Bayesian information criterion \insertCitewang2013calibrating,fan2017highgrasps. "CV": k-fold cross validation with negative log-likelihood loss.
kfold	An integer (default = 5) specifying the number of folds used for crit = "CV".
ebic.tuning	A numeric value in [0,1] (default = 0.5) specifying the tuning parameter to calculate for crit = "EBIC".

Value

An object with S3 class "grasps" containing the following components:

hatOmega

The estimated precision matrix.

lambda

The optimal regularization parameter.

alpha

The optimal mixing parameter.

initial

The initial estimate of hatOmega when a non-convex penalty is chosen via penalty.

gamma

The optimal addtional parameter when a non-convex penalty is chosen via penalty.

iterations

The number of ADMM iterations.

lambda.grid

The actual lambda grid used in the program.

alpha.grid

The actual alpha grid used in the program.

lambda.safe

The conservative upper bound used to initialize the line search for \lambda_{\max}, corresponding to each distinct value in alpha.grid, when lambda = NULL.

loss

The optimal k-fold loss when crit = "CV".

CV.loss

Matrix of CV losses, with rows for parameter combinations and columns for CV folds, when crit = "CV".

score

The optimal information criterion score when crit is set to "AIC", "BIC", "EBIC", or "HBIC".

IC.score

The information criterion score for each parameter combination when crit is set to "AIC", "BIC", "EBIC", or "HBIC".

membership

The group membership.

References

\insertAllCited

Examples

library(grasps)

## reproducibility for everything
set.seed(1234)

## block-structured precision matrix based on SBM
sim <- gen_prec_sbm(p = 30, K = 3,
                    within.prob = 0.25, between.prob = 0.05,
                    weight.dists = list("gamma", "unif"),
                    weight.paras = list(c(shape = 20, rate = 10),
                                        c(min = 0, max = 5)),
                    cond.target = 100)
## ground truth visualization
plot(sim)

## n-by-p data matrix
library(MASS)
X <- mvrnorm(n = 20, mu = rep(0, 30), Sigma = sim$Sigma)

## precision matrix: adaptive lasso; BIC
prec <- grasps(X = X, membership = sim$membership, penalty = "adapt", crit = "BIC")

## precision matrix visualization
plot(prec)

## performance
performance(hatOmega = prec$hatOmega, Omega = sim$Omega)

## adjacency matrix: diagonal = 0; raw partial correlations;
##                   no thresholding; weighted network
adj <- prec_to_adj(prec$hatOmega,
                   diag.zero = TRUE, absolute = FALSE,
                   threshold = NULL, weighted = TRUE)

## adjacency matrix visualization
plot(adj)

grasps documentation built on Sept. 13, 2026, 1:07 a.m.