Fit a network where some edges are constrained to zero?

[mvanaman]

Say I had a theoretical reason to compare two networks on e.g., the iris data: a constrained network where the edge between Sepal.Length and Petal.Width is zero while the rest are freely estimated vs an unconstrained network where all edges are freely estimated. I'm thinking of it as analogous to a nested model comparison in regression.

Is it possible to fit a constrained network of this sort in qgraph?

[mvanaman]

EDIT: it turned out that in this example, everything worked perfectly. However, it was a fluke: trying it out on some other data sets in with some different seeds show that setting the correlation to zero does not guarantee that the edge will be absent in the weights matrix. So this is still an open question.

I figured this out:

1. First, you create two correlation matrices: one - correlations are estimated, and two - some correlations are manually set to zero by the user.
2. Next, you estimate the weights matrices separately for each correlation matrix.
3. Last, you fit the network using qgraph.

Example below doesn't use the iris data set because iris ended up having too few variables for this. So I'm using the bfi data from the psych package instead. In the constrained model, one example edgeweight constrained to zero is the edge from A1 to A2. The rest of the edges are estimated using EBICLasso, which automatically selects the best-fitting model:

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_unconstrained <- cor_auto(network_data, detectOrdinal = TRUE)
#> Variables detected as ordinal: A1; C1; E1; N1; O1; A2; C2; E2; N2; O2
cors_unconstrained %>% 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_unconstrained
cors_constrained[6, 1] <- 0
cors_constrained[7, 2] <- 0
cors_constrained[1, 6] <- 0
cors_constrained[2, 7] <- 0
cors_constrained %>% head(3)
#>             A1           C1          E1          N1          O1         A2
#> A1 1.000000000  0.003380223  0.11767783  0.18407150 -0.00219489  0.0000000
#> 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.000000000 -0.1078603 -0.04690444 -0.14461593
#> E1 0.011974104  0.5153813  0.01297284  0.05267017

# calculate weights matrices (for input into graph)
unconstrained_weights <- EBICglasso(
    S = cors_unconstrained,
    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.
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(unconstrained_weights, edge.labels = TRUE, curveAll = TRUE, title = "unconstrained")

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

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