High Dimensional Data: Must Read!!

Introduction

GGMncv is set up for low-dimensional settings, that is, when the number of observations (n) is much greater than the number of nodes (p). This is perhaps not typical in the Gaussian graphical modeling literature, and is a direct result of my (the author of GGMncv) field encountering low-dimensional data most often (see for example Williams et al. 2019; Williams and Rast 2020). As a result, the defaults are honed in for low-dimensional data!

Of course, GGMncv can readily be used for high-dimensional data. In what follows, I highlight several issues that may arise, in addition to solutions to overcome those issues.

Failure Altogether

By default, GGMncv uses the sample based (inverse) covariance matrix for the initial values, which is needed for employing nonconvex regularization. When p > n, GGMncv will produce an error because the sample based (inverse) covariance matrix cannot be inverted in this situation.

For example, notice the error when running the following code:

library(GGMncv)

# p > n
set.seed(2)
main <- gen_net(p = 50, 
                edge_prob = 0.05)

set.seed(2)
y <- MASS::mvrnorm(n = 49, 
                   mu = rep(0, 50), 
                   Sigma = main$cors)

fit <- ggmncv(cor(y), n = nrow(y))
#> Error in solve.default(R): system is computationally singular: reciprocal condition number = 1.80172e-18

Solution

The solution is to provide an function for the initial matrix. To this end, GGMncv includes the function lediot_wolf which is a shrinkage estimator (Ledoit and Wolf 2004). It is important to note that any function can be used, so long as it return the inverse correlation matrix.

fit <- ggmncv(cor(y), n = nrow(y), 
              penalty = "atan",
              progress = FALSE,
              initial = ledoit_wolf, Y = y)

Notice the Y = y, which is used internally to pass additional arguments via ... to the function provided in initial.

The conditional dependence structure can then be plotted with

plot(get_graph(fit), 
     node_size = 5)

Here is an example of providing a function.

initial_ggmncv <- function(y, ...){
  Rinv <- corpcor::invcor.shrink(y, verbose = FALSE)
  return(Rinv)
}

fit2 <- ggmncv(cor(y), n = nrow(y), 
               penalty = "atan",
               progress = FALSE,
               initial = initial_ggmncv, 
               y = y)

plot(get_graph(fit2), 
     node_size = 5)

Perhaps it is of interest to compare performance, given that different initial values were used.

# ledoit and wolf
score_binary(estimate = fit$adj, 
             true = main$adj, 
             model_name = "lw")
#>       measure     score model_name
#> 1 specificity 0.9836770         lw
#> 2 sensitivity 0.1967213         lw
#> 3   precision 0.3870968         lw
#> 4      recall 0.1967213         lw
#> 5    f1_score 0.2608696         lw
#> 6         mcc 0.2498566         lw

# Shaffer and strimmer
score_binary(estimate = fit2$adj, 
             true = main$adj,  
             model_name = "ss")
#>       measure     score model_name
#> 1 specificity 0.9922680         ss
#> 2 sensitivity 0.1639344         ss
#> 3   precision 0.5263158         ss
#> 4      recall 0.1639344         ss
#> 5    f1_score 0.2500000         ss
#> 6         mcc 0.2749659         ss

Works but Bad Performance

Perhaps a trickier situation is when the covariance matrix can be inverted, but it is still ill-conditioned. This can occur when p approaches but does not exceed n. Here performance can be very bad.

# p -> n
main <- gen_net(p = 50, 
                edge_prob = 0.05)

y <- MASS::mvrnorm(n = 60, 
                   mu = rep(0, 50), 
                   Sigma = main$cors)

fit <- ggmncv(cor(y), n = nrow(y), 
              penalty = "atan",
              progress = FALSE)

score_binary(estimate = fit$adj, 
             true = main$adj)
#>       measure        score
#> 1 specificity  0.055841924
#> 2 sensitivity  0.934426230
#> 3   precision  0.049307958
#> 4      recall  0.934426230
#> 5    f1_score  0.093672966
#> 6         mcc -0.009181913

This is extremely problematic because there was no error, and the performance was terrible (note: 1 - specificity = the false positive rate).

Solution

One solution is again to provide a function to initial.

fit <- ggmncv(cor(y), n = nrow(y), 
              progress = FALSE,
              penalty = "atan",
              initial = ledoit_wolf, Y = y)

score_binary(estimate = fit$adj, 
             true = main$adj)
#>       measure     score
#> 1 specificity 0.9570447
#> 2 sensitivity 0.2131148
#> 3   precision 0.2063492
#> 4      recall 0.2131148
#> 5    f1_score 0.2096774
#> 6         mcc 0.1675807

An additional solution is to use L1-regularization, i.e.,

fit_l1 <- ggmncv(cor(y), n = nrow(y), 
              progress = FALSE, 
              penalty = "lasso")

score_binary(estimate = fit_l1$adj, 
             true = main$adj)
#>       measure      score
#> 1 specificity 0.99226804
#> 2 sensitivity 0.06557377
#> 3   precision 0.30769231
#> 4      recall 0.06557377
#> 5    f1_score 0.10810811
#> 6         mcc 0.12278924

A quick comparison of KL-divergence

# atan
kl_mvn(true = solve(main$cors), estimate = fit$Theta)
#> [1] 1.918541

# lasso
kl_mvn(true = solve(main$cors), estimate = fit_l1$Theta)
#> [1] 2.522567

References

Ledoit, Olivier, and Michael Wolf. 2004. “A Well-Conditioned Estimator for Large-Dimensional Covariance Matrices.” Journal of Multivariate Analysis 88 (2): 365–411.
Williams, Donald R, and Philippe Rast. 2020. “Back to the Basics: Rethinking Partial Correlation Network Methodology.” British Journal of Mathematical and Statistical Psychology 73 (2): 187–212.
Williams, Donald R, Mijke Rhemtulla, Anna C Wysocki, and Philippe Rast. 2019. “On Nonregularized Estimation of Psychological Networks.” Multivariate Behavioral Research 54 (5): 719–50.