Simulation set up

The results below are based on simulations with 900 samples, 5000 variables of which 5 were causal, 10 responses with a per-response proportion of variance explained (PVE) of 0.2. Variables, X, were drawn from MVN(0, Gamma), where Gamma is such that it achieves a correlation between variables of 0 and a scale of 1. Causal effects, B, were drawn from MVN(0, Sigma), Sigma is such that it achieves a correlation between responses of 1. Response 1, 2, 3 have causal effects while the remaining seven responses do not have any causal effect. The responses, Y, were drawn from MN(XB, I, V), where V is such that it achieves a correlation between responses of 0 and a scale defined by PVE, in the case of cuasal responses. Non causal responses are from N(0, 1).

2000 such datasets (i.e., “genes”) were simulated and univariate summary statistics were obtained by simple linear regression in the training data (80% of the data. To mirror what would happen in real data analysis, the indexes of the training-test individuals were the same for all the datasets. However, since these 2000 datasets were simulated independently, I do not think it matters). These regression coefficients and standard errors were used as input in the mash pipeline (from Gao) to compute data-driven covariance matrices (up to the ED step included). In particular, the top variable per dataset was used to define a “strong” set and 4 random variables per dataset were used to define a “random” set. Covariance matrices were estimated using flash, PCA (including the top 3 PCs), and the empirical covariance matrix.

The first 50 datasets were used for the prediction analysis. mr.mash was fitted to the training data, updating V (imposing a diagonal structure) and updating the prior weights using EM updates. The mixture prior consisted of components defined by:

• canonical matrices corresponding to different settings of effect sharing/specificity (i.e., singletons, independent, low heterogeneity, medium heterogeneity, high heterogeneity, shared) plus the spike.

• data-driven matrices estimated as described above plus the spike.

• both canonical and data-driven matrices plus the spike.

The covariance matrices were scaled by a grid of values computed from the univariate summary statistics as in the mash paper. The posterior mean of the regression coefficients were initialized to the estimates of the group-LASSO. The mixture weights were initialized with the proportion of zero-coefficients from the group-LASSO estimate as the weight on the spike and the proportion of non-zero-coefficients split equally among the remaining components.Convergence was declared when the maximum difference in the ELBO between two successive iterations was smaller than 1e-2.

Then, responses were predicted on the test data (20% of the data).

Here, we evaluate the accuracy of prediction assessed by \(r^2\) and bias (slope) from the regression of the true response on the predicted response, and the relative root mean square error (rRMSE) scaled by the standard deviation of the true responses in the test data. The boxplots are across the 50 datasets, for each response separately.

mr.mash vs other methods

Here, we compare mr.mash to the multivariate versions of LASSO as implemented in glmnet. The form of the penalty id the following: \(\lambda[(1-\alpha)/2 ||\mathbf{\beta}_j||^2_2 + \alpha ||\mathbf{\beta}_j||_2]\). \(\lambda\) is chosen by cross-validation in the training set. All the methods were using 4 threads – mr.mash loops over the mixture components in parallel, glmnet loops over folds in parallel.

###Convert from list to data.frame for plotting
dsc_plots <- convert_dsc_to_dataframe(dsc_out)

###Compute rmse score (relative to mr_mash_consec_em) and add it to the data
rrmse_dat <- compute_rrmse(dsc_plots)
dsc_plots <- rbind(dsc_plots, rrmse_dat)

###Remove mse from scores and keep only methods wanted
dsc_plots <- dsc_plots[which(dsc_plots$score_metric!="scaled_mse" ), ]
dsc_plots <- dsc_plots[which(dsc_plots$method %in% c("mr_mash_em_can", "mr_mash_em_data",
                                                          "mr_mash_em_dataAndcan", "mlasso")), ]

###Create factor version of method
dsc_plots$method_fac <- factor(dsc_plots$method, levels=c("mr_mash_em_can", "mr_mash_em_data",
                                                          "mr_mash_em_dataAndcan", "mlasso"),
                                                labels=c("mr_mash_can", "mr_mash_data", "mr_mash_both",
                                                         "mlasso"))

###Build data.frame with best accuracy achievable
hlines <- data.frame(score_metric=c("r2", "bias", "rrmse"), max_val=c(unique(dsc_plots$pve), 1, 1))

###Create factor version of response
dsc_plots$response_fac <- as.factor(dsc_plots$response)

###Plot bias
p_bias <- ggplot(dsc_plots[which(dsc_plots$score_metric=="bias"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  geom_boxplot(color = "black", outlier.size = 1, width = 0.85) +
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "Bias", title = "Bias", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="bias"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_bias)

Past versions of “accuracy vs others indepX-1.png”

###Plot r2
p_r2 <- ggplot(dsc_plots[which(dsc_plots$score_metric=="r2"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  geom_boxplot(color = "black", outlier.size = 1, width = 0.85) +
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "r2", title = "r2", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="r2"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_r2)

Past versions of “accuracy vs others indepX-2.png”

###Plot rrmse
p_rrmse <- ggplot(dsc_plots[which(dsc_plots$score_metric=="rrmse"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  geom_boxplot(color = "black", outlier.size = 1, width = 0.85) +
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "rRMSE", title = "rMSE", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="rrmse"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_rrmse)

Past versions of “accuracy vs others indepX-3.png”

Let’s now remove outliers from the plots to make things a little clearer.

###Plot bias
p_bias <- ggplot(dsc_plots[which(dsc_plots$score_metric=="bias"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  Ipaper::geom_boxplot2(color = "black",width.errorbar = 0, width = 0.85) +  
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "Bias", title = "Bias", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="bias"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_bias)

Past versions of “accuracy no outliers indepX-1.png”

###Plot r2
p_r2 <- ggplot(dsc_plots[which(dsc_plots$score_metric=="r2"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  Ipaper::geom_boxplot2(color = "black",width.errorbar = 0, width = 0.85) +  
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "r2", title = "r2", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="r2"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_r2)

Past versions of “accuracy no outliers indepX-2.png”

###Plot rrmse
p_rrmse <- ggplot(dsc_plots[which(dsc_plots$score_metric=="rrmse"), ], aes_string(x = "response_fac", y = "score_value", fill = "method_fac")) +
  Ipaper::geom_boxplot2(color = "black",width.errorbar = 0, width = 0.85) +
  scale_fill_manual(values = colors) +
  labs(x = "Response", y = "rRMSE", title = "rMSE", fill="Method") +
  geom_hline(data=hlines[which(hlines$score_metric=="rrmse"), ], aes(yintercept = max_val), linetype="dashed", size=1) +
  theme_cowplot(font_size = 20) +
  theme(plot.title = element_text(hjust = 0.5))

print(p_rrmse)

Past versions of “accuracy no outliers indepX-3.png”

Here, we look at the elapsed time (\(log_{10}\) seconds) of each method. Note that the mr.mash run time does not include the run time of group-LASSO (but should be considered since we used it to initialize mr.mash).

dsc_plots_time <- dsc_plots[which(dsc_plots$response==1 & dsc_plots$score_metric=="r2"), 
                          -which(colnames(dsc_plots) %in% c("score_metric", "score_value", "response"))]

p_time <- ggplot(dsc_plots_time, aes_string(x = "method_fac", y = "time", fill = "method_fac")) +
  geom_boxplot(color = "black", outlier.size = 1, width = 0.85) +
  scale_fill_manual(values = colors) +
  scale_y_continuous(trans="log10", breaks = trans_breaks("log10", function(x) 10^x),
                      labels = trans_format("log10", math_format(10^.x))) +
  labs(x = "", y = "Elapsed time (seconds) in log10 scale",title = "Run time", fill="Method") +
  theme_cowplot(font_size = 20) +
  theme(axis.line.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank(),
        plot.title = element_text(hjust = 0.5))

print(p_time)

Past versions of “run time indepX-1.png”

In summary, mr.mash with the data-driven covariance matrices does very well and runs pretty fast. Using only the canonical covariance matrices is not as effective (as expected) in this case, but it’s still good. Using both types of covariance matrices does not add anything in terms of performance to using only the data-driven matrices. However, it makes the method much slower to run.

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
[1] devtools_2.1.0  usethis_1.5.1   magrittr_1.5    scales_1.1.1   
[5] cowplot_1.0.0   ggplot2_3.3.2   dscrutils_0.4.2

loaded via a namespace (and not attached):
 [1] pkgload_1.1.0      jsonlite_1.6       foreach_1.4.4     
 [4] Ipaper_0.1.5       assertthat_0.2.1   sp_1.3-1          
 [7] cellranger_1.1.0   yaml_2.2.1         remotes_2.1.0     
[10] progress_1.2.2     sessioninfo_1.1.1  pillar_1.4.4      
[13] backports_1.1.8    lattice_0.20-38    glue_1.4.1        
[16] digest_0.6.25      RColorBrewer_1.1-2 promises_1.0.1    
[19] colorspace_1.4-1   htmltools_0.3.6    httpuv_1.4.5      
[22] plyr_1.8.6         clipr_0.4.1        pkgconfig_2.0.3   
[25] purrr_0.3.3        processx_3.4.2     whisker_0.3-2     
[28] openxlsx_4.1.0     later_0.7.5        git2r_0.26.1      
[31] tibble_3.0.1       farver_2.0.3       ellipsis_0.3.1    
[34] withr_2.2.0        repr_0.17          cli_2.0.2         
[37] crayon_1.3.4       readxl_1.1.0       memoise_1.1.0     
[40] evaluate_0.14      ps_1.3.3           fs_1.3.1          
[43] fansi_0.4.1        doParallel_1.0.14  xml2_1.2.0        
[46] pkgbuild_1.0.8     tools_3.5.1        data.table_1.12.8 
[49] prettyunits_1.1.1  hms_0.5.3          matrixStats_0.56.0
[52] lifecycle_0.2.0    stringr_1.4.0      munsell_0.5.0     
[55] zip_1.0.0          callr_3.4.3        compiler_3.5.1    
[58] rlang_0.4.6        grid_3.5.1         rstudioapi_0.11   
[61] iterators_1.0.10   base64enc_0.1-3    labeling_0.3      
[64] rmarkdown_1.10     boot_1.3-20        testthat_2.3.2    
[67] gtable_0.3.0       codetools_0.2-15   reshape2_1.4.4    
[70] R6_2.4.1           lubridate_1.7.4    knitr_1.20        
[73] dplyr_0.8.0.1      workflowr_1.6.2    rprojroot_1.3-2   
[76] desc_1.2.0         stringi_1.4.6      parallel_3.5.1    
[79] IRdisplay_0.6.1    Rcpp_1.0.5         vctrs_0.3.1       
[82] tidyselect_0.2.5