Simulation studies partial correlation

Kushal K Dey

5/22/2018

We perform box plot representation of the Frobenius distance between the population correlation matrix and correlation matrix estimators - CorShrink, Shafer Strimmer and the GLASSO methods.

Packages

library(gridExtra)
library(ggplot2)
library(scales)

Hub correlation

We simulate data from a hub structured correlation matrix under different settings of \(n\) and \(p\).

First run the code in

source("../code/Figure5/corshrink_hub_pcor.R")  ## n=50, p=100 

## last line : save(frob_vals, file = paste0("hub_cmd_boot_n_", n, "_P_", P, "_results_pcor.rda"))
## change the file to the destination and file name you prefer.

Change the \(n\) and \(p\) accordingly, we pick \(p=100\) and vary \(n\) from \(30, 50, 100, 1000\).

The final output from these files were saved under shared_output/Pcor/ folders and can be read as follows.

hub_n_10_P_100 <- get(load("../shared_output/Pcor/hub_cmd_boot_n_30_P_100_results_pcor.rda"))
hub_n_50_P_100 <- get(load("../shared_output/Pcor/hub_cmd_boot_n_50_P_100_results_pcor.rda"))
hub_n_100_P_100 <- get(load("../shared_output/Pcor/hub_cmd_boot_n_100_P_100_results_pcor.rda"))
hub_n_1000_P_100 <- get(load("../shared_output/Pcor/hub_cmd_boot_n_1000_P_100_results_pcor.rda"))

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(hub_n_10_P_100[,1], hub_n_10_P_100[,2],
                                    hub_n_10_P_100[,3], hub_n_10_P_100[,4],
                                    hub_n_10_P_100[,5], hub_n_10_P_100[,6],
                                    hub_n_10_P_100[,7], hub_n_10_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p1 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 30, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(hub_n_50_P_100[,1], hub_n_50_P_100[,2],
                                    hub_n_50_P_100[,3], hub_n_50_P_100[,4],
                                    hub_n_50_P_100[,5], hub_n_50_P_100[,6],
                                    hub_n_50_P_100[,7], hub_n_50_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p2 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 50, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())


df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(hub_n_100_P_100[,1], hub_n_100_P_100[,2],
                                    hub_n_100_P_100[,3], hub_n_100_P_100[,4],
                                    hub_n_100_P_100[,5], hub_n_100_P_100[,6],
                                    hub_n_100_P_100[,7], hub_n_100_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p3 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 100, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(hub_n_1000_P_100[,1], hub_n_1000_P_100[,2],
                                    hub_n_1000_P_100[,3], hub_n_1000_P_100[,4],
                                    hub_n_1000_P_100[,5], hub_n_1000_P_100[,6],
                                    hub_n_1000_P_100[,7], hub_n_1000_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p4 <- p + geom_boxplot()+ theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 1000, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

grid.arrange(p1, p2, p3, p4, nrow = 2, ncol = 2) 

Toeplitz matrix

We simulate data from a Toeplitz structured correlation matrix under different settings of \(n\) and \(p\).

First run the code in

source("../code/Figure5/corshrink_toeplitz_pcor.R")  ## n=50, p=100 

## last line : save(frob_vals, file = paste0("toeplitz_cmd_boot_n_", n, "_P_", P, "_results_pcor.rda"))
## change the file to the destination and file name you prefer.

Change the \(n\) and \(p\) accordingly, we pick \(p=100\) and vary \(n\) from \(30, 50, 100, 1000\).

The final output from these files were saved under shared_output/Pcor/ folders and can be read as follows.

toeplitz_n_10_P_100 <- get(load("../shared_output/Pcor/toeplitz_cmd_n_30_P_100_results_pcor.rda"))
toeplitz_n_50_P_100 <- get(load("../shared_output/Pcor/toeplitz_cmd_n_50_P_100_results_pcor.rda"))
toeplitz_n_100_P_100 <- get(load("../shared_output/Pcor/toeplitz_cmd_n_100_P_100_results_pcor.rda"))
toeplitz_n_1000_P_100 <- get(load("../shared_output/Pcor/toeplitz_cmd_n_1000_P_100_results_pcor.rda"))

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(toeplitz_n_10_P_100[,1], toeplitz_n_10_P_100[,2],
                                    toeplitz_n_10_P_100[,3], toeplitz_n_10_P_100[,4],
                                    toeplitz_n_10_P_100[,5], toeplitz_n_10_P_100[,6],
                                    toeplitz_n_10_P_100[,7], toeplitz_n_10_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p5 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 30, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(toeplitz_n_50_P_100[,1], toeplitz_n_50_P_100[,2],
                                    toeplitz_n_50_P_100[,3], toeplitz_n_50_P_100[,4],
                                    toeplitz_n_50_P_100[,5], toeplitz_n_50_P_100[,6],
                                    toeplitz_n_50_P_100[,7], toeplitz_n_50_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p6 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 50, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())


df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(toeplitz_n_100_P_100[,1], toeplitz_n_100_P_100[,2],
                                    toeplitz_n_100_P_100[,3], toeplitz_n_100_P_100[,4],
                                    toeplitz_n_100_P_100[,5], toeplitz_n_100_P_100[,6],
                                    toeplitz_n_100_P_100[,7], toeplitz_n_100_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p7 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 100, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(toeplitz_n_1000_P_100[,1], toeplitz_n_1000_P_100[,2],
                                    toeplitz_n_1000_P_100[,3], toeplitz_n_1000_P_100[,4],
                                    toeplitz_n_1000_P_100[,5], toeplitz_n_1000_P_100[,6],
                                    toeplitz_n_1000_P_100[,7], toeplitz_n_1000_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p8 <- p + geom_boxplot()+ theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 1000, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

grid.arrange(p5, p6, p7, p8, nrow = 2, ncol = 2) 
## Warning: Removed 19 rows containing non-finite values (stat_boxplot).
## Warning: Removed 17 rows containing non-finite values (stat_boxplot).

Banded precision matrix

We simulate data from a Banded Precision structured correlation matrix under different settings of \(n\) and \(p\).

First run the code in

source("../code/Figure5/corshrink_banded_precision_pcor.R")  ## n=50, p=100 

## last line : save(frob_vals, file = paste0("banded_precision_cmd_boot_n_", n, "_P_", P, "_results_pcor.rda"))
## change the file to the destination and file name you prefer.

Change the \(n\) and \(p\) accordingly, we pick \(p=100\) and vary \(n\) from \(30, 50, 100, 1000\).

The final output from these files were saved under shared_output/Pcor/ folders and can be read as follows.

banded_precision_n_10_P_100 <- get(load("../shared_output/Pcor/banded_precision_n_30_P_100_results_pcor.rda"))
banded_precision_n_50_P_100 <- get(load("../shared_output/Pcor/banded_precision_n_50_P_100_results_pcor.rda"))
banded_precision_n_100_P_100 <- get(load("../shared_output/Pcor/banded_precision_n_100_P_100_results_pcor.rda"))
banded_precision_n_1000_P_100 <- get(load("../shared_output/Pcor/banded_precision_n_1000_P_100_results_pcor.rda"))

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(banded_precision_n_10_P_100[,1], banded_precision_n_10_P_100[,2],
                                    banded_precision_n_10_P_100[,3], banded_precision_n_10_P_100[,4],
                                    banded_precision_n_10_P_100[,5], banded_precision_n_10_P_100[,6],
                                    banded_precision_n_10_P_100[,7], banded_precision_n_10_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p9 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 30, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(banded_precision_n_50_P_100[,1], banded_precision_n_50_P_100[,2],
                                    banded_precision_n_50_P_100[,3], banded_precision_n_50_P_100[,4],
                                    banded_precision_n_50_P_100[,5], banded_precision_n_50_P_100[,6],
                                    banded_precision_n_50_P_100[,7], banded_precision_n_50_P_100[,8])))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p10 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 50, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())


df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(banded_precision_n_100_P_100[,1], banded_precision_n_100_P_100[,2],
                                    banded_precision_n_100_P_100[,3], banded_precision_n_100_P_100[,4],
                                    banded_precision_n_100_P_100[,5], banded_precision_n_100_P_100[,6],
                                    banded_precision_n_100_P_100[,7], banded_precision_n_100_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p11 <- p + geom_boxplot() + theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +  ggtitle(paste0("n = ", 100, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

df <- data.frame("method" = factor(c(rep("empirical", 50), rep("corpcor", 50), rep("corshrink", 50),
                              rep("pdsce",50), rep("glasso:1e-02", 50),
                              rep("glasso:0.1", 50), rep("glasso:0.5", 50), rep("glasso:1", 50)), 
                              levels = c("empirical", "corpcor",
                                        "corshrink", "pdsce",
                                        "glasso:1e-02",
                                        "glasso:0.1",
                                        "glasso:0.5",
                                        "glasso:1")),
                 "distance" = log(c(banded_precision_n_1000_P_100[,1], banded_precision_n_1000_P_100[,2],
                                    banded_precision_n_1000_P_100[,3], banded_precision_n_1000_P_100[,4],
                                    banded_precision_n_1000_P_100[,5], banded_precision_n_1000_P_100[,6],
                                    banded_precision_n_1000_P_100[,7], banded_precision_n_1000_P_100[,8]
                                    )))

p <- ggplot(df, aes(method, distance, color = method)) + ylab("log(distance)")
p12 <- p + geom_boxplot()+ theme_bw() + theme(legend.position="none") + 
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + ggtitle(paste0("n = ", 1000, " p = ", 100)) +
scale_y_continuous(breaks= pretty_breaks())

grid.arrange(p9, p10, p11, p12, nrow = 2, ncol = 2) 
## Warning: Removed 14 rows containing non-finite values (stat_boxplot).
## Warning: Removed 17 rows containing non-finite values (stat_boxplot).

p <- ggplot(df, aes(method, distance, color = method)) + xlab("") + ylab("log(distance)")+ ggtitle("n=50, p = 100")
p13 <- p + geom_boxplot() + theme_bw() 
print(p13)

grid.arrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12, nrow = 4, ncol = 3, as.table = FALSE) 
## Warning: Removed 19 rows containing non-finite values (stat_boxplot).
## Warning: Removed 17 rows containing non-finite values (stat_boxplot).
## Warning: Removed 14 rows containing non-finite values (stat_boxplot).
## Warning: Removed 17 rows containing non-finite values (stat_boxplot).


This R Markdown site was created with workflowr