Comparison of linear regression methods in the four scenarios from Zou & Hastie (2005)

Load packages

Load a few packages and custom functions used in the analysis below.

library(dscrutils)
library(ggplot2)
library(cowplot)
source("../code/plots.R")

Import DSC results

Here we use function “dscquery” from the dscrutils package to extract the DSC results we are interested in—the mean squared error in the predictions from each method and in each simulation scenario.

methods <- c("ridge","lasso","elastic_net","susie","varbvs","varbvsmix")
dsc <- dscquery("../dsc/linreg",
                c("simulate","fit","mse.err","simulate.scenario"),
                verbose = FALSE)
dsc <- transform(dsc,
                 simulate          = factor(simulate),
                 simulate.scenario = factor(simulate.scenario),
                 fit               = factor(fit,methods))
names(dsc)[1] <- "seed"
nrow(dsc)
# [1] 7200

After this call, the “dsc” data frame should contain results for 7,200 pipelines—6 methods times 6 scenarios times 200 data sets simulated in each scenario.

nrow(dsc)
# [1] 7200

After these steps, the “dsc” data frame should have five columns: “seed”, the seed used to simulate the data; “simulate”, the data simulation module used; “simulate.scenario”, the particular scenario used in the “zh” simulation module; “fit”, the linear regression method used; and “mse.err”, the mean squared error in the test set predictions.

head(dsc)
#   seed     simulate   fit   mse.err simulate.scenario
# 1    1 null_effects ridge 12.322653              <NA>
# 2    2 null_effects ridge 12.695073              <NA>
# 3    3 null_effects ridge 10.625812              <NA>
# 4    4 null_effects ridge 10.267211              <NA>
# 5    5 null_effects ridge 10.429469              <NA>
# 6    6 null_effects ridge  9.940155              <NA>

Note that you will need to run the DSC before querying the results; see here for instructions on running the DSC. If you did not run the DSC to generate these results, you can replace the dscquery call above by this line to load the pre-extracted results stored in a CSV file:

dsc <- read.csv("../output/linreg_mse.csv")

This is how the CSV file was created:

write.csv(dsc,"../output/linreg_mse.csv",row.names = FALSE,quote = FALSE)

Summarize and discuss simulation results

Compute the mean squared error (MSE) in the predictions relative to ridge regression, so that larger numbers mean greater error relative to predictions from ridge regressions.

rmse <- compute.relative.mse(dsc)
dsc  <- cbind(dsc,rmse)

The boxplots below summarize the prediction errors in each of the simulations. A relative MSE less than 1 indicates an improvement in accuracy over ridge regression, whereas a relative MSE greater than 1 indicates a decrease in accuracy compared to ridge regression. Ridge regression will always have a relative MSE of 1, so the results for ridge regression are not shown.

p1 <- rmse.boxplot(subset(dsc,simulate == "null_effects"),"null")
p2 <- rmse.boxplot(subset(dsc,simulate == "one_effect"),"one effect")
p3 <- rmse.boxplot(subset(dsc,simulate.scenario == 1),"scenario 1")
p4 <- rmse.boxplot(subset(dsc,simulate.scenario == 2),"scenario 2")
p5 <- rmse.boxplot(subset(dsc,simulate.scenario == 3),"scenario 3")
p6 <- rmse.boxplot(subset(dsc,simulate.scenario == 4),"scenario 4")
p  <- plot_grid(p1,p2,p3,p4,p5,p6)
print(p)

Past versions of create-boxplots-1.png

Here are a few initial impressions from these plots.

In most cases, the Elastic Net does at least as well, or better, than the Lasso. This is what we would expect.

Ridge regression actually achieves good accuracy in all cases except Scenario 4 and the “one effect” setting. Ridge regression is expected to do less well in Scenario 4 because the majority of the true coefficients are zero, so a sparse model would be favoured. Similarly, a sparse model should better fit data simulated in the “one effect” scenario.

In Scenario 4, where the predictors are correlated in a structured way, and the effects are sparse, varbvs and varbvsmix perform considerably better than the other methods.

The “varbvsmix” method yields competitive predictions in all scenarios.

Session information

The “Session information” button below gives the version of R and the packages that were used to generate these results. This listing includes the R packages that were also used to run the DSC.

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] C
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] varbvs_2.6-5      susieR_0.9.0.0580 glmnet_2.0-16     foreach_1.4.4    
# [5] Matrix_1.2-15     MASS_7.3-51.1     cowplot_0.9.4     ggplot2_3.1.1    
# [9] dscrutils_0.4.2  
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.2          RColorBrewer_1.1-2  compiler_3.5.1     
#  [4] pillar_1.3.1        later_0.7.5         git2r_0.26.1       
#  [7] plyr_1.8.4          workflowr_1.6.0     iterators_1.0.10   
# [10] prettyunits_1.0.2   tools_3.5.1         progress_1.2.0     
# [13] digest_0.6.18       lattice_0.20-38     evaluate_0.12      
# [16] tibble_2.1.1        gtable_0.2.0        pkgconfig_2.0.2    
# [19] rlang_0.4.0         yaml_2.2.0          withr_2.1.2        
# [22] stringr_1.3.1       dplyr_0.8.0.1       knitr_1.20         
# [25] fs_1.3.1            hms_0.4.2           tidyselect_0.2.5   
# [28] rprojroot_1.3-2     grid_3.5.1          glue_1.3.0         
# [31] data.table_1.12.0   R6_2.3.0            rmarkdown_1.10     
# [34] latticeExtra_0.6-28 purrr_0.3.2         magrittr_1.5       
# [37] whisker_0.3-2       codetools_0.2-15    backports_1.1.2    
# [40] scales_1.0.0        promises_1.0.1      htmltools_0.3.6    
# [43] assertthat_0.2.0    colorspace_1.3-2    httpuv_1.4.5       
# [46] labeling_0.3        nor1mix_1.2-3       stringi_1.2.4      
# [49] lazyeval_0.2.1      munsell_0.5.0       crayon_1.3.4