In this script, we compare a popular sparse inverse correlation estimation method called CLIME and test it against ISEE and the pCorShrink approaches.

Packages

require(graphics)
library(igraph)
library(MASS) 
library(Matrix)
library(corpcor)
library(corrplot)
source("isee_all.R") ## ISEE codes
library(CorShrink)
library(scio)

Banded precision matrix

band.mat <- function(a, p, K=1, permu=c(1:p)){
  ones = rep(1,p)
  Omega0 = a*ones%*%t(ones)
  diag(Omega0) = rep(1,p)
  Omega = 1*band(Omega0,-K,K)
  Sigma = qr.solve(Omega)
  Sigma = Sigma*(abs(Sigma)>1e-4)
  Sigma.half=chol(Sigma)
  Sigma.half = Sigma.half*(abs(Sigma.half)>1e-4)
  Sigma = Sigma[permu,permu]
  Omega = Omega[permu,permu]
  Sigma.half = Sigma.half[permu,permu]
  obj = list(Sigma=Sigma, Omega = Omega, Sigma.half = Sigma.half)
}

make.data <- function(Sigma.half, n, p, seed){
  set.seed(seed)  
  
  X = matrix(rnorm(n*p),n,p)%*%Sigma.half
  return(X)
}

n=50, p=100

n = 50
p = 100
obj = band.mat(a=0.5, p, K = 1)
Sig.half = obj$Sigma.half
Ome.true = obj$Omega
X.mat = make.data(Sig.half, n, p, seed = 1000)

Original

pcorSigma <- -as.matrix(cov2cor(obj$Omega))
diag(pcorSigma) <- rep(1, dim(pcorSigma)[1])
col2 <- c("blue", "white", "red")
corrplot(pcorSigma, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

CLIME

wmat <- scio.cv(X.mat)$w
clime_pcor <- -cov2cor(wmat)
diag(clime_pcor) <- 1
col2 <- c("blue", "white", "red")
corrplot(clime_pcor, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

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.mat, regfactor, npermu, sis.use, bia.cor = 0) 
pseemat <- -cov2cor(obj.n$Omega.isee)
diag(pseemat) <- rep(1, dim(pseemat)[1])
corrplot(pseemat, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

pCorShrink

pcor2 <- pCorShrinkData(X.mat, reg_type = "lm")
## Number of samples less than number of features, switching to glmnet regression
## using the default alpha = 1 for glmnet: same as Lasso. Change glmnet_alpha between 0 and 1
##             for further options, alpha = 0 -> ridge, alpha = 0.5 -> elastic net
corrplot(pcor2, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

n=100, p=100

n = 100
p = 100
obj = band.mat(a=0.5, p, K = 1)
Sig.half = obj$Sigma.half
Ome.true = obj$Omega
X.mat = make.data(Sig.half, n, p, seed = 1000)

Original

pcorSigma <- -as.matrix(cov2cor(obj$Omega))
diag(pcorSigma) <- rep(1, dim(pcorSigma)[1])
col2 <- c("blue", "white", "red")
corrplot(pcorSigma, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

CLIME

wmat <- scio.cv(X.mat)$w
clime_pcor <- -cov2cor(wmat)
diag(clime_pcor) <- 1
col2 <- c("blue", "white", "red")
corrplot(clime_pcor, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

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.mat, regfactor, npermu, sis.use, bia.cor = 0) 
pseemat <- -cov2cor(obj.n$Omega.isee)
diag(pseemat) <- rep(1, dim(pseemat)[1])
corrplot(pseemat, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

pCorShrink

pcor2 <- pCorShrinkData(X.mat, reg_type = "lm")
corrplot(pcor2, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

n=1000, p=100

n = 1000
p = 100
obj = band.mat(a=0.5, p, K = 1)
Sig.half = obj$Sigma.half
Ome.true = obj$Omega
X.mat = make.data(Sig.half, n, p, seed = 1000)

Original

pcorSigma <- -as.matrix(cov2cor(obj$Omega))
diag(pcorSigma) <- rep(1, dim(pcorSigma)[1])
col2 <- c("blue", "white", "red")
corrplot(pcorSigma, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

CLIME

wmat <- scio.cv(X.mat)$w
clime_pcor <- -cov2cor(wmat)
diag(clime_pcor) <- 1
col2 <- c("blue", "white", "red")
corrplot(clime_pcor, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

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.mat, regfactor, npermu, sis.use, bia.cor = 0) 
pseemat <- -cov2cor(obj.n$Omega.isee)
diag(pseemat) <- rep(1, dim(pseemat)[1])
corrplot(pseemat, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")

pCorShrink

pcor2 <- pCorShrinkData(X.mat, reg_type = "lm")
corrplot(pcor2, diag = FALSE, col = colorRampPalette(col2)(200), tl.pos = "td", 
         tl.col = "black", tl.cex = 0.8, rect.col = "white", 
         na.label.col = "white", method = "color", type = "upper")