Arguments in `EBICglasso()` not passing to `glasso()`

[mvanaman]

The documentation for EBICglasso() indicates that you can pass arguments to glasso(). However, it doesn't seem to be working: I tried to pass the zero = argument, but it's complaining the function is not found. This also happened for start =.

Here is an example:

library(tidyverse)
#> Warning: package 'tidyverse' was built under R version 4.1.2
#> Warning: package 'tidyr' was built under R version 4.1.2
#> Warning: package 'readr' was built under R version 4.1.2
#> Warning: package 'dplyr' was built under R version 4.1.2
library(glasso)
library(qgraph)

example.data <- iris[, -5]

cors <- cor_auto(example.data, detectOrdinal = FALSE)
SDs <- apply(select(example.data, colnames(cors)), 2, FUN = "sd") %>% 
  data.frame("sd" = .) %>% 
  rownames_to_column("var")

cor2cov_1 <- function(R,S){
    diag(S) %*% R %*% diag(S)
}
covs <- cor2cov_1(R = cors, S = SDs$sd)
colnames(covs) <- names(as.data.frame(cors))

wm.example <- EBICglasso(S = covs, n = nrow(example.data), zero = c(3, 1))
#> Error in glassopath(S, lambda, trace = 0, penalize.diagonal = penalize.diagonal, : unused argument (zero = c(3, 1))

^{Created on 2022-08-24 by the reprex package (v2.0.1)}

Am I correct that the expected behavior is for zero = to be passed from EBICglasso() to glasso()?

[mvanaman]

This happens because zero = is accepted by glasso(), but not glassopath(). But the ... in the code for EBICglasso() runs to both of them. Also, if penalizeMatrix = is empty, then glassopath gets run. So my very-hacky solution was to copy-and-paste the EBICglasso() function (below) and manually add in the argument zero = zero to every call to glasso(), then assign zero in the environment so that it gets picked up by glasso(). The caveat is that you have to give penalizeMatrix = a matrix, so I created one with all TRUE's (which I assume is the default) and added that argument to the EBICglasso() function call.

So this code will work, after running the function at the bottom. Definitely not an ideal long-term solution, but works for my case:

# load packages
packages <- c("tidyverse", "glasso", "qgraph", "Matrix")
suppressMessages(lapply(packages, require, character.only = TRUE))

# get correlation matrix
example.data <- iris[, -5]
cors <- cor_auto(example.data, detectOrdinal = FALSE)

# fit network
zero <- c(4, 3)
pen_mat <- matrix(rep(TRUE, ncol(covs) * nrow(covs)), ncol = ncol(covs), nrow = nrow(covs))

wm.example <- EBICglasso(S = cors, n = nrow(example.data), penalizeMatrix = pen_mat)
qgraph(wm.example, edge.labels = TRUE, curveAll = TRUE)

And the function:

# need this:
## Computes the EBIC:
EBIC <- function(S,K,n,gamma = 0.5,E,countDiagonal=FALSE)
{
    #   browser()
    L <- logGaus(S, K, n)
    if (missing(E)){
        E <- sum(K[lower.tri(K,diag=countDiagonal)] != 0)
        # E <- sum(abs(K[lower.tri(K,diag=countDiagonal)]) > sqrt(.Machine$double.eps))
    }
    p <- nrow(K)

    # return EBIC:
    -2 * L + E * log(n) + 4 * E * gamma * log(p)
}

# Computes optimal glasso network based on EBIC:
EBICglassoCore <- function(
    S, # Sample covariance matrix
    n, # Sample size
    gamma = 0.5,
    penalize.diagonal = FALSE, # Penalize diagonal?
    nlambda = 100,
    lambda.min.ratio = 0.01,
    returnAllResults = FALSE, # If true, returns a list
    checkPD = TRUE, # Checks if matrix is positive definite and stops if not
    penalizeMatrix, # Optional logical matrix to indicate which elements are penalized
    countDiagonal = FALSE, # Set to TRUE to get old qgraph behavior: conting diagonal elements as parameters in EBIC computation. This is not correct, but is included to replicate older analyses
    refit = TRUE, # If TRUE, network structure is taken and non-penalized version is computed.
    ebicMethod = c("new","old"),
    regularized = TRUE,
    threshold = FALSE,
    verbose = TRUE,
    criterion = "ebic",
    ... # glasso arguments
) {
    ebicMethod <- match.arg(ebicMethod)

    if (checkPD){
        if (any(eigen(S)$values < 0)) stop("'S' is not positive definite")
    }

    # Standardize cov matrix:
    S <- cov2cor(S)

    # Compute lambda sequence (code taken from huge package):
    lambda.max = max(max(S - diag(nrow(S))), -min(S - diag(nrow(S))))
    lambda.min = lambda.min.ratio*lambda.max
    lambda = exp(seq(log(lambda.min), log(lambda.max), length = nlambda))

    # Run glasso path:
    if (missing(penalizeMatrix)){
        glas_path <- glassopath(S, lambda, trace = 0, penalize.diagonal=penalize.diagonal, ...)
    }else{
        glas_path <- list(
            w = array(0, c(ncol(S), ncol(S), length(lambda))),
            wi = array(0, c(ncol(S), ncol(S), length(lambda))),
            rholist = lambda
        )

        for (i in 1:nlambda){
            res <- glasso::glasso(S, penalizeMatrix * lambda[i], trace = 0, penalize.diagonal=penalize.diagonal, zero = zero, ...)
            glas_path$w[,,i] <- res$w
            glas_path$wi[,,i] <- res$wi
        }
    }

    # Threshold:
    if (threshold){
        for (i in 1:nlambda){
            # Degree:
            p <- ncol(glas_path$wi[,,i])
            # D <- max(centrality(ifelse( glas_path$wi[,,i] != 0,1, 0))$OutDegree)
            thresh <- (log(p*(p-1)/2)) / sqrt(n)
            glas_path$wi[,,i] <- ifelse(abs(glas_path$wi[,,i]) < thresh,0,glas_path$wi[,,i])

        }
    }

    # Compute EBICs:
    if (ebicMethod == "old"){
        if (criterion != "ebic") stop("criterion must be 'ebic' when ebicMethod = 'old'")
        EBICs <- sapply(seq_along(lambda),function(i){
            if (!regularized){
                invSigma <- ggmFit(wi2net(glas_path$wi[,,i]), S, sampleSize = n, ebicTuning = gamma, refit = TRUE,verbose = FALSE)$invSigma
            } else {
                invSigma <- glas_path$wi[,,i]
            }
            EBIC(S, invSigma, n, gamma, countDiagonal=countDiagonal)
        })
    } else {
        EBICs <- sapply(seq_along(lambda),function(i){
            fit <- ggmFit(wi2net(glas_path$wi[,,i]), S, n, ebicTuning = gamma,refit = !regularized, verbose = FALSE)
            # print(fit$fitMeasures$ebic)
            # browser()
            fit$fitMeasures[[criterion]]
        })
    }

    # Smallest EBIC:
    opt <- which.min(EBICs)

    # Check if rho is smallest:
    if (opt == 1){
        message("Note: Network with lowest lambda selected as best network: assumption of sparsity might be violated.")
    }

    # Return network:
    net <- as.matrix(forceSymmetric(wi2net(glas_path$wi[,,opt])))
    colnames(net) <- rownames(net) <- colnames(S)

    # Check empty network:
    if (all(net == 0)){
        message("An empty network was selected to be the best fitting network. Possibly set 'lambda.min.ratio' higher to search more sparse networks. You can also change the 'gamma' parameter to improve sensitivity (at the cost of specificity).")
    }

    # Refit network:
    # Refit:
    if (refit){
        if (verbose) message("Refitting network without LASSO regularization")
        if (!all(net[upper.tri(net)]!=0)){
            glassoRes <- suppressWarnings(glasso::glasso(S, 0, zero = which(net == 0 & upper.tri(net), arr.ind=TRUE), trace = 0, penalize.diagonal=penalize.diagonal, zero = zero, ...))
        } else {
            glassoRes <- suppressWarnings(glasso::glasso(S, 0, trace = 0, penalize.diagonal=penalize.diagonal, zero = zero, ...))
        }
        net <- as.matrix(forceSymmetric(wi2net(glassoRes$wi)))
        colnames(net) <- rownames(net) <- colnames(S)
        optwi <- glassoRes$wi
    } else {
        optwi <- glas_path$wi[,,opt]
    }

    # If regularized and low lambda was selected, give warning:
    if (regularized && lambda[opt] < 0.1 * lambda.max && !isTRUE(threshold)){
        warning("A dense regularized network was selected (lambda < 0.1 * lambda.max). Recent work indicates a possible drop in specificity. Interpret the presence of the smallest edges with care. Setting threshold = TRUE will enforce higher specificity, at the cost of sensitivity.")
    }

    # Return
    if (returnAllResults){
        return(list(
            results = glas_path,
            ebic = EBICs,
            # loglik = lik,
            optnet = net,
            lambda = lambda,
            optwi = optwi
        ))
    } else return(net)
}

# Old function:

# Computes optimal glasso network based on EBIC:
EBICglasso <- function(
    S, # Sample covariance matrix
    n, # Sample size
    gamma = 0.5,
    penalize.diagonal = FALSE, # Penalize diagonal?
    nlambda = 100,
    lambda.min.ratio = 0.01,
    returnAllResults = FALSE, # If true, returns a list
    checkPD = TRUE, # Checks if matrix is positive definite and stops if not
    penalizeMatrix, # Optional logical matrix to indicate which elements are penalized
    countDiagonal = FALSE, # Set to TRUE to get old qgraph behavior: conting diagonal elements as parameters in EBIC computation. This is not correct, but is included to replicate older analyses
    refit = FALSE, # If TRUE, network structure is taken and non-penalized version is computed.
    threshold = FALSE,
    verbose = TRUE,
    # ebicMethod = c("new","old"),
    # ebicRefit = FALSE,
    ... # glasso arguments
) {
    EBICglassoCore(S=S, # Sample covariance matrix
                   n=n, # Sample size
                   gamma = gamma,
                   penalize.diagonal = penalize.diagonal, # Penalize diagonal?
                   nlambda = nlambda,
                   lambda.min.ratio = lambda.min.ratio,
                   returnAllResults = returnAllResults, # If true, returns a list
                   checkPD = checkPD, # Checks if matrix is positive definite and stops if not
                   penalizeMatrix = penalizeMatrix, # Optional logical matrix to indicate which elements are penalized
                   countDiagonal = countDiagonal, # Set to TRUE to get old qgraph behavior: conting diagonal elements as parameters in EBIC computation. This is not correct, but is included to replicate older analyses
                   refit = refit, # If TRUE, network structure is taken and non-penalized version is computed.
                   ebicMethod = "old",
                   regularized = TRUE,
                   threshold=threshold,
                   verbose=verbose,
                   ...)
}

^{Created on 2022-08-24 by the reprex package (v2.0.1)}

[SachaEpskamp]

Hi! Why would you want to constrain certain edges to zero exactly? If you have a large dataset you can also fit confirmatory models in the psychonetrics pacgage (see psychonetrics.org) in which you can set certain edges to zero and freely estimate other edges.

[mvanaman]

Hi Sacha, thanks for getting back to me! I previously checked out the psychonetrics package... don't know how I missed that!

Is there a way to test for differences in e.g., node centrality between two networks in psychonetrics?

My conundrum is I'd like to compare the constrained and unconstrained versions of the network, but re-estimating the network with different packages (for example, once in psychonetrics, then again in NetworkComparisonTest or bootnet) makes me nervous given how vastly different the networks can turn out to be.

P.S. The reason I want to set some edges to zero is so that I can compare a model suggested by the literature (some edges constrained to zero) with a purely exploratory model (all edges freely estimated) in node centrality, shortest paths, etc.

Thank you so much for your time!

[SachaEpskamp]

Hi! Sorry for the late reply. I think this can be done easier if you just assign a dummy matrix for the penalizeMatrix argument:

# load packages
packages <- c("tidyverse", "glasso", "qgraph", "Matrix")
suppressMessages(lapply(packages, require, character.only = TRUE))

# get correlation matrix
example.data <- iris[, -5]
cors <- cor(example.data, method = "spearman")

# fit network
zero <- c(4, 3)
pen_mat <- diag(nrow(cors))!=1

wm.example <- EBICglasso(S = cors, n = nrow(example.data), penalizeMatrix = pen_mat, zero = zero)

# Has no edge in 3-4:
wm.example