Last updated: 2019-11-14
Checks: 6 0
Knit directory: scFLASH/
This reproducible R Markdown analysis was created with workflowr (version 1.2.0). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20181103) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: code/initialization/
Ignored: data/Ensembl2Reactome.txt
Ignored: data/droplet.rds
Ignored: data/mus_pathways.rds
Ignored: output/backfit/
Ignored: output/final_montoro/
Ignored: output/lowrank/
Ignored: output/prior_type/
Ignored: output/pseudocount/
Ignored: output/pseudocount_redux/
Ignored: output/size_factors/
Ignored: output/var_type/
Untracked files:
Untracked: analysis/NBapprox.Rmd
Untracked: analysis/trachea4.Rmd
Untracked: code/alt_montoro/
Untracked: code/missing_data.R
Untracked: code/pulseseq/
Untracked: code/trachea4.R
Untracked: output/alt_montoro/
Untracked: output/pulseseq_fit.rds
Unstaged changes:
Modified: analysis/index.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 6a9bd7a | Jason Willwerscheid | 2019-11-14 | wflow_publish(“analysis/elbo_problems.Rmd”) |
I’ve argued that when fitting scRNA-seq data, the log likelihood of the EBMF model can be made infinitely large by taking the pseudocount \(\alpha \to 0\). Here I give a simpler explanation for this phenomenon.
As usual, let \(X\) be a \(n \times p\) matrix of raw counts and let \(Y\) be a \(n \times p\) matrix of transformed counts \[ Y = \log(X + \alpha). \]
As before, I begin with the observation that as \(\alpha \to 0\), \(Y\) is approximately binary, with \(Y_{ij} = \log \alpha\) where \(X_{ij} = 0\) and \(Y_{ij} \approx 0\) where \(X_{ij} > 0\).
Now, how does the (unadjusted) ELBO change as \(\alpha \to 0\)? Fix \(\alpha\) to something small and consider the even smaller pseudocount \(\alpha^c\), where \(c > 1\). The transformed matrix is \[ \tilde{Y} = \log(X + \alpha^c) \approx c Y \] since \(\tilde{Y}_{ij} = c \log \alpha\) where \(X_{ij} = 0\) and \(\tilde{Y}_{ij} \approx 0\) where \(X_{ij} > 0\). That is, \(c\) simply scales the transformed data matrix \(Y\), which doesn’t affect the KL-divergence between priors and posteriors but increases the data log likelihood by \(np \log c\). (This can be calculated directly by noting that scaling the data by a factor of \(c\) scales both the estimated residual standard deviations and expected squared residuals by a factor of \(c\). And indeed, \(np \log c\) is the adjustment that would be obtained using the change-of-variables formula with the scaling transformation.)
Thus the ELBO increases logarithmically as the logged pseudocount goes to \(-\infty\).
But how does the ELBO adjustment obtained using the change-of-variables formula \(\log |f'(x)|\) behave as \(\alpha \to 0\)? Fix \(\alpha\) as above, but now consider \(\alpha / k\) for some \(k > 1\). The ELBO adjustment is: \[ - \sum_{i, j} \log(X_{ij} + \alpha / k) + \sum_{i, j} \log(X_{ij} + \alpha) \approx \sum_{\{i, j: X_{ij} = 0\}} \log k = snp \log k, \] where \(s\) is the sparsity of \(X\) (i.e., the proportion of counts that are nonzero).
Thus the ELBO adjustment increases logarithmically as the (non-logged) pseudocount goes to \(-\infty\). In other words, even though the ELBO actually increases logarithmically as the logged pseudocount goes to \(-\infty\), it is adjusted as if it were increasing linearly, which drives the adjusted ELBO to \(\infty\).
There is, of course, a fundamental problem with using a change-of-variables adjustment with binary data: since the tranformation only really transforms two values, the derivative of the transformation is not defined. My argument shows, I think, that this problem can confound attempts to use transformations on sparse data as well. (I hesitate to generalize to discrete data, because it might be possible to obtain good results by requiring transformations to be smooth and monotonic.)
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
#>
#> loaded via a namespace (and not attached):
#> [1] workflowr_1.2.0 Rcpp_1.0.1 digest_0.6.18 rprojroot_1.3-2
#> [5] backports_1.1.3 git2r_0.25.2 magrittr_1.5 evaluate_0.13
#> [9] stringi_1.4.3 fs_1.2.7 whisker_0.3-2 rmarkdown_1.12
#> [13] tools_3.5.3 stringr_1.4.0 glue_1.3.1 xfun_0.6
#> [17] yaml_2.2.0 compiler_3.5.3 htmltools_0.3.6 knitr_1.22