Grades of membership model: choose random seed

Summary

Motivation: when running CountClut, sample cluster assignment may differ between different random seeds, despite similar cluster characteristics in feature abundancce. This issue is not unique to CountClust. In model-based unsupervised clustering methods, the parameter space is usually complex, and as a result, difficult to find global maximum of the model likelihood.

I run CountClust under 100 different random seeds and then compare Bayes Factors across different random seeds and sample cluster membership. I will use the random seed that generates the largest Bayes Factor in all subsequent CountClust analysis.

This document shows the analysis steps for when K=2. Results show that for when K=2, across 100 random seeds 409 produced the largest BF and is thus chosen as our random seed. This chosen random seed (best seed) is then used in our final results (https://jhsiao999.github.io/nasalmicrobiome/step3_gradesOfMembership.html).

The fitting results for K=2 to 5 under the same 100 random seeds are saved in output/gom-seeds-k-2.rds, output/gom-seeds-k-3.rds, output/gom-seeds-k-4.rds, and output/gom-seeds-k-5.rds.

Package

library(knitr)
library(kableExtra)
library(dplyr)
library(metagenomeSeq)
library(CountClust)
library(parallel)

Loading data

Read in filtered data.

MRobj = readRDS("../data/nasal_filtered_normed_batchcorrected.rds")
MRobj

MRexperiment (storageMode: environment)
assayData: 553 features, 197 samples 
  element names: counts 
protocolData: none
phenoData
  sampleNames: EM0042 EM0047 ... E0168 (197 total)
  varLabels: StudyID age ... infnone (94 total)
  varMetadata: labelDescription
featureData
  featureNames:
    k__Bacteria;p__[Thermi];c__Deinococci;o__Deinococcales;f__Deinococcaceae;g__CM44;s__
    k__Bacteria;p__[Thermi];c__Deinococci;o__Deinococcales;f__Deinococcaceae;g__Deinococcus;s__
    ... s__zeae (553 total)
  fvarLabels: Kingdom Phylum ... Species (7 total)
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:  

K=2

100 random seeds

Apply CountClust to fit grades of membership model.

counts <- MRcounts(MRobj,norm=FALSE,log=FALSE)

# get 100 random seeds
source("../code/print-prime.R")
seeds <- prime(543)

res <- vector("list", length(seeds))

for (i in 1:100) {
  set.seed(seeds[i])
  res[[i]] <- FitGoM(t(counts), K=c(2), tol=1e-6)[[1]]
}

saveRDS(res, file = "../data/res.rds")

res <- readRDS(file="../data/res.rds")

Bayes factor differs between the seeds.

summary(do.call(c,lapply(res, "[[", 4)))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
1085545 1379648 1483781 1489853 1601809 1738210 

Choose the seed with maximum BF: 67

which.max(lapply(res, "[[", 4))

[1] 19

saveRDS(res[[19]], file = "../output/gom-seeds.Rmd/res_best_seed_67.rds")

\(~\)

Results using the seed with the highest BF.

res_best <- readRDS(file = "../output/gom-seeds.Rmd/res_best_seed_67.rds")

cluster <- apply(res_best$omega, 1, which.max)

table(cluster)

cluster
  1   2 
167  30 

kable(res_best$omega)

K=3 to 5

Code in code/gom-seeds.R.

clust_2 <- readRDS("../output/gom-seeds-k-2.rds")
clust_3 <- readRDS("../output/gom-seeds-k-3.rds")
clust_4 <- readRDS("../output/gom-seeds-k-4.rds")
clust_5 <- readRDS("../output/gom-seeds-k-5.rds")

Identify the best seed.

K=2: 67; K=3: 67, K=4: 59; K=5: 307.

seeds[which.max(lapply(clust_2, "[[", 4))]

[1] 67

seeds[which.max(lapply(clust_3, "[[", 4))]

[1] 67

seeds[which.max(lapply(clust_4, "[[", 4))]

[1] 59

seeds[which.max(lapply(clust_5, "[[", 4))]

[1] 307

save the best seeds.

saveRDS(clust_2[which.max(lapply(clust_2, "[[", 4))],
        "../output/gom-k2-best-seed-67.rds")
saveRDS(clust_3[which.max(lapply(clust_3, "[[", 4))],
        "../output/gom-k3-best-seed-67.rds")
saveRDS(clust_4[which.max(lapply(clust_4, "[[", 4))],
        "../output/gom-k4-best-seed-59.rds")
saveRDS(clust_5[which.max(lapply(clust_5, "[[", 4))],
        "../output/gom-k5-best-seed-307.rds")

Session Info

R version 3.4.1 (2017-06-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13

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] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] CountClust_1.4.1     ggplot2_2.2.1        metagenomeSeq_1.21.1
 [4] RColorBrewer_1.1-2   glmnet_2.0-16        foreach_1.4.4       
 [7] Matrix_1.2-14        limma_3.34.9         Biobase_2.38.0      
[10] BiocGenerics_0.24.0  dplyr_0.7.4          kableExtra_0.8.0    
[13] knitr_1.20          

loaded via a namespace (and not attached):
 [1] httr_1.3.1         viridisLite_0.3.0  gtools_3.5.0      
 [4] assertthat_0.2.0   highr_0.6          stats4_3.4.1      
 [7] yaml_2.1.19        slam_0.1-43        pillar_1.2.2      
[10] backports_1.1.2    lattice_0.20-35    glue_1.2.0        
[13] digest_0.6.15      rvest_0.3.2        colorspace_1.3-2  
[16] picante_1.6-2      cowplot_0.9.2      htmltools_0.3.6   
[19] plyr_1.8.4         pkgconfig_2.0.1    scales_0.5.0      
[22] gdata_2.18.0       git2r_0.21.0       tibble_1.4.2      
[25] mgcv_1.8-23        nnet_7.3-12        lazyeval_0.2.1    
[28] magrittr_1.5       evaluate_0.10.1    nlme_3.1-137      
[31] MASS_7.3-50        gplots_3.0.1       xml2_1.2.0        
[34] vegan_2.5-1        tools_3.4.1        hms_0.4.2         
[37] matrixStats_0.53.1 stringr_1.3.0      munsell_0.4.3     
[40] cluster_2.0.7-1    bindrcpp_0.2.2     compiler_3.4.1    
[43] caTools_1.17.1     rlang_0.2.0        grid_3.4.1        
[46] iterators_1.0.9    rstudioapi_0.7     bitops_1.0-6      
[49] rmarkdown_1.9      boot_1.3-20        gtable_0.2.0      
[52] codetools_0.2-15   flexmix_2.3-14     reshape2_1.4.3    
[55] R6_2.2.2           bindr_0.1.1        rprojroot_1.3-2   
[58] maptpx_1.9-3       KernSmooth_2.23-15 permute_0.9-4     
[61] readr_1.1.1        ape_5.1            modeltools_0.2-21 
[64] stringi_1.2.2      SQUAREM_2017.10-1  Rcpp_0.12.16