Introduction

workflowr

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/GSE103354_Trachea_droplet_UMIcounts.txt

Untracked files:
    Untracked:  analysis/NBapprox.Rmd
    Untracked:  analysis/pseudocount2.Rmd
    Untracked:  analysis/trachea4.Rmd
    Untracked:  code/missing_data.R
    Untracked:  code/pseudocounts.R
    Untracked:  code/trachea4.R
    Untracked:  data/hard_bimodal1.txt
    Untracked:  data/hard_bimodal2.txt
    Untracked:  data/hard_bimodal3.txt
    Untracked:  docs/figure/pseudocount2.Rmd/
    Untracked:  output/pseudocount/

Unstaged changes:
    Modified:   analysis/pseudocount.Rmd
    Modified:   code/sc_comparisons.R

My primary goal in these analyses is to compare different flashier fits using scRNA datasets. I hope to develop methods and metrics that I can extend to other domains (in particular, population genetics and linguistics). Among the questions I’d like to answer are:

• Which variance structure should be assumed? Should a residual variance parameter be estimated for each gene, or is it enough to estimate a single variance parameter that’s shared across all genes?
• How should the data be transformed? The usual practice is to add a pseudocount and then take logarithms. How does the choice of pseudocount affect the fit?
• How should the data be normalized? Does it suffice to take, for example, some fraction of the library size as a scale factor or should more sophisticated methods be used? Should counts be normalized before or after adding a pseudocount?
• Which prior families should be used? Does the flexibility of normal-mixture priors outweigh the speed of point-normal priors? Or should one prefer the interpretability of semi-nonnegative factorizations?
• Does the usual stopping criterion work when choosing the number of factors to add? If not, when to stop?
• How does backfitting affect the fit? Does it result in “nicer” factors or better predictive accuracy?

Clearly, I will need some metrics for evaluating fits. First, though, I want to take a step back and think about what I want these factor analyses to accomplish. In general, I’d argue that the value of a factor analysis consists in:

• The factors themselves. A given factor tells us something about which genes tend to co-express.
• The factors in so far as they decompose individual samples. A cell’s factor loadings give us a view of that cell’s expression profile in a glance while providing a more nuanced understanding than, for example, cluster membership would.
• The factors as a low-rank representation of the complete data. A factor analysis should be able to do anything that any other matrix factorization (e.g., an SVD) can do.
• The model as a whole. Goodness-of-fit measures should allow us to say something about whether the model assumptions are reasonable.

Since any particular metric will likely favor one of these perspectives over the others, it’d be desirable to have a stable of metrics that acknowledges each perspective. Together, the metrics should be able to assess whether a fit achieves the following goals:

• Factors should be easily interpretable. Usually, this means that they should be both sparse and coherent. Sparsity can be measured either via the estimated priors on the factor loadings or via the local false sign rates of the loadings. “Coherence,” a term that I’ve borrowed from the NLP community (by way of Zihao Wang), is supposed to measure whether it makes sense for a given set of genes to appear together. This kind of metric might be established by comparing gene loadings to known pathways. Usually, it’s possible to align factors from different fits (by simply computing cross-correlations), so the sparsity and coherence of a factor from one fit could be compared to the sparsity and coherence of the “same” factor from another fit.
• Factor loadings should make sense given what we already know about samples. If, for example, cell types are “known” in advance, then the factors should be able to discriminate among those types. (For example, it’d be good to find one or more factors for each cell type that are primarily loaded on that type.) Further, factors should be able to successfully “decompose” samples we haven’t seen yet. A good factor model should be able to explain a large proportion of variance in a hold-out set of samples.
• The “de-noising” implied by representing the data as a low-rank matrix should be sensible. In particular, the distribution of residuals should appear reasonable (this can be checked by plotting \(p\)-values). Since flashier can deal with missing data, one can also compare fits via data imputation tasks. The latter can be tricky for data transformations, however, since data must be imputed on a different scale for each fit.
• The model fit should be good. The variational lower bound gives some measure of goodness of fit, but it’s not perfect. For nested models (comparing, for example, constant and gene-wise variances), some sort of information criterion (AIC or BIC) might be useful.
  
  

  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.5
#> 
#> 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