Robocov application on 1-band precision matrix

Banded precision matrix for population.

banded_prec_sim = function(N, P){
  diags <- list()
  diags[[1]] <- rep(1, P)
  diags[[2]] <- rep(-0.5, P)
  Kinv <- bandSparse(P, k = -(0:1), diag = diags, symm = TRUE)
  K <- solve(Kinv)
  corSigma <- cov2cor(K)
  data <- MASS::mvrnorm(N,rep(0,P),corSigma)
  ll = list("dat" = data, "cor" = corSigma)
  return(ll)
}

nloglik = function(data, cormat){
  llik = 0
  for(m in 1:nrow(data)){
    idx = which(!is.na(data[m,]))
    if(length(idx) > 2){
      llik = llik + emdbook::dmvnorm(data[m, idx], rep(0, length(idx)), cormat[idx, idx], log = T)
    }
  }
  return(-llik)
}

angle_norm = function(S, Sigma){
  dist = 1 - (tr(as.matrix(cov2cor(S)%*%cov2cor(Sigma))))/(norm(cov2cor(S), type = "F")* norm(Sigma, type = "F"))
  return(dist)
}

Data generation

N=500
P=50
prop_missing = 0.25

corSigma = as.matrix(banded_prec_sim(N, P)$cor)

data = banded_prec_sim(N, P)$dat

  #######################   Turn some of the entries to NA   ###################################

data_missing = t(apply(data, 1, function(x){
  if(prop_missing > 0){
    rand = sample(1:length(x), floor(prop_missing*length(x)), replace = F)
    y = x
    y[rand] = NA
    return(y)
  }else{
    return(x)
  }}))

Pop

corrplot(corSigma,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

p-Pop

pcorSigma = -cov2cor(as.matrix(solve(corSigma)))
diag(pcorSigma) = 1
corrplot(pcorSigma,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.cex = 1.5, tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

Standard

standard_cor = cor(data_missing, use = "pairwise.complete.obs")
corrplot(standard_cor,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.cex = 1.5, tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

Robocov

robocov_box_cor = Robocov_cor(data_with_missing = data_missing, loss = "lasso")
corrplot(robocov_box_cor,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.cex = 1.5, tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

p-Robocov

robo_prec = Robocov_precision(data_with_missing = data_missing, alpha = 0.1, lambda=1)
corrplot(robo_prec,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.cex = 1.5, tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

CorShrink

cov_sample_ML <-  CorShrinkData(data_missing, sd_boot = FALSE,
                                  ash.control = list())
corshrink_cor = cov2cor(cov_sample_ML$cor)
corrplot(corshrink_cor,  diag = TRUE,
         col = colorRampPalette(c("lightblue4", "lightblue2", "white", "indianred1", "indianred3"))(200),
         tl.pos = "n", tl.cex = 1.5, tl.col = "black",
         rect.col = "white",na.label.col = "white",
         method = "color", type = "lower", tl.srt=45)

Eigenvalue trends

NUM=20
df = cbind(eigen(corSigma)$values[1:NUM], eigen(standard_cor)$values[1:NUM],
eigen(robocov_box_cor)$values[1:NUM], eigen(corshrink_cor)$values[1:NUM])
colnames(df) = c("Pop-corr", "Sample-Corr", "Robocov-Corr", "CorShrink")
df

        Pop-corr Sample-Corr Robocov-Corr  CorShrink
 [1,] 25.5000000  24.9004541   19.5078745 24.5178502
 [2,]  8.5000000   8.9987425    9.8740390  8.8991064
 [3,]  4.2500000   4.2498613    4.9544810  4.2075595
 [4,]  2.5500000   2.4760684    2.9898691  2.4658236
 [5,]  1.7000000   1.8546298    2.1500979  1.8463956
 [6,]  1.2142857   1.3430415    1.5686767  1.3414157
 [7,]  0.9107143   1.0061640    1.3243405  1.0059649
 [8,]  0.7083333   0.6876399    0.9554003  0.6852039
 [9,]  0.5666667   0.5825276    0.7351148  0.5837400
[10,]  0.4636364   0.4943548    0.6139360  0.4954514
[11,]  0.3863636   0.4349342    0.5333996  0.4343129
[12,]  0.3269231   0.3541334    0.4553435  0.3543637
[13,]  0.2802198   0.3238786    0.4315094  0.3232409
[14,]  0.2428571   0.3156612    0.4086670  0.3153196
[15,]  0.2125000   0.2681948    0.3647196  0.2694940
[16,]  0.1875000   0.2580605    0.3298700  0.2592631
[17,]  0.1666667   0.2388775    0.2771693  0.2376500
[18,]  0.1491228   0.2016679    0.2645459  0.2020885
[19,]  0.1342105   0.1971960    0.2271823  0.1974092
[20,]  0.1214286   0.1619391    0.2148577  0.1625495

eigs.df <- data.frame ("x" = rep(1:NUM, 4),
                       "y" = c(df[,1], df[,2], df[,3], df[,4]),
                       "color" = factor(c(rep("Pop", NUM), 
                                   rep("Sample",  NUM), 
                                   rep("Robocov", NUM),
                                   rep("Corshrink", NUM)),
                                   levels = c("Pop", "Robocov", "Corshrink", "Sample")))

p <- ggplot(eigs.df, aes(x=x, y=y, colour=color, linetype = color)) + geom_line(lty = 1, lwd = 0.7) +
  scale_linetype_manual(values = c(rep("solid", 4), rep("dashed", 1))) +
  scale_colour_manual(values=c("#000000", "blue", "green", "gold",
                               "red", "#0072B2", "#CC79A7", "#F0E442")) +
  ggtitle(paste0("n=", N, ", p=", P, " \u03C0", "=", prop_missing)) + xlab("Index") + ylab("eigenvalues")+
  theme_bw()
p

NUM=20
df = cbind(eigen(pcorSigma)$values[1:NUM], eigen(robo_prec)$values[1:NUM])
colnames(df) = c("Pop", "Robocov")
df

           Pop  Robocov
 [1,] 1.998103 1.977503
 [2,] 1.992421 1.942069
 [3,] 1.982973 1.897682
 [4,] 1.969797 1.840530
 [5,] 1.952942 1.781324
 [6,] 1.932472 1.724520
 [7,] 1.908465 1.674409
 [8,] 1.881012 1.606774
 [9,] 1.850217 1.537470
[10,] 1.816197 1.461318
[11,] 1.779081 1.415708
[12,] 1.739009 1.364062
[13,] 1.696134 1.290637
[14,] 1.650618 1.260008
[15,] 1.602635 1.200976
[16,] 1.552365 1.149770
[17,] 1.500000 1.104155
[18,] 1.445738 1.091344
[19,] 1.389786 1.066381
[20,] 1.332355 1.022802

eigs.df <- data.frame ("x" = rep(1:NUM, 4),
                       "y" = c(df[,1], df[,2]),
                       "color" = factor(c(rep("Pop", NUM),  
                                   rep("Robocov", NUM)),
                                   levels = c("Pop", "Robocov")))



p <- ggplot(eigs.df, aes(x=x, y=y, colour=color, linetype = color)) + geom_line(lty = 1, lwd = 0.7) +
  scale_linetype_manual(values = c(rep("solid", 4), rep("dashed", 1))) +
  scale_colour_manual(values=c("#000000", "blue", "green", "gold",
                               "red", "#0072B2", "#CC79A7", "#F0E442")) +
  ggtitle(paste0("n=", N, ", p=", P, " \u03C0", "=", prop_missing)) + xlab("Index") + ylab("eigenvalues")+
  theme_bw()
p

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] corrplot_0.84   ggplot2_3.1.1   CorShrink_0.1-6 Robocov_0.1-6  
[5] CVXR_0.99-2     psych_1.8.12    Matrix_1.2-14   corpcor_1.6.9  
[9] glasso_1.10    

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5  purrr_0.3.2       ashr_2.2-38      
 [4] reshape2_1.4.3    lattice_0.20-35   colorspace_1.4-1 
 [7] htmltools_0.3.6   yaml_2.2.0        gmp_0.5-13.2     
[10] rlang_0.4.2       pillar_1.3.1      R.oo_1.22.0      
[13] mixsqp_0.1-97     withr_2.1.2       glue_1.3.1       
[16] foreign_0.8-70    Rmpfr_0.7-1       R.utils_2.7.0    
[19] bit64_0.9-7       scs_1.1-1         foreach_1.4.4    
[22] plyr_1.8.4        stringr_1.4.0     munsell_0.5.0    
[25] gtable_0.3.0      workflowr_1.1.1   R.methodsS3_1.7.1
[28] codetools_0.2-15  evaluate_0.12     labeling_0.3     
[31] knitr_1.20        doParallel_1.0.14 pscl_1.5.2       
[34] parallel_3.5.1    Rcpp_1.0.1        backports_1.1.4  
[37] scales_1.0.0      truncnorm_1.0-8   bit_1.1-14       
[40] gridExtra_2.3     mnormt_1.5-5      digest_0.6.19    
[43] stringi_1.4.3     dplyr_0.8.0.1     grid_3.5.1       
[46] rprojroot_1.3-2   ECOSolveR_0.4     tools_3.5.1      
[49] magrittr_1.5      tibble_2.1.1      lazyeval_0.2.2   
[52] glmnet_2.0-18     pkgconfig_2.0.2   crayon_1.3.4     
[55] whisker_0.3-2     MASS_7.3-50       SQUAREM_2017.10-1
[58] assertthat_0.2.1  rmarkdown_1.10    iterators_1.0.10 
[61] R6_2.4.0          nlme_3.1-137      git2r_0.23.0     
[64] compiler_3.5.1