CorShrink on hub matrices : part 2

library(MASS)
library(Matrix)
library(glasso)
library(CorShrink)
library(corpcor)

set up

P <- 200
block <- 20
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)

Sparse correlation matrix

col=c(rep(rgb(0,1,0), 10000),rev(rgb(seq(1,0,length=1000),1,seq(1,0,length=1000))),
      rgb(1,seq(1,0,length=1000),seq(1,0,length=1000)), rep(rgb(1,0,0),10000))

image(as.matrix(cov2cor(Sigma)),
      col=col, main=paste0("pop corr matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

Sparse precision matrix

col=c(rep(rgb(0,1,0), 10000),rev(rgb(seq(1,0,length=1000),1,seq(1,0,length=1000))),
      rgb(1,seq(1,0,length=1000),seq(1,0,length=1000)), rep(rgb(1,0,0),10000))

image(as.matrix(solve(cov2cor(Sigma))),
      col=col, main=paste0("precision matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

## n = 10, p = 200

set.seed(100)
P <- dim(Sigma)[1]
n=10
data <- MASS::mvrnorm(n,rep(0,P),Sigma)

nr  <- 100
rho <- seq(0.1,10,length=nr)
bic <- rho
S <- cov(data)
for(j in 1:nr){
  a       <- glasso::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)

a <- glasso::glasso(S,rho[best])

## Strimmer Shafer 

strimmer_sample <- corpcor::cor.shrink(data)

## Estimating optimal shrinkage intensity lambda (correlation matrix): 0.6896

## CorShrink

cov_sample_ML <-  CorShrinkData(data, sd_boot = FALSE, optmethod = "mixEM",
                     image.control = list(x.cex = 0.3, y.cex = 0.3),
                     ash.control = list())

## Due to absence of package REBayes, switching to EM algorithm

## ash cor only and ash cor PD matrices are different

col=c(rep(rgb(0,1,0), 1),rev(rgb(seq(1,0,length=10000),1,seq(1,0,length=10000))),
      rgb(1,seq(1,0,length=10000),seq(1,0,length=10000)), rep(rgb(1,0,0),1))

par(mfrow = c(2,2))

image(as.matrix(cov2cor(S)),
      col=col, main=paste0("sample corr matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov2cor(a$w)),
      col=col, main=paste0("GLASSO"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(strimmer_sample,
      col=col, main=paste0("shafer strimmer"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov_sample_ML$ash_cor_PD),
      col=col, main=paste0("CorShrink"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

n = 25, p = 200

set.seed(100)
P <- dim(Sigma)[1]
n=25
data <- MASS::mvrnorm(n,rep(0,P),Sigma)

nr  <- 100
rho <- seq(0.1,10,length=nr)
bic <- rho
S <- cov(data)
for(j in 1:nr){
  a       <- glasso::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)

a <- glasso::glasso(S,rho[best])

## Strimmer Shafer 

strimmer_sample <- corpcor::cor.shrink(data)

## Estimating optimal shrinkage intensity lambda (correlation matrix): 0.4923

## CorShrink

cov_sample_ML <-  CorShrinkData(data, sd_boot = FALSE, optmethod = "mixEM",
                     image.control = list(x.cex = 0.3, y.cex = 0.3),
                     ash.control = list())

## Due to absence of package REBayes, switching to EM algorithm

## ash cor only and ash cor PD matrices are different

col=c(rep(rgb(0,1,0), 1),rev(rgb(seq(1,0,length=10000),1,seq(1,0,length=10000))),
      rgb(1,seq(1,0,length=10000),seq(1,0,length=10000)), rep(rgb(1,0,0),1))

par(mfrow = c(2,2))

image(as.matrix(cov2cor(S)),
      col=col, main=paste0("sample corr matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov2cor(a$w)),
      col=col, main=paste0("GLASSO"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(strimmer_sample,
      col=col, main=paste0("shafer strimmer"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov_sample_ML$ash_cor_PD),
      col=col, main=paste0("CorShrink"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

n = 50, p = 200

set.seed(100)
P <- dim(Sigma)[1]
n=50
data <- MASS::mvrnorm(n,rep(0,P),Sigma)

nr  <- 100
rho <- seq(0.1,10,length=nr)
bic <- rho
S <- cov(data)
for(j in 1:nr){
  a       <- glasso::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)

a <- glasso::glasso(S,rho[best])

## Strimmer Shafer 

strimmer_sample <- corpcor::cor.shrink(data)

## Estimating optimal shrinkage intensity lambda (correlation matrix): 0.3211

## CorShrink

cov_sample_ML <-  CorShrinkData(data, sd_boot = FALSE, optmethod = "mixEM",
                     image.control = list(x.cex = 0.3, y.cex = 0.3),
                     ash.control = list())

## Due to absence of package REBayes, switching to EM algorithm

## ash cor only and ash cor PD matrices are different

col=c(rep(rgb(0,1,0), 1),rev(rgb(seq(1,0,length=10000),1,seq(1,0,length=10000))),
      rgb(1,seq(1,0,length=10000),seq(1,0,length=10000)), rep(rgb(1,0,0),1))

par(mfrow = c(2,2))

image(as.matrix(cov2cor(S)),
      col=col, main=paste0("sample corr matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov2cor(a$w)),
      col=col, main=paste0("GLASSO"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(strimmer_sample,
      col=col, main=paste0("shafer strimmer"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov_sample_ML$ash_cor_PD),
      col=col, main=paste0("CorShrink"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

n = 75, p = 200

set.seed(100)
P <- dim(Sigma)[1]
n=75
data <- MASS::mvrnorm(n,rep(0,P),Sigma)

nr  <- 100
rho <- seq(0.1,10,length=nr)
bic <- rho
S <- cov(data)
for(j in 1:nr){
  a       <- glasso::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)

a <- glasso::glasso(S,rho[best])

## Strimmer Shafer 

strimmer_sample <- corpcor::cor.shrink(data)

## Estimating optimal shrinkage intensity lambda (correlation matrix): 0.2222

## CorShrink

cov_sample_ML <-  CorShrinkData(data, sd_boot = FALSE, optmethod = "mixEM",
                     image.control = list(x.cex = 0.3, y.cex = 0.3),
                     ash.control = list())

## Due to absence of package REBayes, switching to EM algorithm

## ash cor only and ash cor PD matrices are different

col=c(rep(rgb(0,1,0), 1),rev(rgb(seq(1,0,length=10000),1,seq(1,0,length=10000))),
      rgb(1,seq(1,0,length=10000),seq(1,0,length=10000)), rep(rgb(1,0,0),1))

par(mfrow = c(2,2))

image(as.matrix(cov2cor(S)),
      col=col, main=paste0("sample corr matrix"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov2cor(a$w)),
      col=col, main=paste0("GLASSO"), cex.main=1,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(strimmer_sample,
      col=col, main=paste0("shafer strimmer"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))

image(as.matrix(cov_sample_ML$ash_cor_PD),
      col=col, main=paste0("CorShrink"), cex.main=2,
      xaxt = "n", yaxt = "n", zlim=c(-1,1))