Possible to use qgraph or weights matrix as input in bootnet?

[mvanaman]

Hi Sacha,

I have an example below where I have estimated a network where some edge weights are constrained to be zero. I accomplished this by manually setting some correlations in my correlation matrix to zero, and then estimating the weight matrix from this correlation matrix before fitting the network with qgraph. Everything seem to work OK.

However, I want to get confidence intervals for this network, but it seems that with the bootnet package, you have to feed the raw data to the estimateNetwork() function which estimates its own new network.

Is there a way to use the bootnet package to get CIs for this constrained model?

For example, by forcing bootnet to accept a weights matrix, correlation matrix, or qgraph object? Or do I need to write some kind of custom function to tell estimateNetwork() to set specific correlations in the correlation matrix to zero before estimating the network? Is that possible?

Thank you for your time,
Matthew

library(psych)
library(qgraph)
library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.1.2
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

# get data
network_data <- bfi %>% select(contains("1") | contains("2"))

# get correlation matrix (for input into unconstrained model)
cors <- cor_auto(network_data, detectOrdinal = TRUE)
#> Variables detected as ordinal: A1; C1; E1; N1; O1; A2; C2; E2; N2; O2
cors %>% head(3)
#>             A1           C1          E1          N1          O1         A2
#> A1 1.000000000  0.003380223  0.11767783  0.18407150 -0.00219489 -0.4084507
#> C1 0.003380223  1.000000000 -0.02859722 -0.08329691  0.21056996  0.1153827
#> E1 0.117677827 -0.028597217  1.00000000  0.02187842 -0.11836930 -0.2395297
#>             C2         E2          N2          O2
#> A1 0.008238801  0.1112004  0.15759982  0.08369205
#> C1 0.483160927 -0.1078603 -0.04690444 -0.14461593
#> E1 0.011974104  0.5153813  0.01297284  0.05267017

# set correlations to zero (for input into constrained model)
cors_constrained <- cors
cors_constrained[6, 1] <- 0
cors_constrained[7, 2] <- 0
cors_constrained[1, 6] <- 0
cors_constrained[2, 7] <- 0

constrained_weights <- EBICglasso(
    S = cors_constrained,
    n = nrow(network_data),
    # default arguments:
    gamma = 0.5,
    penalize.diagonal = FALSE,
    nlambda = 100,
    lambda.min.ratio = 0.01,
    returnAllResults = FALSE,
    checkPD = TRUE,
    countDiagonal = FALSE,
    refit = FALSE,
    threshold = FALSE,
    verbose = TRUE
)
#> Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
#> penalize.diagonal, : 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.

# fit network
qgraph(constrained_weights, edge.labels = TRUE, curveAll = TRUE, title = "constrained")

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

[mvanaman]

It turns out that the example above actually does not work: that is, setting a correlation to zero does not guarantee that its edges will be absent in the estimated weights matrix.

However, with respect to the original question, I finally found an answer to this question in the documentation. It turns out that as long as your function returns a results object i.e. a weights matrix, then you can just use this function with the estimateNetwork() function with the additional argument fun = my_function.

Again: in my case, the actual rationale for modifying this correlation matrix does not make sense. However, for those who might have have a case where the process of estimating a weights matrix is too idiosyncratic to automate with bootnet, perhaps you will find it helpful to know that you can stipulate how you want your weight matrix to be calculated by intervening with that fun argument.

Here is how I accomplished this for my specific situation:

library(psych)
library(qgraph)
library(dplyr)
#> Warning: package 'dplyr' was built under R version 4.1.2
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(bootnet)
#> Loading required package: ggplot2
#> 
#> Attaching package: 'ggplot2'
#> The following objects are masked from 'package:psych':
#> 
#>     %+%, alpha
#> This is bootnet 1.5
#> For questions and issues, please see github.com/SachaEpskamp/bootnet.

# get data
network_data <- bfi %>% select(contains("1") | contains("2"))

# get correlation matrix (for input into unconstrained model)
cors <- cor_auto(network_data, detectOrdinal = TRUE)
#> Variables detected as ordinal: A1; C1; E1; N1; O1; A2; C2; E2; N2; O2

# set correlations to zero (for input into constrained model)
cors_constrained <- cors
cors_constrained[6, 1] <- 0
cors_constrained[7, 2] <- 0
cors_constrained[1, 6] <- 0
cors_constrained[2, 7] <- 0

constrained_weights <- EBICglasso(S = cors_constrained, n = nrow(network_data))
#> Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
#> penalize.diagonal, : 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.

# plot network with qgraph
qgraph(constrained_weights, edge.labels = TRUE, curveAll = TRUE, title = "constrained", layout = "spring")

# achieve same network with bootnet
my_fun <- function(data){
  set.seed(352)
  cor_example <- cor_auto(data)
  n <- nrow(data)
  cors_constrained <- cor_example
  cors_constrained[6, 1] <- 0
  cors_constrained[7, 2] <- 0
  cors_constrained[1, 6] <- 0
  cors_constrained[2, 7] <- 0
  
  # estimate weights matrix
  graph <- EBICglasso(S = cors_constrained, n = n)

  return(graph = graph)
}
estimateNetwork(network_data, fun = my_fun) %>% plot(edge.labels = TRUE, layout = "spring")
#> Variables detected as ordinal: A1; C1; E1; N1; O1; A2; C2; E2; N2; O2
#> Warning in EBICglassoCore(S = S, n = n, gamma = gamma, penalize.diagonal =
#> penalize.diagonal, : 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.

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

[SachaEpskamp]

Hi! Indeed you can use the fun argument for specifying a custom estimation function! Please make sure that the input to EBICglasso is the sample variance-covariance matrix and not a matrix with artificially changed values - such a matrix would no longer lead to the correct Gaussian likelihood of the data making the EBICglasso function not work at all.

[mvanaman]

Thank you for replying, and for the heads up about the likelihood!!