library(glasso)
library(corpcor)
library(CorShrink)
library(network)
source("isee_all.R")
library(scales)
library(statnet)

In this script, we consider a simple network, and conisder the data from this network structure and try to estimate from the data, the underlying network using three approaches to graph inference - ISEE, pCorShrink and GLASSO. This example has been taken from this lecture note.

 n <- 100
 p <- 10

 # Construct a sparse precision matrix

 trueP <- diag(p)
 for(j in 5:p){
  trueP[j,j-1] <- -.4
  trueP[j-1,j] <- -.4
 }
 trueS <- solve(trueP) ## sample covariance matrix 
 trueA <- ifelse(trueP!=0 & row(trueP)!=col(trueP),1,0) ## graph structure

 trueP
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1    0    0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
##  [2,]    0    1    0  0.0  0.0  0.0  0.0  0.0  0.0   0.0
##  [3,]    0    0    1  0.0  0.0  0.0  0.0  0.0  0.0   0.0
##  [4,]    0    0    0  1.0 -0.4  0.0  0.0  0.0  0.0   0.0
##  [5,]    0    0    0 -0.4  1.0 -0.4  0.0  0.0  0.0   0.0
##  [6,]    0    0    0  0.0 -0.4  1.0 -0.4  0.0  0.0   0.0
##  [7,]    0    0    0  0.0  0.0 -0.4  1.0 -0.4  0.0   0.0
##  [8,]    0    0    0  0.0  0.0  0.0 -0.4  1.0 -0.4   0.0
##  [9,]    0    0    0  0.0  0.0  0.0  0.0 -0.4  1.0  -0.4
## [10,]    0    0    0  0.0  0.0  0.0  0.0  0.0 -0.4   1.0
 library(network)
 g <- network(trueA)
 plot(g,label=1:p,main="True network")

## Generate Data

x <- matrix(rnorm(n*p),n,p)%*%chol(trueS)
S <- cov(x)
round(S,2)
##        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
##  [1,]  0.97  0.05  0.07 -0.11 -0.07 -0.01 -0.04 -0.06  0.14  0.15
##  [2,]  0.05  1.02 -0.16  0.00  0.05 -0.05 -0.08  0.02  0.01  0.06
##  [3,]  0.07 -0.16  1.05  0.09 -0.04  0.24  0.03  0.22  0.11 -0.11
##  [4,] -0.11  0.00  0.09  1.06  0.54  0.15  0.02 -0.11 -0.17 -0.03
##  [5,] -0.07  0.05 -0.04  0.54  1.70  0.64  0.32  0.18  0.05  0.04
##  [6,] -0.01 -0.05  0.24  0.15  0.64  1.60  0.74  0.58  0.16 -0.04
##  [7,] -0.04 -0.08  0.03  0.02  0.32  0.74  1.55  0.92  0.54  0.28
##  [8,] -0.06  0.02  0.22 -0.11  0.18  0.58  0.92  1.62  0.92  0.16
##  [9,]  0.14  0.01  0.11 -0.17  0.05  0.16  0.54  0.92  1.75  0.63
## [10,]  0.15  0.06 -0.11 -0.03  0.04 -0.04  0.28  0.16  0.63  1.28
nr  <- 100
 rho <- seq(0.1,1,length=nr)
 bic <- rho
 for(j in 1:nr){
  a       <- glasso(S,rho[j])
  p_off_d <- sum(a$wi!=0 & col(S)<row(S))
  bic[j]  <- -2*(a$loglik) + p_off_d*log(n)
 }

 best <- which.min(bic)
 plot(rho,bic)
 points(rho[best],bic[best],pch=19)

glasso network

a <- glasso(S,rho[best])
names(a)
## [1] "w"       "wi"      "loglik"  "errflag" "approx"  "del"     "niter"
P <- a$wi
A <- ifelse(P!=0 & row(P)!=col(P),1,0)

 g <- network(A)
 plot(g,label=1:p,main="Estimated network")

ISEE

regfactor = "log"  # or "one", "sqrt"
npermu = 1         # or >= 2
sis.use = 0        # or 1, whether to use SIS for screening
bia.cor = 0        # or 1, whether to apply bias correction for ISEE
obj.n = isee(x, regfactor, npermu, sis.use, bia.cor = 0) 
pcor2 <- -cov2cor(obj.n$Omega.isee)
diag(pcor2) <- 1

A <- ifelse(pcor2!=0 & row(pcor2)!=col(pcor2),1,0)

g <- network(A)
plot(g,label=1:p,main="Estimated network by ISEE")

pCorShrink

pcor4 <- pCorShrinkData(x, reg_type = "lm",
                        maxiter = 1000,
                        ash.control = list(mixcompdist = "halfuniform"))
A <- ifelse(pcor4 > 0.1 & row(pcor4)!=col(pcor4),1,0)

g <- network(A)
plot(g,label=1:p,main="Estimated network by pCorShrink+LM")

IGraph plot for pCorShrink

g = network(5*pcor4, matrix.type = "adjacency", directed = FALSE,
               names.eval="weight", ignore.eval=FALSE)
#get.edge.value(g, "weight")
#plot(g,label=1:p,main="Estimated network by pCorShrink+LM")
col2 <- c("blue", "white", "red")
gplot(g,gmode="graph",label = 1:p, edge.lwd=g%e%'weight', 
      edge.col = colorRampPalette(col2)(200)[rescale(g%e%'weight', 
                                          to = c(1, 200), from = c(-5,5))],
      vertex.col="black")