Skip to content
{{ message }}

Instantly share code, notes, and snippets.

# SachaEpskamp/WilliamsRast.R

Created Apr 7, 2018
 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)){ 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)
to join this conversation on GitHub. Already have an account? Sign in to comment