Create a gist now

Instantly share code, notes, and snippets.

Embed
library("BDgraph")
library("glasso")
# Seed to reproduce results:
set.seed(1)
# Simulate 20% sparse network and data:
bdres <- bdgraph.sim(20, n = 10000, prob = 0.8, b = 20, D = diag(20, 20))
# Sparsify inverse:
Kappa <- bdres$K * bdres$G + diag(diag(bdres$K))
# 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(bdres$G[upper.tri(bdres$G)] == 0 & path$wi[,,i][upper.tri( path$wi[,,i])] != 0)/sum(bdres$G[upper.tri(bdres$G)] == 0),
falseNegatives = sum(bdres$G[upper.tri(bdres$G)] != 0 & path$wi[,,i][upper.tri( path$wi[,,i])] == 0)/sum(bdres$G[upper.tri(bdres$G)] != 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