Compare_ebpmf_nmf2

Description

I compare ebpmf with nmf algorithms on 10xgenomics dataset cd14_monocytes (\(\text{n-sample} = 2611, \text{n-feature} = 359\)). Training set \(X_{ij} = Binomial(data_{ij}, 0.5)\), and validation set is \(Y_{ij} = data_{ij} - X{ij}\).

Summary of results:

• It is disconcerting that among ebpmf methods, the ELBO does not correspond to the validation loglikelihood: ebpmf-point-gamma has much higher ELBO than ebpmf-exponential-mixture, but the validation loglikelihood is worse. Will investigaet into that
  
• ebpmf does better than nnmf on validation set.
  
• I accidently made a mistake that gives iteresting result in using NNLM::nnmf: I initialized it using init = list(W0 = ..., H0 = ...) and got the best training and validation loglikelihood among all algorithms. This specification means fitting \[ A \approx WH + W_0 H_1 + W_1 H_0 \] in https://cran.r-project.org/web/packages/NNLM/vignettes/Fast-And-Versatile-NMF.pdf section Content deconvolution and designable factorization

devtools::load_all("../ebpmf")
devtools::load_all("../ebpm")
library(ebpmf)
library(gtools)
library(NNLM)
library(ggplot2)

dataset

10X genomics dataset

X = read.csv("data/10xgenomics/cd14_monocytes/filtered_matrices_mex/hg19/Y.csv")
Y = read.csv("data/10xgenomics/cd14_monocytes/filtered_matrices_mex/hg19/Yhat.csv")
X = as.matrix(X)
Y = as.matrix(Y)
rownames(X)  = NULL
colnames(X)  = NULL
rownames(Y)  = NULL
colnames(Y)  = NULL

real = list(X = as.matrix(X), Y = as.matrix(Y))
print(dim(real$X))

[1] 2611  359

hist(real$X, breaks = 100, main = "hist for Y_train")

Past versions of unnamed-chunk-2-1.png

algorithm setup

K = 3
maxiter = 1000
method = c()
ll_train = c()
ll_val = c()
out = list()

run experiments

out = readRDS("data/Compare_ebpmf_nmf2_out.Rds")

## Run ebpmf_exponential_mixture
method_ = "ebpmf_exponential_mixture"
# res = ebpmf::ebpmf_exponential_mixture(X, K, m = 2^0.25, maxiter.out = maxiter)
res = out$ebpmf_exponential_mixture
Lam = res$qg$qls_mean %*% t(res$qg$qfs_mean)
method = c(method, method_)
ll_train = c(ll_train, sum(dpois(X, Lam, log = T)))
ll_val   = c(ll_val, sum(dpois(Y, Lam, log = T)))
out[[method_]] = res
#plot(1:maxiter, res$ELBO, main = sprintf("(maximized) objective for %s", method_), xlab = "iter", ylab = "ELBO")

## Run ebpmf_point_gamma
method_ = "ebpmf_point_gamma"
# res = ebpmf::ebpmf_point_gamma(X, K, maxiter.out = maxiter)
res = out$ebpmf_point_gamma
Lam = res$qg$qls_mean %*% t(res$qg$qfs_mean)
method = c(method, method_)
ll_train = c(ll_train, sum(dpois(X, Lam, log = T)))
ll_val   = c(ll_val, sum(dpois(Y, Lam, log = T)))
out[[method_]] = res
#plot(1:maxiter, res$ELBO, main = sprintf("(maximized) objective for %s", method_), xlab = "iter", ylab = "ELBO")

## Run nnmf
method_ = "nnmf"
#res = NNLM::nnmf(A = X, k = K, init = list(W0 = res$qg$qls_mean, H0 = t(res$qg$qfs_mean)), loss = "mkl", method = "lee", max.iter = maxiter)
res = NNLM::nnmf(A = X, k = K, loss = "mkl", method = "lee", max.iter = maxiter)
Lam = res$W %*% res$H
method = c(method, method_)
ll_train = c(ll_train, sum(dpois(X, Lam, log = T)))
ll_val   = c(ll_val, sum(dpois(Y, Lam, log = T)))
out[[method_]] = res

## Run nnmf with  initialization from ebpmf_point_gamma
method_ = "nnmf_init_wh_from_point_gamma"
res = out$ebpmf_point_gamma
res = NNLM::nnmf(A = X, k = K, init = list(W = res$qg$qls_mean, H = t(res$qg$qfs_mean)), loss = "mkl", method = "lee", max.iter = maxiter)
Lam = res$W %*% res$H
method = c(method, method_)
ll_train = c(ll_train, sum(dpois(X, Lam, log = T)))
ll_val   = c(ll_val, sum(dpois(Y, Lam, log = T)))
out[[method_]] = res

## Run nnmf with  initialization from ebpmf_point_gamma
method_ = "nnmf_init_w0h0_from_point_gamma"
res = out$ebpmf_point_gamma
res = NNLM::nnmf(A = X, k = K, init = list(W0 = res$qg$qls_mean, H0 = t(res$qg$qfs_mean)), loss = "mkl", method = "lee", max.iter = maxiter)
Lam = res$W %*% res$H
method = c(method, method_)
ll_train = c(ll_train, sum(dpois(X, Lam, log = T)))
ll_val   = c(ll_val, sum(dpois(Y, Lam, log = T)))
out[[method_]] = res

Compare loglikelihood in training and validation dataset.

data.frame(method = method, ll_train = ll_train, ll_val = ll_val)

                           method  ll_train    ll_val
1       ebpmf_exponential_mixture -968660.9 -978453.3
2               ebpmf_point_gamma -969986.0 -979553.7
3                            nnmf -969532.4 -982661.4
4   nnmf_init_wh_from_point_gamma -968021.5 -981134.6
5 nnmf_init_w0h0_from_point_gamma -954207.7 -972323.3

show objectives

elbos = data.frame(iters = 1:length(out$ebpmf_exponential_mixture$ELBO),
                   ebpm_exponential_mixture = out$ebpmf_exponential_mixture$ELBO, 
                   ebpm_point_gamma = out$ebpmf_point_gamma$ELBO)
ggplot(elbos)+
  geom_line(aes(x = iters, y = ebpm_exponential_mixture, color = "ebpm_exponential_mixture"), show.legend = T)+
  geom_line(aes(x = iters, y = ebpm_point_gamma, color = "ebpm_point_gamma"), show.legend = T)+
  xlab("iter")+
  ylab("ELBO")+
  theme(legend.position="top")

Past versions of unnamed-chunk-6-1.png

Save results

saveRDS(out, "data/Compare_ebpmf_nmf2_out_ver2.Rds")

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

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] ggplot2_3.2.1   NNLM_0.4.2      gtools_3.8.1    ebpm_0.0.0.9001
[5] ebpmf_0.1.0     testthat_2.2.1 

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5    xfun_0.8            remotes_2.1.0      
 [4] purrr_0.3.2         colorspace_1.4-1    usethis_1.5.1      
 [7] htmltools_0.3.6     yaml_2.2.0          rlang_0.4.0        
[10] pkgbuild_1.0.3      mixsqp_0.1-121      pillar_1.4.2       
[13] glue_1.3.1          withr_2.1.2         sessioninfo_1.1.1  
[16] stringr_1.4.0       munsell_0.5.0       gtable_0.3.0       
[19] workflowr_1.4.0     devtools_2.2.1.9000 memoise_1.1.0      
[22] evaluate_0.14       labeling_0.3        knitr_1.25         
[25] callr_3.2.0         ps_1.3.0            Rcpp_1.0.2         
[28] backports_1.1.5     scales_1.0.0        desc_1.2.0         
[31] pkgload_1.0.2       fs_1.3.1            digest_0.6.22      
[34] stringi_1.4.3       processx_3.3.1      dplyr_0.8.1        
[37] rprojroot_1.3-2     grid_3.5.1          cli_1.1.0          
[40] tools_3.5.1         magrittr_1.5        lazyeval_0.2.2     
[43] tibble_2.1.3        crayon_1.3.4        whisker_0.3-2      
[46] pkgconfig_2.0.3     ellipsis_0.3.0      prettyunits_1.0.2  
[49] assertthat_0.2.1    rmarkdown_1.13      rstudioapi_0.10    
[52] R6_2.4.0            git2r_0.25.2        compiler_3.5.1