Comparing flashier fits of GWAS data

Variance structures

The most general EBMF model is: \[ Y = LF' + S + E \] where \(Y \in \mathbb{R}^{n \times p}\), \(S_{ij} \sim N(0, \sigma_{ij}^2)\) and \(E_{ij} \sim N(0, 1 / \tau_{ij})\), with the \(\sigma_{ij}^2\)s known and the \(\tau_{ij}\)s to be estimated (I ignore the priors on \(L\) and \(F\) here). One can either include \(S\) or not; and one can make different assumptions about how \(\tau\) is structured. I test eight variance structures:

• “constant”: \(S = 0\) and \(\tau_{ij} = \tau\) (1 variance parameter needs to be estimated)
• “by-row”: \(S = 0\) and \(\tau_{ij} = \tau_i\) (\(n\) variance parameters)
• “by-column”: \(S = 0\) and \(\tau_{ij} = \tau_j\) (\(p\) variance parameters)
• “Kronecker”: \(S = 0\) and \(\tau_{ij} = \tau_i^{(1)} \tau_j^{(2)}\) (\(n + p\) variance parameters)
• “zero”: only \(S\) is used; no additional variance is estimated
• “noisy (constant)”: \(S\) is used and \(\tau_{ij} = \tau\) (1 variance parameter)
• “noisy (by-row)”: \(S\) is used and \(\tau_{ij} = \tau_i\) (\(n\) variance parameters)
• “noisy (by-column)”: \(S\) is used and \(\tau_{ij} = \tau_j\) (\(p\) variance parameters)

Initialization methods

• “flashier”: the default method adds factors one at a time. Each factor is initialized by finding the best rank-one approximation to the matrix of residuals.
• “softImpute”: uses package softImpute to initialize factors. This gives different results when there is missing data.
• “initialize from data”: adds a bunch of factors all at once using softImpute and then refines the fit via backfitting.

Code

The code used to produce the fits can be viewed here.

Results

For each fit, I give the ELBO relative to the best overall ELBO attained and (in parentheses) the number of factors added.

res <- readRDS("./output/jean/flashier_res.rds")

calls <- lapply(res, `[[`, "call")

var.type <- sapply(calls, function(call) {
  if (is.null(call$S)) {
    return(switch(paste(as.character(call$var.type), collapse = " "),
                  "0" = "constant",
                  "1" = "by.row",
                  "2" = "by.col",
                  "1 2" = "kronecker"))
  } else {
    return(switch(paste(as.character(call$var.type), collapse = " "),
                  "0" = "noisy.const",
                  "1" = "noisy.byrow",
                  "2" = "noisy.bycol",
                  "1 2" = "noisy.kronecker",
                  "zero"))
  }
})

init.type <- sapply(calls, function(call) {
  if (!is.null(call$init.fn))
    return("soft.impute")
  else if (!is.null(call$EF.init))
    return("init.from.data")
  else
    return("flashier")
})

best.elbo <- max(sapply(res, `[[`, "elbo"))

display.str <- sapply(res, function(fl) {
  return(paste0(formatC(fl$elbo - best.elbo, format = "f", digits = 0), " (",
                sum(fl$pve > 0), ")"))
})

df <- data.frame(var.type = factor(var.type,
                                   levels = c("constant", "by.row", "by.col",
                                              "kronecker", "zero", "noisy.const",
                                              "noisy.byrow", "noisy.bycol")),
                 init.type = factor(init.type,
                                    levels = c("flashier", "soft.impute",
                                               "init.from.data")),
                 display.str = display.str)

knitr::kable(df %>% spread(init.type, display.str))

Some observations:

• softImpute typically does “better” (judging by the ELBO) than the default method.
• “Initialize from data” does better than either, but leaves a lot of factors in place.
  
  

  Session information

sessionInfo()

#> R version 3.5.3 (2019-03-11)
#> Platform: x86_64-apple-darwin15.6.0 (64-bit)
#> Running under: macOS Mojave 10.14.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] mixsqp_0.1-119  ashr_2.2-38     ebnm_0.1-24     flashier_0.1.15
#> [5] tidyr_0.8.3    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.1        highr_0.8         pillar_1.3.1     
#>  [4] compiler_3.5.3    git2r_0.25.2      workflowr_1.2.0  
#>  [7] iterators_1.0.10  tools_3.5.3       digest_0.6.18    
#> [10] evaluate_0.13     tibble_2.1.1      lattice_0.20-38  
#> [13] pkgconfig_2.0.2   rlang_0.3.1       Matrix_1.2-15    
#> [16] foreach_1.4.4     yaml_2.2.0        parallel_3.5.3   
#> [19] xfun_0.6          stringr_1.4.0     knitr_1.22       
#> [22] fs_1.2.7          tidyselect_0.2.5  rprojroot_1.3-2  
#> [25] grid_3.5.3        glue_1.3.1        rmarkdown_1.12   
#> [28] purrr_0.3.2       magrittr_1.5      whisker_0.3-2    
#> [31] backports_1.1.3   codetools_0.2-16  htmltools_0.3.6  
#> [34] MASS_7.3-51.1     stringi_1.4.3     doParallel_1.0.14
#> [37] pscl_1.5.2        truncnorm_1.0-8   SQUAREM_2017.10-1
#> [40] crayon_1.3.4