Instantly share code, notes, and snippets.

Embed
What would you like to do?
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment