How to choose network analysis estimation for application research

Compare multiple R packages for psychological network analysis

Bayesian

Tutorial

Network Analysis

Author

Jihong Zhang

Published

April 4, 2024

Background

Multiple estimation methods for network analysis have been proposed, from regularized to unregularized, from frequentist to Bayesian approach, from one-step to multiple steps. Isvoranu & Epskamp (2023) provides an throughout illustration of current network estimation methods and simulation study. This post tries to reproduce the procedures and coding illustrated by the paper and compare the results for packages with up-to-date version. Five packages will be used for network modeling: (1) qgraph, (2)psychonetrics, (3)MGM, (4)BGGM, (5)GGMnonreg. Note that the tutorial of each R package is not the scope of this post. The specific variants of estimation is described in Figure 1. Please see here for BGGM R package for more detailed example.

Isvoranu, A.-M., & Epskamp, S. (2023). Which estimation method to choose in network psychometrics? Deriving guidelines for applied researchers. Psychological Methods, 28(4), 925–946. https://doi.org/10.1037/met0000439

Figure 1: Estimation methods used in Isvoranu and Epskamp (2023)

Data Generation

For the sake of simplicity, I will only pick one condition from the simulation settings. To simulate data, I used the same code from the paper. The function is a wrapper of bootnet::ggmGenerator function

dataGenerator() function

library("bootnet")
dataGenerator <- function(
  trueNet,  # true network models: DASS21, PTSD, BFI
  sampleSize, # sample size
  data = c("normal","skewed","uniform ordered","skewed ordered"), # type of data
  nLevels = 4,
  missing = 0
){
  data <- match.arg(data)
  nNode <- ncol(trueNet)
  
  if (data == "normal" || data == "skewed"){
    # Generator function:
    gen <- ggmGenerator()
    
    # Generate data:
    Data <- gen(sampleSize,trueNet)
    
    # Generate replication data:
    Data2 <- gen(sampleSize,trueNet)
    
    if (data == "skewed"){ ## exponential transformation to make data skewed
      for (i in 1:ncol(trueNet)){
        Data[,i] <- exp(Data[,i])
        Data2[,i] <- exp(Data2[,i])
      }
    }
  
  } else {
    # Skew factor:
    skewFactor <- switch(data,
                         "uniform ordered" = 1,
                         "skewed ordered" = 2,
                         "very skewed ordered" = 4)
    
    # Generator function:
    gen <- ggmGenerator(ordinal = TRUE, nLevels = nLevels)
    
    # Make thresholds:
    thresholds <- lapply(seq_len(nNode),function(x)qnorm(seq(0,1,length=nLevels + 1)[-c(1,nLevels+1)]^(1/skewFactor)))
    
    # Generate data:
    Data <- gen(sampleSize, list(
      graph = trueNet,
      thresholds = thresholds
    ))
    
    # Generate replication data:
    Data2 <- gen(sampleSize, list(
      graph = trueNet,
      thresholds = thresholds
    ))
  }
  
  # Add missings:
  if (missing > 0){
    for (i in 1:ncol(Data)){
      Data[runif(sampleSize) < missing,i] <- NA
      Data2[runif(sampleSize) < missing,i] <- NA
    }
  }
  
  return(list(
    data1 = Data,
    data2 = Data2
  ))
}

BFI contains 26 columns: first columns are the precision matrix with 25 \times 25 and second columns as clusters. We will use the structure of BFI for the simulation.

Click to see the code

library(here)
bfi_truenetwork <- read.csv(here("notes", "2024-04-04-Network-Estimation-Methods", 
                                 "weights matrices", "BFI.csv"))
bfi_truenetwork <- bfi_truenetwork[,1:25]
head(psychTools::bfi)

      A1 A2 A3 A4 A5 C1 C2 C3 C4 C5 E1 E2 E3 E4 E5 N1 N2 N3 N4 N5 O1 O2 O3 O4
61617  2  4  3  4  4  2  3  3  4  4  3  3  3  4  4  3  4  2  2  3  3  6  3  4
61618  2  4  5  2  5  5  4  4  3  4  1  1  6  4  3  3  3  3  5  5  4  2  4  3
61620  5  4  5  4  4  4  5  4  2  5  2  4  4  4  5  4  5  4  2  3  4  2  5  5
61621  4  4  6  5  5  4  4  3  5  5  5  3  4  4  4  2  5  2  4  1  3  3  4  3
61622  2  3  3  4  5  4  4  5  3  2  2  2  5  4  5  2  3  4  4  3  3  3  4  3
61623  6  6  5  6  5  6  6  6  1  3  2  1  6  5  6  3  5  2  2  3  4  3  5  6
      O5 gender education age
61617  3      1        NA  16
61618  3      2        NA  18
61620  2      2        NA  17
61621  5      2        NA  17
61622  3      1        NA  17
61623  1      2         3  21

To generate data from bootnet::ggmGenerator(), we can generate a function gen from it with first argument as sample size and second argument as partial correlation matrix.

Click to see the code

gen <- bootnet::ggmGenerator()
data <- gen(2800, as.matrix(bfi_truenetwork))

Measures

Overview of Centrality Measures

Centrality measures are used as one characteristic of network model to summarize the importance of nodes/items given a network. Some ideas of centrality measures come from social network. Some may not make much sense in psychometric network. I also listed the centrality measures that have been used in literatures of psychological or psychometric modeling. Following the terminology of Christensen (2018), centrality measures include: (1) betweenness centrality (BC); (2) the randomized shortest paths betweeness centrality (RSPBC); (3) closeness centrality (CC); (4)node degree (ND); (5) node strength (NS); (6) expected influence (EI); (7) engenvector centrality (EC); (8) Eccentricity (k; Hage & Harary (1995)); (9) hybrid centrality (HC; Christensen et al. (2018a)).

Christensen, A. P. (2018). NetworkToolbox: Methods and measures for brain, cognitive, and psychometric network analysis in r. R J., 10(2), 422. https://files.osf.io/v1/resources/6kmav/providers/osfstorage/5aec6ad0b80ad4001054a820?action=download&version=1&displayName=NetworkToolbox_Preprint_psyarxiv-2018-05-04T14:14:40.992Z.pdf&direct

Hage, P., & Harary, F. (1995). Eccentricity and centrality in networks. Social Networks, 17(1), 57–63. https://doi.org/10.1016/0378-8733(94)00248-9

Christensen, A. P., Kenett, Y. N., Aste, T., Silvia, P. J., & Kwapil, T. R. (2018a). Network structure of the Wisconsin Schizotypy ScalesShort Forms: Examining psychometric network filtering approaches. Behavior Research Methods, 50(6), 2531–2550. https://doi.org/10.3758/s13428-018-1032-9

They can be grouped into three categories:

Type 1—Number of path relevant metrics: BC measures how often a node is on the shortest path from one node to another. RSPBC is a adjusted BC measure which measures how oftern a node is on the random shortest path from one to another; CC measures the average number of paths from one node to all other nodes in the network. ND measures how many connections are connected to each node. Eccentricity (k) measures the maximum “distance” between one node with other nodes with lower values suggests higher centrality, where one node’s distance is the length of a shortest path between this node to other nodes.

Type 2—Path strength relevant metrics: NS measures the sum of the absolute values of edge weights connected to a single node. EI measures the sum of the positive or negative values of edge weights connected to a single node.

Type 3—Composite of multiple centrality measures: HC measures nodes on the basis of their centrality values across multiple measures of centrality (BC, LC, k, EC, and NS) which describes highly central nodes with large values and highly peripheral nodes with small values (Christensen et al., 2018b).

Christensen, A. P., Kenett, Y. N., Aste, T., Silvia, P. J., & Kwapil, T. R. (2018b). Network structure of the Wisconsin Schizotypy ScalesShort Forms: Examining psychometric network filtering approaches. Behavior Research Methods, 50(6), 2531–2550. https://doi.org/10.3758/s13428-018-1032-9

Hallquist, M. N., Wright, A. G. C., & Molenaar, P. C. M. (2021). Problems with centrality measures in psychopathology symptom networks: Why network psychometrics cannot escape psychometric theory. Multivariate Behavioral Research, 56(2), 199–223. https://doi.org/10.1080/00273171.2019.1640103

Neal, Z. P., Forbes, M. K., Neal, J. W., Brusco, M. J., Krueger, R., Markon, K., Steinley, D., Wasserman, S., & Wright, A. G. C. (2022). Critiques of network analysis of multivariate data in psychological science. Nature Reviews Methods Primers, 2(1), 1–2. https://doi.org/10.1038/s43586-022-00177-9

Bringmann, L. F., Elmer, T., Epskamp, S., Krause, R. W., Schoch, D., Wichers, M., Wigman, J. T. W., & Snippe, E. (2019). What do centrality measures measure in psychological networks? Journal of Abnormal Psychology, 128(8), 892–903. https://doi.org/10.1037/abn0000446

Some literature of network psychometrics criticized the use of Type 1 centrality measures and some Type II centrality measures (Hallquist et al., 2021; Neal et al., 2022). For example, Hallquist et al. (2021) argued that some Type I centrality measures such as BC and CC derive from the concept of distance, which builds on the physical nature of many traditional graph theory applications, including railways and compute networks. In association networks, the idea of network distance does not have physical reference. Bringmann et al. (2019) also think distance-based centrality like closeness centrality (CC) and betweenness centrality (BC) do not have “shortest path”, or “node exchangeability” assumptions in psychological networks so they are not suitable in the context of psychological network. In addition, ND (degree), CC (closeness), and BC (betweenness) do not take negative values of edge weights into account that may lose important information in the context of psychology.

Network level accuracy metrics

Structure accuracy measures such as sensitivity, specificity, and precision investigate if an true edge is include or not or an false edge is include or not. Several edge weight accuracy measures include: (1) correlation between absolute values of estimated edge weights and true edge weights, (2) average absolute bias between estimated edge weights and true weights, (3) average absolute bias for true edges.

Click to see the code

cor0 <- function(x,y,...){
  if (sum(!is.na(x)) < 2 || sum(!is.na(y)) < 2 || sd(x,na.rm=TRUE)==0 | sd(y,na.rm=TRUE) == 0){
    return(0)
  } else {
    return(cor(x,y,...))
  }
}

bias <- function(x,y) mean(abs(x-y),na.rm=TRUE)


### Inner function:
comparison_metrics <- function(real, est, name = "full"){
  
  # Output list:
  out <- list()
  
  # True positives:
  TruePos <- sum(est != 0 &  real != 0)
  
  # False pos:
  FalsePos <- sum(est != 0 & real == 0)
  
  # True Neg:
  TrueNeg <- sum(est == 0 & real == 0)
  
  # False Neg:
  FalseNeg <- sum(est == 0 & real != 0)
  
  # Sensitivity:
  out$sensitivity <- TruePos / (TruePos + FalseNeg)
  
  # Sensitivity top 50%:
  top50 <- which(abs(real) > median(abs(real[real!=0])))
  out[["sensitivity_top50"]] <- sum(est[top50]!=0 & real[top50] != 0) / sum(real[top50] != 0)
  
  # Sensitivity top 25%:
  top25 <- which(abs(real) > quantile(abs(real[real!=0]), 0.75))
  out[["sensitivity_top25"]] <- sum(est[top25]!=0 & real[top25] != 0) / sum(real[top25] != 0)
  
  # Sensitivity top 10%:
  top10 <- which(abs(real) > quantile(abs(real[real!=0]), 0.90))
  out[["sensitivity_top10"]] <- sum(est[top10]!=0 & real[top10] != 0) / sum(real[top10] != 0)
  
  # Specificity:
  out$specificity <- TrueNeg / (TrueNeg + FalsePos)
  
  # Precision (1 - FDR):
  out$precision <- TruePos / (FalsePos + TruePos)
  
  # precision top 50% (of estimated edges):
  top50 <- which(abs(est) > median(abs(est[est!=0])))
  out[["precision_top50"]] <- sum(est[top50]!=0 & real[top50] != 0) / sum(est[top50] != 0)
  
  # precision top 25%:
  top25 <- which(abs(est) > quantile(abs(est[est!=0]), 0.75))
  out[["precision_top25"]] <- sum(est[top25]!=0 & real[top25] != 0) / sum(est[top25] != 0)
  
  # precision top 10%:
  top10 <- which(abs(est) > quantile(abs(est[est!=0]), 0.90))
  out[["precision_top10"]] <- sum(est[top10]!=0 & real[top10] != 0) / sum(est[top10] != 0)
  
  # Signed sensitivity:
  TruePos_signed <- sum(est != 0 &  real != 0 & sign(est) == sign(real))
  out$sensitivity_signed <- TruePos_signed / (TruePos + FalseNeg)
  
  # Correlation:
  out$correlation <- cor0(est,real)
  
  # Correlation between absolute edges:
  out$abs_cor <- cor0(abs(est),abs(real))
  
  #
  out$bias <- bias(est,real)
  
  ## Some measures for true edges only:
  if (TruePos > 0){
    
    trueEdges <- est != 0 & real != 0
    
    out$bias_true_edges <- bias(est[trueEdges],real[trueEdges])
    out$abs_cor_true_edges <- cor0(abs(est[trueEdges]),abs(real[trueEdges]))
  } else {
    out$bias_true_edges <- NA
    out$abs_cor_true_edges <- NA
  }
  
  out$truePos <- TruePos
  out$falsePos <- FalsePos
  out$trueNeg <- TrueNeg
  out$falseNeg <- FalseNeg
  
  # Mean absolute weight false positives:
  false_edges <- (est != 0 &  real == 0) | (est != 0 & real != 0 & sign(est) != sign(real) )
  out$mean_false_edge_weight <- mean(abs(est[false_edges]))
  out$SD_false_edge_weight <- sd(abs(est[false_edges]))
  
  # Fading:
  out$maxfade_false_edge <- max(abs(est[false_edges])) / max(abs(est))
  out$meanfade_false_edge <- mean(abs(est[false_edges])) / max(abs(est))
  
  
  # Set naname
  if (name != ""){
    names(out) <- paste0(names(out),"_",name)  
  }
  out
}

Estimation

Method 1 - EBICglasso in `bootnet`

This method uses the bootnet::estimateNetwork function.

EBICtuning <- 0.5 sets up the hyperparameter (\gamma) as .5 that controls how much the EBIC prefers simpler models (fewer edges), which by default is .5.

\text{EBIC} =-2L+E(\log(N))+4\gamma E(\log(P))

Another important setting is the tuning parameter (\lambda) of EBICglasso that controls the sparsity level in the penalized likelihood function as:

\log \det(\boldsymbol{K}) - \text{trace}(\boldsymbol{SK})-\lambda \Sigma|\kappa_{ij}|

where \boldsymbol{S} represents the sample variance-covariance matrix and \boldsymbol{K} represents the precision matrix that lasso aims to estimate. Overall, the regularization limits the sum of absolute partial correlation coefficients.

sampleSize="pairwise_average" will set the sample size to the average of sample sizes used for each individual correlation in EBICglasso.

Figure 2 is the estimated network structure with BFI-25 data.

Method 1: EBICglasso

library("qgraph")
library("bootnet")
sampleAdjust  <-  "pairwise_average"
EBICtuning <- 0.5
transformation <- "polychoric/categorical" # EBICglasso use cor_auto
res_m1 <- estimateNetwork(data, # input as polychoric correlation
                          default = "EBICglasso",
                          sampleSize = sampleAdjust,
                          tuning = EBICtuning,
                          corMethod = ifelse(transformation == "polychoric/categorical",
                                             "cor_auto",
                                             "cor"))
    
estnet_m1 <- res_m1$graph

qgraph(estnet_m1)

Figure 2: Method 1: Network Plot for BFI-25 with EBICglasso

Click to see the code

performance <- function(network) {
  res_comparison <- comparison_metrics(real = as.matrix(bfi_truenetwork), est = as.matrix(network))
  c(
    sensitivity = res_comparison$sensitivity_full,
    specificity = res_comparison$specificity_full,
    precision = res_comparison$precision_full,
    
    abs_corr = res_comparison$abs_cor_true_edges_full, # The Pearson correlation between the absolute edge weights
    bias = res_comparison$bias_full, # The average absolute deviation
    bias_true = res_comparison$bias_true_edges_full # The average absolute deviation between the true edge weight and the estimated edge weight
  )
}
performance(estnet_m1)

sensitivity specificity   precision    abs_corr        bias   bias_true 
0.982142857 0.840399002 0.774647887 0.963820333 0.009005281 0.021625848

Method 2 - ggmModSelect in `bootnet`

The ggmModSelect algorithm is a non-regularized method.

Method 2: ggmModSelect

# With stepwise
res_m2_1 <- estimateNetwork(data, default = "ggmModSelect",
                           sampleSize = sampleAdjust, stepwise = TRUE,
                           tuning = EBICtuning,
                           corMethod = ifelse(transformation == "polychoric/categorical",
                                              "cor_auto",
                                              "cor"))

estnet_m2_1 <- res_m2_1$graph
qgraph(estnet_m2_1)

Figure 3: Method 2.1: Network Plot for BFI-25 with ggmModSelect and Stepwise

Method 2: ggmModSelect

# Without stepwise
res_m2_2 <- estimateNetwork(data, default = "ggmModSelect",
                           sampleSize = sampleAdjust, stepwise = FALSE,
                           tuning = EBICtuning,
                           corMethod = ifelse(transformation == "polychoric/categorical",
                                              "cor_auto",
                                              "cor"))
estnet_m2_2 <- res_m2_2$graph
qgraph(estnet_m2_2)

Figure 4: Method 2.1: Network Plot for BFI-25 with ggmModSelect and No Stepwise

Click to see the code

performance_compr <- Reduce(rbind, lapply(list(estnet_m1, estnet_m2_1, estnet_m2_2), performance))
rownames(performance_compr) <- c("EBICglasso", "ggmModSelect_stepwise", "ggmModSelect_nostepwise")
kableExtra::kable(performance_compr)

	sensitivity	specificity	precision	abs_corr	bias	bias_true
EBICglasso	0.9821429	0.8403990	0.7746479	0.9638203	0.0090053	0.0216258
ggmModSelect_stepwise	0.8839286	1.0000000	1.0000000	0.9588360	0.0081832	0.0170543
ggmModSelect_nostepwise	0.9375000	0.9002494	0.8400000	0.9669727	0.0078012	0.0148425

Method 3 - FIML in `psychonetrics`

Method 3: FIML

library(psychonetrics)
library(magrittr)
# With stepwise
res_m3_1 <- ggm(data, estimator = "FIML", standardize = "z") %>% 
  prune(alpha = .01, recursive = FALSE) 
    
estnet_m3_1 <- getmatrix(res_m3_1, "omega")

qgraph(estnet_m3_1)

Figure 5: Method 3: Network Plot for BFI-25 with FIML and Prue

Method 3: FIML

# With stepwise
res_m3_2 <- ggm(data, estimator = "FIML", standardize = "z") %>% 
  prune(alpha = .01, recursive = FALSE) |> 
  stepup(alpha = .01)
    
estnet_m3_2 <- getmatrix(res_m3_2, "omega")

qgraph(estnet_m3_2)

Figure 6: Method 3: Network Plot for BFI-25 with FIML and Stepup

Click to see the code

performance_compr <- Reduce(rbind, lapply(list(estnet_m1, estnet_m2_1, estnet_m2_2, 
                                               estnet_m3_1, estnet_m3_2), performance))
rownames(performance_compr) <- c("EBICglasso", "ggmModSelect_stepwise", "ggmModSelect_nostepwise",
                                 "FIML_prune", "FIML_stepup")
kableExtra::kable(performance_compr, digits = 3)

	sensitivity	specificity	precision	abs_corr	bias	bias_true
EBICglasso	0.982	0.84	0.775	0.964	0.009	0.022
ggmModSelect_stepwise	0.884	1.00	1.000	0.959	0.008	0.017
ggmModSelect_nostepwise	0.938	0.90	0.840	0.967	0.008	0.015
FIML_prune	0.929	0.99	0.981	0.966	0.007	0.015
FIML_stepup	0.982	0.99	0.982	0.972	0.006	0.014

Method 4 - WLS in `psychonetrics`

Method 4: WLS

# With stepwise
res_m4_1 <- ggm(data, estimator = "WLS", standardize = "z") %>% 
  prune(alpha = .01, recursive = FALSE) 
    
estnet_m4_1 <- getmatrix(res_m4_1, "omega")

qgraph(estnet_m4_1)

Figure 7: Method 4: Network Plot for BFI-25 with WLS and prune

Method 4: WLS

res_m4_2 <- ggm(data, estimator = "WLS", standardize = "z") %>% 
  prune(alpha = .01, recursive = FALSE) |> 
  modelsearch(prunealpha = .01, addalpha = .01)
    
estnet_m4_2 <- getmatrix(res_m4_2, "omega")

qgraph(estnet_m4_2)

Figure 8: Method 4: Network Plot for BFI-25 with WLS, prune, and modelsearch

Click to see the code

performance_compr <- Reduce(rbind, lapply(list(estnet_m1, estnet_m2_1, estnet_m2_2, 
                                               estnet_m3_1, estnet_m3_2,
                                               estnet_m4_1, estnet_m4_2), performance))
rownames(performance_compr) <- c("EBICglasso", "ggmModSelect_stepwise", "ggmModSelect_nostepwise",
                                 "FIML_prune", "FIML_stepup",
                                 "WLS_prune", "WLS_prune_modelsearch")
kableExtra::kable(performance_compr, digits = 3)

	sensitivity	specificity	precision	abs_corr	bias	bias_true
EBICglasso	0.982	0.840	0.775	0.964	0.009	0.022
ggmModSelect_stepwise	0.884	1.000	1.000	0.959	0.008	0.017
ggmModSelect_nostepwise	0.938	0.900	0.840	0.967	0.008	0.015
FIML_prune	0.929	0.990	0.981	0.966	0.007	0.015
FIML_stepup	0.982	0.990	0.982	0.972	0.006	0.014
WLS_prune	0.938	0.995	0.991	0.973	0.006	0.014
WLS_prune_modelsearch	0.938	0.995	0.991	0.973	0.006	0.014

Method 5 - mgm in `mgm`

mgm package is used for estimating mixed graphical models. Node type can be “g” for Gaussian or “c” for categorical. For continuous variables, set level to 1 and type to g.

Method 5: mgm

library(mgm)
## mgm with CV
res_m5_1 <- mgm(na.omit(data), type = rep("g", ncol(data)), level = rep(1, ncol(data)), 
           lambdaFolds = 20,
           lambdaSel = "CV", lambdaGam = EBICtuning, 
           pbar = FALSE)

Note that the sign of parameter estimates is stored separately; see ?mgm

Method 5: mgm

getmatrix_mgm <- function(res) {
  sign = ifelse(is.na(res$pairwise$signs), 1, res$pairwise$signs)
  sign * res$pairwise$wadj
}
estnet_m5_1 <- getmatrix_mgm(res_m5_1)
qgraph(estnet_m5_1)

Figure 9: Method 5: Network Plot for BFI-25 with mgm and CV

Method 5: mgm

## mgm with EBIC
res_m5_2 <- mgm(na.omit(data), type = rep("g", ncol(data)), level = rep(1, ncol(data)), 
           lambdaFolds = 20,
           lambdaSel = "EBIC", lambdaGam = EBICtuning,
           pbar = FALSE)

Note that the sign of parameter estimates is stored separately; see ?mgm

Method 5: mgm

estnet_m5_2 <- getmatrix_mgm(res_m5_2)
qgraph(estnet_m5_2)

Figure 10: Method 5: Network Plot for BFI-25 with mgm and EBIC

Click to see the code

performance_compr <- Reduce(rbind, lapply(list(estnet_m1, estnet_m2_1, estnet_m2_2, 
                                               estnet_m3_1, estnet_m3_2,
                                               estnet_m4_1, estnet_m4_2,
                                               estnet_m5_1, estnet_m5_2), performance))
rownames(performance_compr) <- c("EBICglasso", "ggmModSelect_stepwise", "ggmModSelect_nostepwise",
                                 "FIML_prune", "FIML_stepup",
                                 "WLS_prune", "WLS_prune_modelsearch",
                                 "mgm_CV", "mgm_EBIC")
kableExtra::kable(performance_compr, digits = 3)

	sensitivity	specificity	precision	abs_corr	bias	bias_true
EBICglasso	0.982	0.840	0.775	0.964	0.009	0.022
ggmModSelect_stepwise	0.884	1.000	1.000	0.959	0.008	0.017
ggmModSelect_nostepwise	0.938	0.900	0.840	0.967	0.008	0.015
FIML_prune	0.929	0.990	0.981	0.966	0.007	0.015
FIML_stepup	0.982	0.990	0.982	0.972	0.006	0.014
WLS_prune	0.938	0.995	0.991	0.973	0.006	0.014
WLS_prune_modelsearch	0.938	0.995	0.991	0.973	0.006	0.014
mgm_CV	0.964	0.865	0.800	0.973	0.009	0.017
mgm_EBIC	0.821	0.980	0.958	0.965	0.012	0.026

Method 6 - BGGM in `BGGM`

Method 6: BGGM

res_m6_1 <- BGGM::explore(data, type = "continuous") |> 
  BGGM:::select.explore(BF_cut = 3)
estnet_m6_1 <- res_m6_1$pcor_mat_zero
qgraph(estnet_m6_1)

Figure 11: Method 6: Network Plot for BFI-25 with BGGM and explore

Method 6: BGGM

res_m6_2 <- BGGM::explore(data, type = "continuous") |> 
  BGGM:::select.estimate(cred = 0.95)
estnet_m6_2 <- res_m6_2$pcor_adj
qgraph(estnet_m6_2)

Figure 12: Method 6: Network Plot for BFI-25 with BGGM and estimate

Click to see the code

performance_compr <- Reduce(rbind, lapply(list(estnet_m1, estnet_m2_1, estnet_m2_2, 
                                               estnet_m3_1, estnet_m3_2,
                                               estnet_m4_1, estnet_m4_2,
                                               estnet_m5_1, estnet_m5_2,
                                               estnet_m6_1, estnet_m6_2), performance))
rownames(performance_compr) <- c("EBICglasso", "ggmModSelect_stepwise", "ggmModSelect_nostepwise",
                                 "FIML_prune", "FIML_stepup",
                                 "WLS_prune", "WLS_prune_modelsearch",
                                 "mgm_CV", "mgm_EBIC",
                                 "BGGM_explore", "BGGM_estimate")
kableExtra::kable(performance_compr, digits = 3)

	sensitivity	specificity	precision	abs_corr	bias	bias_true
EBICglasso	0.982	0.840	0.775	0.964	0.009	0.022
ggmModSelect_stepwise	0.884	1.000	1.000	0.959	0.008	0.017
ggmModSelect_nostepwise	0.938	0.900	0.840	0.967	0.008	0.015
FIML_prune	0.929	0.990	0.981	0.966	0.007	0.015
FIML_stepup	0.982	0.990	0.982	0.972	0.006	0.014
WLS_prune	0.938	0.995	0.991	0.973	0.006	0.014
WLS_prune_modelsearch	0.938	0.995	0.991	0.973	0.006	0.014
mgm_CV	0.964	0.865	0.800	0.973	0.009	0.017
mgm_EBIC	0.821	0.980	0.958	0.965	0.012	0.026
BGGM_explore	0.920	1.000	1.000	0.970	0.007	0.015
BGGM_estimate	0.991	0.980	0.965	0.969	0.006	0.015

Summary

Click to see the code

library(dplyr)
library(tidyr)
library(ggplot2)
performance_tbl <- performance_compr |> 
  as.data.frame() |> 
  mutate(method = rownames(performance_compr)) |> 
  mutate(method = factor(method, levels = rownames(performance_compr)))
ggplot(performance_tbl) +
  geom_point(aes(x = sensitivity, y = specificity, col = method), size = 3) +
  theme_bw() +
  theme(legend.position = 'bottom')

Figure 13: Specificity and Sensitivity Balance

FLML_prune and BGGM_estimation seem to perform best among all methods regarding balance between sensitivity and specificity of network weights.

Click to see the code

ggplot(performance_tbl) +
  geom_col(aes(y = forcats::fct_reorder(method, precision), x = precision, fill = method)) +
  theme_bw() +
  labs(y = "") +
  theme(legend.position = 'bottom')

ggmMoSelect with the stepwise procedure appears to have highest precision, followed by BGGM explore.

Background

Data Generation

Measures

Overview of Centrality Measures

Network level accuracy metrics

Estimation

Method 1 - EBICglasso in bootnet

Method 2 - ggmModSelect in bootnet

Method 3 - FIML in psychonetrics

Method 4 - WLS in psychonetrics

Method 5 - mgm in mgm

Method 6 - BGGM in BGGM

Summary

Method 1 - EBICglasso in `bootnet`

Method 2 - ggmModSelect in `bootnet`

Method 3 - FIML in `psychonetrics`

Method 4 - WLS in `psychonetrics`

Method 5 - mgm in `mgm`

Method 6 - BGGM in `BGGM`