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