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

2 Simulation Study 1

2.1 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")
library("qgraph")

mycolors = c("#4682B4", "#B4464B", "#752021",  
             "#1B9E77", "#66A61E", "#D95F02", "#7570B3",
             "#E7298A", "#E75981", "#B4AF46", "#B4726A")

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("posts", "2024", "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))

2.2 Measures

2.2.1 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.

2.2.2 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
}

2.3 Estimation

This method uses the bootnet::estimateNetwork function.

This method makes used of polychoric correlation as input, graphical LASSO as regularlization, and EBIC as model selection method. The polychoric correlation is a technique for estimating the correlation between two hypothesised normally distributed continuous latent variables, from two observed ordinal variables.

The figure below plot the relationships between two continous variables vs. relationships between two ordinal variables (using floor function). It appears that they did not deviate a lot given the values of pearson correlation (-.32) and polychoric correlation (-.3).

Click to see the code

dt <- as.data.frame(data[, 1:2])
dt_ordinal <- as.data.frame(floor(data[, 1:2]))

 # polychoric correlation
pearson_corr <- round(cor_auto(dt[, 1:2])[1, 2], 2)
polychoric_corr <- round(cor_auto(floor(data[, 1:2]))[1, 2], 2)
ggplot(dt) +
  geom_point(aes(x = V1, y = V2), data = dt, color = mycolors[1], alpha = .4) +
  geom_smooth(aes(x = V1, y = V2), data = dt, color = mycolors[1], method = "lm", se = FALSE) +
  geom_point(aes(x = V1, y = V2), data = dt_ordinal, color = mycolors[2]) +
  geom_smooth(aes(x = V1, y = V2), data = dt_ordinal, color = mycolors[2], method = "lm", se = FALSE) +
  geom_text(aes(x = 2, y = 2, label = paste0("Pearson: ", pearson_corr))) +
  geom_text(aes(x = 2, y = 1.7, label = paste0("Polychoric: ", polychoric_corr))) +
  theme_classic()

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+\log(N)+4\gamma \log(P)$

$\text{AIC} = -2L + 2P$

$\text{BIC} = -2L + P\log(N)$

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
summary(res_m1)


=== Estimated network ===
Number of nodes: 25 
Number of non-zero edges: 137 / 300 
Mean weight: 0.01972212 
Network stored in object$graph 
 
Default set used: EBICglasso 
 
Use plot(object) to plot estimated network 
Use bootnet(object) to bootstrap edge weights and centrality indices 

Relevant references:

    Friedman, J. H., Hastie, T., & Tibshirani, R. (2008). Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9 (3), 432-441.
    Foygel, R., & Drton, M. (2010). Extended Bayesian information criteria for Gaussian graphical models. 
    Friedman, J. H., Hastie, T., & Tibshirani, R. (2014). glasso: Graphical lasso estimation of gaussian graphical models. Retrieved from https://CRAN.R-project.org/package=glasso
    Epskamp, S., Cramer, A., Waldorp, L., Schmittmann, V. D., & Borsboom, D. (2012). qgraph: Network visualizations of relationships in psychometric data. Journal of Statistical Software, 48 (1), 1-18.
    Epskamp, S., Borsboom, D., & Fried, E. I. (2016). Estimating psychological networks and their accuracy: a tutorial paper. arXiv preprint, arXiv:1604.08462.

Method 1: EBICglasso

str(res_m1$results)

List of 5
 $ results:List of 5
  ..$ w      : num [1:25, 1:25, 1:100] 1 -0.342 -0.285 -0.139 -0.177 ...
  ..$ wi     : num [1:25, 1:25, 1:100] 1.22117 0.32844 0.19109 -0.00906 0 ...
  ..$ approx : logi FALSE
  ..$ rholist: num [1:100] 0.00702 0.00736 0.00771 0.00808 0.00846 ...
  ..$ errflag: int [1:100] 0 0 0 0 0 0 0 0 0 0 ...
 $ ebic   : num [1:100] 53234 53210 53158 53091 53054 ...
 $ optnet : num [1:25, 1:25] 0 -0.215 -0.118 0 0 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:25] "V1" "V2" "V3" "V4" ...
  .. ..$ : chr [1:25] "V1" "V2" "V3" "V4" ...
 $ lambda : num [1:100] 0.00702 0.00736 0.00771 0.00808 0.00846 ...
 $ optwi  : num [1:25, 1:25] 1.165 0.287 0.161 0 0 ...

Method 1: EBICglasso

EBIC = mean(res_m1$results$ebic) # EBIC of the final model
(AIC = EBIC - 2 * log(nrow(data))) # AIC =  52287.32

[1] 56560.88

Method 1: EBICglasso

qgraph(estnet_m1)

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

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

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")
fit(res_m3_1)

           Measure     Value
              logl -89477.10
 unrestricted.logl -89347.11
     baseline.logl -99313.20
              nvar        25
              nobs       350
              npar       154
                df       196
         objective     17.97
             chisq    259.96
            pvalue    0.0015
    baseline.chisq  19932.16
     baseline.npar        50
       baseline.df       300
   baseline.pvalue       ~ 0
               nfi      0.99
              pnfi      0.64
               tli       1.0
              nnfi       1.0
               rfi      0.98
               ifi       1.0
               rni       1.0
               cfi       1.0
             rmsea     0.011
    rmsea.ci.lower    0.0069
    rmsea.ci.upper     0.014
      rmsea.pvalue         1
            aic.ll 179262.19
           aic.ll2 179280.24
             aic.x   -132.04
            aic.x2    567.96
               bic 180176.55
              bic2 179687.24
           ebic.25 180672.25
            ebic.5 181167.96
           ebic.75 181564.53
             ebic1 182159.38

Method 3: FIML

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

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

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

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

2.4 Performance Comparison

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

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	1.000	0.875	0.818	0.958	0.009	0.023
ggmModSelect_stepwise	0.902	1.000	1.000	0.966	0.007	0.015
ggmModSelect_nostepwise	0.964	0.920	0.871	0.969	0.007	0.014
FIML_prune	0.929	1.000	1.000	0.967	0.006	0.014
FIML_stepup	0.973	1.000	1.000	0.973	0.005	0.013
WLS_prune	0.929	1.000	1.000	0.963	0.007	0.015
WLS_prune_modelsearch	0.929	1.000	1.000	0.963	0.007	0.015
mgm_CV	0.991	0.925	0.881	0.969	0.007	0.016
mgm_EBIC	0.911	0.995	0.990	0.964	0.010	0.023
BGGM_explore	0.893	1.000	1.000	0.969	0.007	0.015
BGGM_estimate	0.973	0.955	0.924	0.968	0.007	0.015

2.5 Summary

Click to see the code

library(dplyr)
library(tidyr)
library(ggplot2)
mycolors = c("#4682B4", "#B4464B", "#752021",  
             "#1B9E77", "#66A61E", "#D95F02", "#7570B3",
             "#E7298A", "#E75981", "#B4AF46", "#B4726A")
performance_tbl <- performance_compr |> 
  as.data.frame() |> 
  mutate(method = rownames(performance_compr)) |> 
  mutate(method = factor(method, levels = rownames(performance_compr))) |> 
  filter(
    method != "WLS_prune_modelsearch"
  )
ggplot(performance_tbl) +
  geom_point(aes(x = sensitivity, y = specificity, col = method), size = 3) +
  geom_text(aes(x = sensitivity, y = specificity, label = method), size = 3,
            nudge_x = .005, nudge_y = .005) +
  geom_smooth(aes(x = sensitivity, y = specificity), se = FALSE, method = lm) +
  theme_bw() +
  scale_color_manual(values = mycolors) +
  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.