WilliamsRast.R

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)