genGGMexample.R

library("BDgraph")
library("glasso")
library("bootnet") # Install from github: devtools::install_github("sachaepskamp/bootnet)

# Seed to reproduce results:
set.seed(1)

# Simulate 20% sparse network:
bnres <- genGGM(20, 0.8, graph = "random")

# Sparsify inverse:
Kappa <- diag(20) - bnres

# True covariance matrix:
Sigma <- solve(Kappa)

# Estimate the glassopath:
lambdaseq <- exp(seq(log(0.0001), log(1), length = 10000))
path <- glassopath(Sigma, penalize.diagonal = FALSE, rholist = lambdaseq)

# Compare to true networks:
comparison <- list()
for (i in 1:dim(path$wi)[3]){
  comparison[[i]] <- data.frame(
    lambda = lambdaseq[i],
    falsePositives = sum(Kappa[upper.tri(Kappa)] == 0 & path$wi[,,i][upper.tri( path$wi[,,i])] != 0)/sum(Kappa[upper.tri(Kappa)] == 0),
    falseNegatives = sum(Kappa[upper.tri(Kappa)] != 0 & path$wi[,,i][upper.tri( path$wi[,,i])] == 0)/sum(Kappa[upper.tri(Kappa)] != 0)
  )
}

# Combine:
Results <- as.data.frame(do.call(rbind,comparison))

# Plot:
options(scipen=5)
plot(Results$lambda, Results$falsePositives, type = "l", log = "x",
     ylab = "", xlab = "GLASSO tuning parameter", bty = "n", col = "red", ylim=c(0,1))
points(Results$lambda, Results$falseNegatives, type = "l", col = "blue")
legend("topleft",legend=c("False positive rate","False negative rate"),col=c("red", "blue"),
       lty = 1, bty = "n", cex = 1)