A general function to estimate Gaussian graphical models using L0 penalty approximations. All penalties are implemented using either single-pass or full Local Linear Approximation (LLA: Fan & Li, 2001; Zou & Li, 2008)
network_estimation(
data,
n = NULL,
corr = c("auto", "pearson", "spearman"),
na_data = c("pairwise", "listwise"),
penalty = c("atan", "exp", "gumbel", "log", "weibull"),
gamma = NULL,
adaptive = TRUE,
nlambda = 50,
lambda_min_ratio = 0.01,
penalize_diagonal = TRUE,
ic = c("AIC", "AICc", "BIC", "BIC0", "EBIC", "MBIC"),
ebic_gamma = 0.5,
fast = TRUE,
LLA = TRUE,
LLA_threshold = 0.001,
LLA_iter = 10000,
network_only = TRUE,
verbose = FALSE,
...
)When network_only = TRUE (default), returns a \(p \times p\)
numeric matrix of partial correlations representing the estimated Gaussian
graphical model. Off-diagonal entry \([i,j]\) is the partial correlation
between variables \(i\) and \(j\) controlling for all other variables,
with values in \([-1, 1]\); a value of zero indicates the absence of an
edge. Diagonal entries are zero. Row and column names are inherited from
data.
When network_only = FALSE, returns a named list with the following
elements:
networkThe \(p \times p\) partial correlation matrix described above
KThe \(p \times p\) estimated inverse covariance (precision) matrix at the optimal lambda. Diagonal entries are the conditional precisions; off-diagonal entries are proportional to partial covariances
RThe \(p \times p\) regularized covariance matrix
returned by GLASSO at the optimal lambda (the w component of the
GLASSO solution)
penaltyCharacter string naming the penalty function used
(one of "atan", "exp", "gumbel", "log",
"weibull")
lambdaNumeric scalar giving the regularization parameter value selected by the information criterion. Larger values correspond to sparser networks
gammaNumeric scalar giving the shape parameter of the
penalty actually used. For adaptive penalties (adaptive = TRUE),
this is the data-derived value; otherwise it is the default or
user-supplied value
correlationThe \(p \times p\) empirical correlation
matrix computed from data, used as input to the GLASSO
criterionCharacter string naming the information criterion
used for model selection (e.g., "bic")
ICNumeric scalar giving the value of the information criterion at the optimal lambda
MLEList containing the MLE estimated parameters if
adaptive = TRUE and one of "exp", "gumbel", or
"weibull" was used
Matrix or data frame. Should consist only of variables to be used in the analysis
Numeric (length = 1).
Sample size must be provided if data provided is a correlation matrix
Character (length = 1).
Method to compute correlations.
Defaults to "auto".
Available options:
"auto" --- Automatically computes appropriate correlations for
the data using Pearson's for continuous, polychoric for ordinal,
tetrachoric for binary, and polyserial/biserial for ordinal/binary with
continuous. To change the number of categories that are considered
ordinal, use ordinal_categories
(see polychoric_matrix for more details)
"pearson" --- Pearson's correlation is computed for all
variables regardless of categories
"spearman" --- Spearman's rank-order correlation is computed
for all variables regardless of categories
For other similarity measures, compute them first and input them
into data with the sample size (n)
Character (length = 1).
How should missing data be handled?
Defaults to "pairwise".
Available options:
"pairwise" --- Computes correlation for all available cases between
two variables
"listwise" --- Computes correlation for all complete cases in the dataset
Character (length = 1).
Defaults to "weibull".
Available options:
"atan" --- Arctangent (Wang & Zhu, 2016)
$$\lambda \cdot (\gamma + \frac{2}{\pi}) \cdot \arctan\left(\frac{|x|}{\gamma}\right)$$
"exp" --- EXP (Wang, Fan, & Zhu, 2018)
$$\lambda \cdot (1 - e^{-\frac{|x|}{\gamma}})$$
"gumbel" --- Gumbel
$$\frac{\lambda}{1 - e^{-1}} \cdot \left(e^{-e^{-\frac{|x|}{\gamma}}} - e^{-1}\right)$$
"log" --- Log (Candes, Wakin, & Boyd, 2008)
$$\frac{\lambda \cdot \log\left(1 + \frac{|x|}{\gamma}\right)}{\log\left(1 + \frac{1}{\gamma}\right)}$$
"weibull" --- Weibull
$$\lambda \cdot \left(1 - e^{-\left(\frac{|x|}{\gamma}\right)^k}\right)$$
Numeric (length = 1). Adjusts the shape of the penalty. Defaults:
"atan" = 0.01
"exp" = 0.01
"gumbel" = 0.01
"log" = 0.10
"weibull" = 0.01
Boolean (length = 1).
Whether data-adaptive (gamma) parameters should be used.
Defaults to TRUE.
Set to FALSE to apply default gamma parameters for
adaptive penalties.
Available options:
"exp"
"gumbel"
"weibull"
When adaptive = TRUE, gamma is set to the 10th percentile of the
distribution fitted to the empirical partial correlations
Numeric (length = 1).
Number of lambda values to test.
Defaults to 50
Numeric (length = 1).
Ratio of lowest lambda value compared to maximal lambda.
Defaults to 0.01
Boolean (length = 1).
Should the diagonal be penalized?
Defaults to TRUE
Character (length = 1). What information criterion should be used for model selection? Available options include:
"AIC" --- Akaike's information criterion: \(-2L + 2E\)
"AICc" --- AIC corrected: \(AIC + \frac{2E^2 + 2E}{n - E - 1}\)
"BIC" --- Bayesian information criterion: \(-2L + E \cdot \log{(n)}\)
"BIC0" --- Bayesian information criterion not (Dicker et al., 2013): \(\log{\large(\frac{D}{n - E}\large)} + \large(\frac{\log{(n)}}{n}\large) \cdot E\)
"EBIC" --- Extended BIC: \(BIC + 4E \cdot \gamma \cdot \log{(p)}\)
"MBIC" --- Modified Bayesian information criterion (Wang et al., 2018): \(\log{\large(\frac{D}{n - E}\large)} + \large(\frac{\log{(n)} \cdot E}{n}\large) \cdot \log{(\log{(p)}})\)
Term definitions:
\(n\) --- sample size
\(p\) --- number of variables
\(E\) --- edges
\(S\) --- empirical correlation matrix
\(K\) --- estimated inverse covariance matrix (network)
\(L = \frac{n}{2} \cdot \log \text{det} K - \sum_{i=1}^p (SK)_{ii}\)
\(D = n \cdot \sum_{i=1}^p (SK)_{ii} - \log \text{det} K\)
Defaults to "BIC"
Numeric (length = 1)
Value to set gamma parameter in EBIC (see above).
Defaults to 0.50
Only used if ic = "EBIC"
Boolean (length = 1).
Whether the glassoFast version should be used
to estimate the GLASSO.
Defaults to TRUE.
The fast results may differ by less than floating point of the original
GLASSO implemented by glasso and should not impact reproducibility much (set to FALSE if concerned)
Boolean (length = 1).
Should Local Linear Approximation be used to find optimal minimum?
Defaults to TRUE to find global minimum based on
convergence (LLA_threshold).
Set to FALSE for a single-pass approximation, which can be
significantly faster (Zou & Li, 2008)
Numeric (length = 1).
When performing the Local Linear Approximation, the maximum threshold
until convergence is met.
Defaults to 1e-03
Numeric (length = 1).
Maximum number of iterations to perform to reach convergence.
Defaults to 10000
Boolean (length = 1).
Whether the network only should be output.
Defaults to TRUE.
Set to FALSE to obtain all output for the
network estimation method
Boolean (length = 1).
Whether messages and (insignificant) warnings should be output.
Defaults to FALSE (silent calls).
Set to TRUE to see all messages and warnings for every function call
Additional arguments to be passed on to auto_correlate
Alexander P. Christensen <alexpaulchristensen@gmail.com>
Log penalty
Candes, E. J., Wakin, M. B., & Boyd, S. P. (2008).
Enhancing sparsity by reweighted l1 minimization.
Journal of Fourier Analysis and Applications, 14(5), 877--905.
Adaptive regularization
Christensen, A. P., & Choi, J. (2026).
Adaptive regularization via extreme value distributions for Gaussian graphical models.
PsyArXiv.
BIC0
Dicker, L., Huang, B., & Lin, X. (2013).
Variable selection and estimation with the seamless-L0 penalty.
Statistica Sinica, 23(2), 929--962.
Local Linear Approximation
Fan, J., & Li, R. (2001).
Variable selection via nonconcave penalized likelihood and its oracle properties.
Journal of the American Statistical Association, 96(456), 1348--1360.
EXP penalty
Wang, Y., Fan, Q., & Zhu, L. (2018).
Variable selection and estimation using a continuous approximation to the L0 penalty.
Annals of the Institute of Statistical Mathematics, 70(1), 191--214.
Atan penalty
Wang, Y., & Zhu, L. (2016).
Variable selection and parameter estimation with the Atan regularization method.
Journal of Probability and Statistics, 2016, 1--12.
Seminal simulation in network psychometrics
Williams, D. R. (2020).
Beyond lasso: A survey of nonconvex regularization in Gaussian graphical models.
PsyArXiv.
One-step Local Linear Approximation
Zou, H., & Li, R. (2008).
One-step sparse estimates in nonconcave penalized likelihood models.
Annals of Statistics, 36(4), 1509--1533.
# Obtain default estimator (adaptive Weibull)
weibull_network <- network_estimation(
data = basic_smallworld, LLA = TRUE
)
# Obtain Atan network
atan_network <- network_estimation(
data = basic_smallworld, penalty = "atan", LLA = TRUE
)
# Obtain static EXP network
exp_network <- network_estimation(
data = basic_smallworld, penalty = "exp",
adaptive = FALSE, LLA = TRUE
)
Run the code above in your browser using DataLab