Packages

library(MASS) 
library(Matrix)
library(corrplot)
library(CorShrink)
library(corpcor)
library(scales)
library(gridExtra)
library(ggplot2)

We look at the Frobenius distance between the estimated partial correlation matrix using ISEE and pCorSrhink agaisnt standard methods.

Hub matrix

Original partial correlation matrix

P <- 100
block <- 10
mat <- 0.3*diag(1,block) + 0.7*rep(1,block) %*% t(rep(1, block))
Sigma <-   bdiag(mat, mat, mat, mat, mat, mat, mat, mat, mat, mat)
corSigma <- cov2cor(Sigma)
pcorSigma <- as.matrix(cor2pcor(corSigma))
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")

frobs_list <- get(load("R/Hub/frob_dist_hub.rda"))
n_vec = c(30, 50, 100, 110, 500, 1000)
M = 30

p_list <- list()
for(n in 1:length(n_vec)){
  df <- data.frame("method" = factor(c(rep("empirical", M), rep("isee", M), 
                                     rep("isee_X", M), 
                                     rep("pcorshrink_lasso", M)), 
                              levels = c("empirical", "isee",
                                        "isee_X",
                                        "pcorshrink_lasso")),
                 "distance" = log(c(frobs_list$empirical[n,],
                                    frobs_list$isee[n,],
                                    frobs_list$isee_X[n,],
                                    frobs_list$pcorshrink_lasso[n,])))
p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p_list[[n]] <- p + geom_boxplot() + theme_bw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", n_vec[n], " p = ", 100)) +scale_y_continuous(breaks= pretty_breaks())
}
grid.arrange(p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]],
             p_list[[5]],p_list[[6]], nrow = 3,  ncol=2, as.table=TRUE)

Toeplitz matrix

Original partial correlation matrix

DM_toeplitz = function(n,P){
  library("MASS")
  index1=sort(sample(seq(1:n),(n/2)))
  index2=seq(1:n)[-index1]
  Sigmatp=function(P){
    a=array(0,dim=c(P,P))
    for(i in 1:P){
      for(j in 1:P){
        a[i,j]=max(1-0.1*(abs(i-j)),0)
      }
    }
    return(a)
  }
  Sigma = Sigmatp(P)
  data = mvrnorm(n,rep(0,P),Sigma)
  Xtest = data[index2,]
  Xtrain = data[index1,]
  Omega = solve(Sigma)
  return(list(Xtrain = Xtrain, Xtest = Xtest, Sigma = Sigma))
}
ll <- DM_toeplitz(n=10, P=100)
Sigma <- ll$Sigma
corSigma <- cov2cor(Sigma)
pcorSigma <- as.matrix(cor2pcor(corSigma))
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")

frobs_list <- get(load("R/Toeplitz/frob_dist_toeplitz.rda"))
n_vec = c(30, 50, 100, 110, 500, 1000)
M = 30

p_list <- list()
for(n in 1:length(n_vec)){
  df <- data.frame("method" = factor(c(rep("empirical", M), rep("isee", M), 
                                     rep("isee_X", M), 
                                     rep("pcorshrink_lasso", M)), 
                              levels = c("empirical", "isee",
                                        "isee_X",
                                        "pcorshrink_lasso")),
                 "distance" = log(c(frobs_list$empirical[n,],
                                    frobs_list$isee[n,],
                                    frobs_list$isee_X[n,],
                                    frobs_list$pcorshrink_lasso[n,])))
p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p_list[[n]] <- p + geom_boxplot() + theme_bw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", n_vec[n], " p = ", 100)) +scale_y_continuous(breaks= pretty_breaks())
}
grid.arrange(p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]],
             p_list[[5]], p_list[[6]], nrow = 3,  ncol=2, as.table=TRUE)

Banded precision 1

Original partial correlation matrix

P <- 100
diags <- list()
diags[[1]] <- rep(1, 100)
diags[[2]] <- rep(-0.5, 100)
Kinv <- bandSparse(100, k = -(0:1), diag = diags, symm = TRUE)
K <- solve(Kinv)
corSigma <- cov2cor(K)
pcorSigma <- as.matrix(cor2pcor(corSigma))
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")

frobs_list <- get(load("R/Banded_Precision_1/frob_dist_banded_precision.rda"))
n_vec = c(30, 50, 100, 110, 500, 1000)
M = 30

p_list <- list()
for(n in 1:length(n_vec)){
  df <- data.frame("method" = factor(c(rep("empirical", M), rep("isee", M), 
                                     rep("isee_X", M), 
                                     rep("pcorshrink_lasso", M)), 
                              levels = c("empirical", "isee",
                                        "isee_X",
                                        "pcorshrink_lasso")),
                 "distance" = log(c(frobs_list$empirical[n,],
                                    frobs_list$isee[n,],
                                    frobs_list$isee_X[n,],
                                    frobs_list$pcorshrink_lasso[n,])))
p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p_list[[n]] <- p + geom_boxplot() + theme_bw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", n_vec[n], " p = ", 100)) +scale_y_continuous(breaks= pretty_breaks())
}
grid.arrange(p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]],
             p_list[[5]], p_list[[6]], nrow = 3,  ncol=2, as.table=TRUE)

Banded precision 2

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)
}
P <- 100
obj = band.mat(a=0.5, P, K = 1)
Sig.true = obj$Sigma
Ome.true = obj$Omega
corSigma = cov2cor(Sig.true)
pcorSigma = as.matrix(cor2pcor(corSigma))
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")

frobs_list <- get(load("R/Banded_Precision_2/frob_dist_banded_precision_2.rda"))
n_vec = c(30, 50, 100, 110, 500, 1000)
M = 30

p_list <- list()
for(n in 1:length(n_vec)){
  df <- data.frame("method" = factor(c(rep("empirical", M), rep("isee", M), 
                                     rep("isee_X", M), 
                                     rep("pcorshrink_lasso", M)), 
                              levels = c("empirical", "isee",
                                        "isee_X",
                                        "pcorshrink_lasso")),
                 "distance" = log(c(frobs_list$empirical[n,],
                                    frobs_list$isee[n,],
                                    frobs_list$isee_X[n,],
                                    frobs_list$pcorshrink_lasso[n,])))
p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p_list[[n]] <- p + geom_boxplot() + theme_bw() +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", n_vec[n], " p = ", 100)) +scale_y_continuous(breaks= pretty_breaks())
}
grid.arrange(p_list[[1]], p_list[[2]], p_list[[3]], p_list[[4]],
             p_list[[5]], p_list[[6]], nrow = 3,  ncol=2, as.table=TRUE)