Set up

# Load required packages
library(mashr); library(ExtremeDeconvolution); library(flashr2); library(mclust)

Loading required package: ashr

Package 'mclust' version 5.4
Type 'citation("mclust")' for citing this R package in publications.

Attaching package: 'mclust'

The following object is masked from 'package:ashr':

    dens

FlashResult = readRDS('../output/Immune.flash2.center.greedy.K10.rds')

Flash Loadings

From Flash, we have \[\tilde{Z} = LF' + E\] where F is \(7 \times 5\), L is \(n \times 5\), E is \(n\times7\).

Suppose the rows of L come from a mixture of multivariate normals with covariances \(U_{1}, \cdots, U_{n}\). \[l_{i} \sim \sum_{m=1}^{n} \pi_{m} N_{5}(\mu_{m}, U_{m})\]

\[F l_{i} \sim \sum_{m=1}^{n} \pi_{m} N_{7}(F \mu_{m}, F U_{m}F')\]

loading = FlashResult$L_flash[,1:5]
mod = Mclust(loading, G=1:50)
saveRDS(mod, '../output/Immune.flash.loading.reduce.rds')

mod = readRDS('../output/Immune.flash.loading.reduce.rds')
summary(mod$BIC)

Best BIC values:
            VVV,28      VVV,31       VVV,25
BIC      -273340.4 -273515.681 -273550.2564
BIC diff       0.0    -175.311    -209.8863

plot(mod, what = "BIC",legendArgs = list(x = "bottomleft"))

Using the model fitted above, we have 28 components.

\[l_{i} \sim \sum_{m=1}^{28} \hat{\pi}_{m} N_{5}(\hat{\mu}_{m}, \hat{U}_{m})\] The estimated covariance structures for LF’ are \(F\hat{U}_{1}F', \cdots, F\hat{U}_{28}F'\).

library(plyr)
U_list = alply(mod$parameters$variance$sigma,3)
Factors = FlashResult$F_flash[,1:5]
U.loading = lapply(U_list, function(U){Factors %*% (U %*% t(Factors))})
names(U.loading) = paste0('Load', "_", (1:length(U.loading)))

saveRDS(U.loading, '../output/Immune.flash.load.reduce.cov.rds')

Session information

sessionInfo()

R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.2

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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] plyr_1.8.4               mclust_5.4              
[3] flashr2_0.4-0            ExtremeDeconvolution_1.3
[5] mashr_0.2-4              ashr_2.1-27             

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.14      knitr_1.17        magrittr_1.5     
 [4] MASS_7.3-47       doParallel_1.0.11 pscl_1.5.2       
 [7] SQUAREM_2017.10-1 lattice_0.20-35   foreach_1.4.4    
[10] stringr_1.2.0     tools_3.4.3       parallel_3.4.3   
[13] grid_3.4.3        rmeta_2.16        git2r_0.20.0     
[16] htmltools_0.3.6   iterators_1.0.9   assertthat_0.2.0 
[19] yaml_2.1.16       rprojroot_1.2     digest_0.6.13    
[22] Matrix_1.2-12     codetools_0.2-15  evaluate_0.10.1  
[25] rmarkdown_1.8     stringi_1.1.6     compiler_3.4.3   
[28] backports_1.1.2   mvtnorm_1.0-6     truncnorm_1.0-7