Package 'L0ggm'

L0ggm-package

Description

Implements L0 norm regularization penalties for Gaussian graphical models

Author(s)

Alexander P. Christensen <[email protected]>

See Also

Useful links:

• Report bugs at https://github.com/AlexChristensen/L0ggm/issues

---

Automatic correlations

Description

Automatically computes the proper correlations between continuous and categorical variables. NA values are not treated as categories

Usage

auto_correlate(
  data,
  corr = c("kendall", "pearson", "spearman"),
  ordinal_categories = 7,
  forcePD = TRUE,
  na_data = c("pairwise", "listwise"),
  empty_method = c("none", "zero", "all"),
  empty_value = c("none", "point_five", "one_over"),
  forceReturn = FALSE,
  verbose = FALSE,
  ...
)

Arguments

data	Matrix or data frame. Should consist only of variables to be used in the analysis
corr	Character (length = 1). The standard correlation method to be used. Defaults to "pearson". Using "pearson" will compute polychoric, tetrachoric, polyserial, and biserial correlations for categorical and categorical/continuous correlations by default. To obtain "pearson" correlations regardless, use cor. Other options of "kendall" and "spearman" are provided for completeness and use cor
ordinal_categories	Numeric (length = 1). Up to the number of categories before a variable is considered continuous. Defaults to 7 categories before 8 is considered continuous
forcePD	Boolean (length = 1). Whether positive definite matrix should be enforced. Defaults to TRUE
na_data	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
empty_method	Character (length = 1). Method for empty cell correction in polychoric_matrix. Defaults to "none" Available options: "none" — Adds no value (empty_value = "none") to the empirical joint frequency table between two variables "zero" — Adds empty_value to the cells with zero in the joint frequency table between two variables "all" — Adds empty_value to all in the joint frequency table between two variables
empty_value	Character (length = 1). Value to add to the joint frequency table cells in polychoric_matrix. Defaults to "none". Accepts numeric values between 0 and 1 or specific methods: "none" — Adds no value (0) to the empirical joint frequency table between two variables "point_five" — Adds 0.5 to the cells defined by empty_method "one_over" — Adds 1 / n where n equals the number of cells based on empty_method. For empty_method = "zero", n equals the number of zero cells
forceReturn	Boolean (length = 1). Whether correlation matrix should be forced to return. Defaults to FALSE. Set to TRUE to receive the correlation matrix "as is"
verbose	Boolean (length = 1). Whether messages should be printed. Defaults to FALSE
...	Not actually used but makes it easier for general functionality in the package

Value

A symmetric numeric matrix of dimension $p \times p$, where $p$ is the number of variables in data, with values in $[-1, 1]$ and ones on the diagonal. The correlation method used for each pair of variables depends on their types when corr = "pearson" (the default): polychoric for ordinal-ordinal pairs, polyserial for ordinal-continuous pairs, and Pearson's for continuous-continuous pairs, where ordinal means having at most ordinal_categories unique values. When corr is "kendall" or "spearman", cor is used for all pairs regardless of variable type. Row and column names are inherited from data. When forcePD = TRUE (default) and the resulting matrix is not positive definite, the nearest positive definite matrix is returned via nearPD.

Author(s)

Alexander P. Christensen <[email protected]>

Examples

# Obtain correlations
R <- auto_correlate(basic_smallworld)

---

Toy Small-world Network Data Example

Description

A toy small-world network example (n = 500)

Usage

data(basic_smallworld)

Format

A 500 $\times$ 20 continuous response matrix

Examples

# Generated using the {psychNetwoRk} package:
# basic_smallworld <- psychNetwoRk::simulate_smallworld(
#  nodes = 20, # 20 nodes in the network
#  density = 0.30, # moderate initial lattice connectivity
#  rewire = 0.20, # 20% rewiring probability
#  sample_size = 500 # number of cases = 500
# )$data
data("basic_smallworld")

---

Confusion Matrix Metrics for Edge Comparison and Recovery

Description

Computes many commonly used confusion matrix metrics

Usage

edge_confusion(
  base,
  comparison,
  metric = c("all", "sen", "spec", "ppv", "npv", "fdr", "fom", "ba", "f1", "csi", "mcc"),
  full.names = FALSE
)

Arguments

base	Matrix or data frame. Network that will be treated as the "ground truth" such that a false positive represents an edge that is present in comparison but not in this network
comparison	Matrix or data frame. Network that will be treated as the estimator such that a false positive represents an edge that is present in this network but not in base
metric	Character vector. Defaults to "all" metrics. Available options: "all" — All available metrics (default) "sen" — Sensitivity (True Positive Rate): $\frac{TP}{TP + FN}$ "spec" — Specificity (True Negative Rate): $\frac{TN}{TN + FP}$ "ppv" — Positive Predictive Value (Precision): $\frac{TP}{TP + FP}$ "npv" — Negative Predictive Value: $\frac{TN}{TN + FN}$ "fdr" — False Discovery Rate: $1 - PPV = \frac{FP}{TP + FP}$ "fom" — False Omission Rate: $1 - NPV = \frac{FN}{TN + FN}$ "ba" — Balanced Accuracy: $\frac{Sensitivity + Specificity}{2}$ "f1" — F1 Score (harmonic mean of PPV and Sensitivity): $\frac{2TP}{2TP + FP + FN}$ "csi" — Critical Success Index (Jaccard / Threat Score): $\frac{TP}{TP + FP + FN}$ "mcc" — Matthews Correlation Coefficient: $\frac{TP \times TN - FP \times FN}{\sqrt{(TP+FP)(TP+FN)(TN+FP)(TN+FN)}}$
full.names	Boolean (length = 1). Whether full or abbreviated names should be used. Defaults to FALSE. Set to TRUE for full names

Value

A named numeric vector of the requested confusion matrix metrics. Values are generally in $[0, 1]$, with the exception of "mcc" (Matthews Correlation Coefficient), which ranges from $-1$ to $1$ where 1 indicates perfect agreement, 0 indicates chance-level performance, and -1 indicates perfect disagreement. The vector contains only the elements specified by metric (all ten metrics when metric = "all", which is the default). Names are abbreviated (e.g., "sen") when full.names = FALSE (default), or expanded (e.g., "Sensitivity") when full.names = TRUE. Any metric whose denominator is zero (e.g., sensitivity when no edges exist in base) is returned as NA.

Author(s)

Alexander P. Christensen <[email protected]>

Examples

# Set split
split <- sample(
  1:nrow(basic_smallworld),
  round(nrow(basic_smallworld) / 2)
)

# Estimate networks
split1 <- network_estimation(basic_smallworld[split,])
split2 <- network_estimation(basic_smallworld[-split,])

# Estimate metrics
edge_confusion(split1, split2)

---

L0 Norm Regularized Network Estimation

Description

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)

Usage

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,
  ...
)

Arguments

data	Matrix or data frame. Should consist only of variables to be used in the analysis
n	Numeric (length = 1). Sample size must be provided if data provided is a correlation matrix
corr	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)
na_data	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
penalty	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)$
gamma	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
adaptive	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
nlambda	Numeric (length = 1). Number of lambda values to test. Defaults to 50
lambda_min_ratio	Numeric (length = 1). Ratio of lowest lambda value compared to maximal lambda. Defaults to 0.01
penalize_diagonal	Boolean (length = 1). Should the diagonal be penalized? Defaults to TRUE
ic	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"
ebic_gamma	Numeric (length = 1) Value to set gamma parameter in EBIC (see above). Defaults to 0.50 Only used if ic = "EBIC"
fast	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)
LLA	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)
LLA_threshold	Numeric (length = 1). When performing the Local Linear Approximation, the maximum threshold until convergence is met. Defaults to 1e-03
LLA_iter	Numeric (length = 1). Maximum number of iterations to perform to reach convergence. Defaults to 10000
network_only	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
verbose	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

Value

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:

network	The $p \times p$ partial correlation matrix described above
K	The $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
R	The $p \times p$ regularized covariance matrix returned by GLASSO at the optimal lambda (the w component of the GLASSO solution)
penalty	Character string naming the penalty function used (one of "atan", "exp", "gumbel", "log", "weibull")
lambda	Numeric scalar giving the regularization parameter value selected by the information criterion. Larger values correspond to sparser networks
gamma	Numeric 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
correlation	The $p \times p$ empirical correlation matrix computed from data, used as input to the GLASSO
criterion	Character string naming the information criterion used for model selection (e.g., "bic")
IC	Numeric scalar giving the value of the information criterion at the optimal lambda
MLE	List containing the MLE estimated parameters if adaptive = TRUE and one of "exp", "gumbel", or "weibull" was used

Author(s)

Alexander P. Christensen <[email protected]>

References

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.

Examples

# 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
)

---

Traditional Fit Metrics for Networks

Description

Computes several traditional fit metrics for networks including

• chi-square ($\chi^2$)

• root mean square error of approximation (RMSEA) with confidence intervals

• confirmatory fit index (CFI)

• Tucker-Lewis index (TLI)

• standardized root mean residual (SRMR)

• log-likelihood

• Akaike's information criterion (AIC)

• Bayesian information criterion (BIC)

Usage

network_fit(network, n, S, ci = 0.95)

Arguments

network	Matrix or data frame. A p by p square network matrix
n	Numeric (length = 1). Sample size
S	Matrix or data frame. A p by p square zero-order correlation matrix corresponding with the input network
ci	Numeric (length = 1). Confidence interval for RMSEA. Defaults to 0.95

Value

A named numeric vector of traditional and likelihood-based fit indices. The vector always contains the following elements:

chisq

Chi-square statistic ($\chi^2 = n \cdot F_{ML}$), where $F_{ML}$ is the maximum likelihood discrepancy function between the model-implied and empirical correlation matrices

df

Degrees of freedom: total number of unique off-diagonal correlations minus the number of non-zero edges in network

chisq.p.value

p-value for the chi-square test of exact fit (H0: model-implied covariance equals the population covariance)

RMSEA

Root mean square error of approximation. Values $\leq$ 0.05 indicate close fit; values $\leq$ 0.08 indicate acceptable fit

RMSEA.XX.lower, RMSEA.XX.upper

Lower and upper bounds of the ci-level confidence interval for RMSEA, where XX is the integer percentage (e.g., RMSEA.95.lower and RMSEA.95.upper for a 95% CI)

RMSEA.p.value

p-value for the one-sided test of close fit (H0: RMSEA $\leq$ 0.05)

CFI

Comparative fit index, comparing the target model to an independence (null) baseline. Values $\geq$ 0.95 indicate acceptable fit

TLI

Tucker-Lewis index (non-normed fit index). Values $\geq$ 0.95 indicate acceptable fit; can fall outside $[0, 1]$ for severely misspecified models

SRMR

Standardized root mean residual: the root mean squared difference between the model-implied and observed correlation matrices. Values $\leq$ 0.08 indicate acceptable fit

logLik

Gaussian log-likelihood of the model-implied correlation matrix, assuming zero mean structure (means are not estimated)

AIC

Akaike's information criterion: $-2 \cdot \log L + 2 \cdot E$, where $E$ is the number of non-zero edges in network

BIC

Bayesian information criterion: $-2 \cdot \log L + E \cdot \log(n)$

Author(s)

Alexander P. Christensen <[email protected]>

References

Epskamp, S., Rhemtulla, M., & Borsboom, D. (2017). Generalized network psychometrics: Combining network and latent variable models. Psychometrika, 82(4), 904–927.

Examples

# Obtain correlation matrix
S <- auto_correlate(basic_smallworld)

# Obtain Weibull network
weibull_network <- network_estimation(data = basic_smallworld, LLA = TRUE)

# Obtain fit (expects continuous variables!)
network_fit(network = weibull_network, n = nrow(basic_smallworld), S = S)
# Scaled metrics are not yet available for
# dichotomous or polytomous data!

---

Computes Polychoric Correlations

Description

A fast implementation of polychoric correlations in C. Uses the Beasley-Springer-Moro algorithm (Boro & Springer, 1977; Moro, 1995) to estimate the inverse univariate normal CDF, Genz's approximation (Genz, 2004), a double-precision extension of Drezner and Wesolosky (1990) developed by Genz and Ge, to estimate the bivariate normal CDF, and Newton-Raphson (Fisher scoring, with an exact-Hessian arcsin-reparameterized fallback) for optimization of rho

Usage

polychoric_matrix(
  data,
  na_data = c("pairwise", "listwise"),
  empty_method = c("none", "zero", "all"),
  empty_value = c("none", "point_five", "one_over"),
  ...
)

Arguments

data	Matrix or data frame. A dataset with all ordinal values (rows = cases, columns = variables). Data are required to be between 0 and 11. Proper adjustments should be made prior to analysis (e.g., scales from -3 to 3 in increments of 1 should be shifted by added 4 to all values)
na_data	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
empty_method	Character (length = 1). Method for empty cell correction. Defaults to "none". Available options: "none" — Adds no value (empty_value = "none") to the empirical joint frequency table between two variables "zero" — Adds empty_value to the cells with zero in the joint frequency table between two variables "all" — Adds empty_value to all in the joint frequency table between two variables
empty_value	Character (length = 1). Value to add to the joint frequency table cells. Defaults to "none". Accepts numeric values between 0 and 1 or specific methods: "none" — Adds no value (0) to the empirical joint frequency table between two variables "point_five" — Adds 0.5 to the cells defined by empty_method "one_over" — Adds 1 / n where n equals the number of cells based on empty_method. For empty_method = "zero", n equals the number of zero cells
...	Not used but made available for easier argument passing

Value

A symmetric numeric matrix of dimension $p \times p$, where $p$ is the number of variables (columns) in data. Each off-diagonal entry $[i,j]$ is the polychoric correlation between variables $i$ and $j$ — the estimated Pearson correlation between the latent continuous variables assumed to underlie the observed ordinal categories — with values in $[-1, 1]$. Diagonal entries are $1$. Row and column names are inherited from data. If any variable has zero variance, its corresponding row and column are set to NA and a warning is issued.

Author(s)

Alexander P. Christensen <[email protected]> with assistance from GPT-4

References

Beasley-Moro-Springer algorithm
Beasley, J. D., & Springer, S. G. (1977). Algorithm AS 111: The percentage points of the normal distribution. Journal of the Royal Statistical Society. Series C (Applied Statistics), 26(1), 118-121.

Moro, B. (1995). The full monte. Risk 8 (February), 57-58.

Newton-Raphson (Fisher scoring) optimization
Olsson, U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika, 44(4), 443-460.

Genz bivariate normal approximation
Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing, 14(3), 251-260.

Drezner, Z., & Wesolowsky, G. O. (1990). On the computation of the bivariate normal integral. Journal of Statistical Computation and Simulation, 35(1-2), 101-107.

Examples

# Categorize data (5 equal-frequency categories)
basic_categorized <- apply(
  basic_smallworld, 2, function(x){
    as.numeric(cut(x, breaks = quantile(x, probs = 0:5 / 5), include.lowest = TRUE))
  }
)

# Compute polychoric correlation matrix
correlations <- polychoric_matrix(basic_categorized)

# Randomly assign missing data
basic_categorized[sample(1:length(basic_categorized), 1000)] <- NA

# Compute polychoric correlation matrix with pairwise missing
pairwise_correlations <- polychoric_matrix(
  basic_categorized, na_data = "pairwise"
)

# Compute polychoric correlation matrix with listwise missing
pairwise_correlations <- polychoric_matrix(
  basic_categorized, na_data = "listwise"
)

---