Behaviour of Poisson NMF applied to data generated from the Poisson log1p NMF model
Peter Carbonetto
Last updated: 2026-07-21
Checks: 7 0
Knit directory: pcarbo/analysis/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks 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(1) 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 job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 0465dcc. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
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: analysis/fit-scd-ex-k=10.rds
Ignored: analysis/train.mat
Ignored: sbatch/makefile_demo/sims/sim0001.csv
Ignored: sbatch/makefile_demo/sims/sim0002.csv
Ignored: sbatch/makefile_demo/sims/sim0003.csv
Ignored: sbatch/makefile_demo/sims/sim0004.csv
Ignored: sbatch/makefile_demo/sims/sim0005.csv
Ignored: sbatch/makefile_demo/sims/sim0006.csv
Ignored: sbatch/makefile_demo/sims/sim0007.csv
Ignored: sbatch/makefile_demo/sims/sim0008.csv
Ignored: sbatch/makefile_demo/sims/sim0009.csv
Ignored: sbatch/makefile_demo/sims/sim0010.csv
Ignored: sbatch/makefile_demo/sims/sim0011.csv
Ignored: sbatch/makefile_demo/sims/sim0012.csv
Ignored: sbatch/makefile_demo/sims/sim0013.csv
Ignored: sbatch/makefile_demo/sims/sim0014.csv
Ignored: sbatch/makefile_demo/sims/sim0015.csv
Ignored: sbatch/makefile_demo/sims/sim0016.csv
Ignored: sbatch/makefile_demo/sims/sim0017.csv
Ignored: sbatch/makefile_demo/sims/sim0018.csv
Ignored: sbatch/makefile_demo/sims/sim0019.csv
Ignored: sbatch/makefile_demo/sims/sim0020.csv
Ignored: sbatch/makefile_demo/sims/sim0021.csv
Ignored: sbatch/makefile_demo/sims/sim0022.csv
Ignored: sbatch/makefile_demo/sims/sim0023.csv
Ignored: sbatch/makefile_demo/sims/sim0024.csv
Ignored: sbatch/makefile_demo/sims/sim0025.csv
Ignored: sbatch/makefile_demo/sims/sim0026.csv
Ignored: sbatch/makefile_demo/sims/sim0027.csv
Ignored: sbatch/makefile_demo/sims/sim0028.csv
Ignored: sbatch/makefile_demo/sims/sim0029.csv
Ignored: sbatch/makefile_demo/sims/sim0030.csv
Ignored: sbatch/makefile_demo/sims/sim0031.csv
Ignored: sbatch/makefile_demo/sims/sim0032.csv
Ignored: sbatch/makefile_demo/sims/sim0033.csv
Ignored: sbatch/makefile_demo/sims/sim0034.csv
Ignored: sbatch/makefile_demo/sims/sim0035.csv
Ignored: sbatch/makefile_demo/sims/sim0036.csv
Ignored: sbatch/makefile_demo/sims/sim0037.csv
Ignored: sbatch/makefile_demo/sims/sim0038.csv
Ignored: sbatch/makefile_demo/sims/sim0039.csv
Ignored: sbatch/makefile_demo/sims/sim0040.csv
Ignored: sbatch/makefile_demo/sims/sim0041.csv
Ignored: sbatch/makefile_demo/sims/sim0042.csv
Ignored: sbatch/makefile_demo/sims/sim0043.csv
Ignored: sbatch/makefile_demo/sims/sim0044.csv
Ignored: sbatch/makefile_demo/sims/sim0045.csv
Ignored: sbatch/makefile_demo/sims/sim0046.csv
Ignored: sbatch/makefile_demo/sims/sim0047.csv
Ignored: sbatch/makefile_demo/sims/sim0048.csv
Ignored: sbatch/makefile_demo/sims/sim0049.csv
Ignored: sbatch/makefile_demo/sims/sim0050.csv
Ignored: sbatch/makefile_demo/sims/sim0051.csv
Ignored: sbatch/makefile_demo/sims/sim0052.csv
Ignored: sbatch/makefile_demo/sims/sim0053.csv
Ignored: sbatch/makefile_demo/sims/sim0054.csv
Ignored: sbatch/makefile_demo/sims/sim0055.csv
Ignored: sbatch/makefile_demo/sims/sim0056.csv
Ignored: sbatch/makefile_demo/sims/sim0057.csv
Ignored: sbatch/makefile_demo/sims/sim0058.csv
Ignored: sbatch/makefile_demo/sims/sim0059.csv
Ignored: sbatch/makefile_demo/sims/sim0060.csv
Ignored: sbatch/makefile_demo/sims/sim0061.csv
Ignored: sbatch/makefile_demo/sims/sim0062.csv
Ignored: sbatch/makefile_demo/sims/sim0063.csv
Ignored: sbatch/makefile_demo/sims/sim0064.csv
Ignored: sbatch/makefile_demo/sims/sim0065.csv
Ignored: sbatch/makefile_demo/sims/sim0066.csv
Ignored: sbatch/makefile_demo/sims/sim0067.csv
Ignored: sbatch/makefile_demo/sims/sim0068.csv
Ignored: sbatch/makefile_demo/sims/sim0069.csv
Ignored: sbatch/makefile_demo/sims/sim0070.csv
Ignored: sbatch/makefile_demo/sims/sim0071.csv
Ignored: sbatch/makefile_demo/sims/sim0072.csv
Ignored: sbatch/makefile_demo/sims/sim0073.csv
Ignored: sbatch/makefile_demo/sims/sim0074.csv
Ignored: sbatch/makefile_demo/sims/sim0075.csv
Ignored: sbatch/makefile_demo/sims/sim0076.csv
Ignored: sbatch/makefile_demo/sims/sim0077.csv
Ignored: sbatch/makefile_demo/sims/sim0078.csv
Ignored: sbatch/makefile_demo/sims/sim0079.csv
Ignored: sbatch/makefile_demo/sims/sim0080.csv
Ignored: sbatch/makefile_demo/sims/sim0081.csv
Ignored: sbatch/makefile_demo/sims/sim0082.csv
Ignored: sbatch/makefile_demo/sims/sim0083.csv
Ignored: sbatch/makefile_demo/sims/sim0084.csv
Ignored: sbatch/makefile_demo/sims/sim0085.csv
Ignored: sbatch/makefile_demo/sims/sim0086.csv
Ignored: sbatch/makefile_demo/sims/sim0087.csv
Ignored: sbatch/makefile_demo/sims/sim0088.csv
Ignored: sbatch/makefile_demo/sims/sim0089.csv
Ignored: sbatch/makefile_demo/sims/sim0090.csv
Ignored: sbatch/makefile_demo/sims/sim0091.csv
Ignored: sbatch/makefile_demo/sims/sim0092.csv
Ignored: sbatch/makefile_demo/sims/sim0093.csv
Ignored: sbatch/makefile_demo/sims/sim0094.csv
Ignored: sbatch/makefile_demo/sims/sim0095.csv
Ignored: sbatch/makefile_demo/sims/sim0096.csv
Ignored: sbatch/makefile_demo/sims/sim0097.csv
Ignored: sbatch/makefile_demo/sims/sim0098.csv
Ignored: sbatch/makefile_demo/sims/sim0099.csv
Ignored: sbatch/makefile_demo/sims/sim0100.csv
Ignored: sbatch/makefile_demo/sims/sim0101.csv
Ignored: sbatch/makefile_demo/sims/sim0102.csv
Ignored: sbatch/makefile_demo/sims/sim0103.csv
Ignored: sbatch/makefile_demo/sims/sim0104.csv
Ignored: sbatch/makefile_demo/sims/sim0105.csv
Ignored: sbatch/makefile_demo/sims/sim0106.csv
Ignored: sbatch/makefile_demo/sims/sim0107.csv
Ignored: sbatch/makefile_demo/sims/sim0108.csv
Ignored: sbatch/makefile_demo/sims/sim0109.csv
Ignored: sbatch/makefile_demo/sims/sim0110.csv
Ignored: sbatch/makefile_demo/sims/sim0111.csv
Ignored: sbatch/makefile_demo/sims/sim0112.csv
Ignored: sbatch/makefile_demo/sims/sim0113.csv
Ignored: sbatch/makefile_demo/sims/sim0114.csv
Ignored: sbatch/makefile_demo/sims/sim0115.csv
Ignored: sbatch/makefile_demo/sims/sim0116.csv
Ignored: sbatch/makefile_demo/sims/sim0117.csv
Ignored: sbatch/makefile_demo/sims/sim0118.csv
Ignored: sbatch/makefile_demo/sims/sim0119.csv
Ignored: sbatch/makefile_demo/sims/sim0120.csv
Ignored: sbatch/makefile_demo/sims/sim0121.csv
Ignored: sbatch/makefile_demo/sims/sim0122.csv
Ignored: sbatch/makefile_demo/sims/sim0123.csv
Ignored: sbatch/makefile_demo/sims/sim0124.csv
Ignored: sbatch/makefile_demo/sims/sim0125.csv
Ignored: sbatch/makefile_demo/sims/sim0126.csv
Ignored: sbatch/makefile_demo/sims/sim0127.csv
Ignored: sbatch/makefile_demo/sims/sim0128.csv
Ignored: sbatch/makefile_demo/sims/sim0129.csv
Ignored: sbatch/makefile_demo/sims/sim0130.csv
Ignored: sbatch/makefile_demo/sims/sim0131.csv
Ignored: sbatch/makefile_demo/sims/sim0132.csv
Ignored: sbatch/makefile_demo/sims/sim0133.csv
Ignored: sbatch/makefile_demo/sims/sim0134.csv
Ignored: sbatch/makefile_demo/sims/sim0135.csv
Ignored: sbatch/makefile_demo/sims/sim0136.csv
Ignored: sbatch/makefile_demo/sims/sim0137.csv
Ignored: sbatch/makefile_demo/sims/sim0138.csv
Ignored: sbatch/makefile_demo/sims/sim0139.csv
Ignored: sbatch/makefile_demo/sims/sim0140.csv
Ignored: sbatch/makefile_demo/sims/sim0141.csv
Ignored: sbatch/makefile_demo/sims/sim0142.csv
Ignored: sbatch/makefile_demo/sims/sim0143.csv
Ignored: sbatch/makefile_demo/sims/sim0144.csv
Ignored: sbatch/makefile_demo/sims/sim0145.csv
Ignored: sbatch/makefile_demo/sims/sim0146.csv
Ignored: sbatch/makefile_demo/sims/sim0147.csv
Ignored: sbatch/makefile_demo/sims/sim0148.csv
Ignored: sbatch/makefile_demo/sims/sim0149.csv
Ignored: sbatch/makefile_demo/sims/sim0150.csv
Ignored: sbatch/makefile_demo/sims/sim0151.csv
Ignored: sbatch/makefile_demo/sims/sim0152.csv
Ignored: sbatch/makefile_demo/sims/sim0153.csv
Ignored: sbatch/makefile_demo/sims/sim0154.csv
Ignored: sbatch/makefile_demo/sims/sim0155.csv
Ignored: sbatch/makefile_demo/sims/sim0156.csv
Ignored: sbatch/makefile_demo/sims/sim0157.csv
Ignored: sbatch/makefile_demo/sims/sim0158.csv
Ignored: sbatch/makefile_demo/sims/sim0159.csv
Ignored: sbatch/makefile_demo/sims/sim0160.csv
Ignored: sbatch/makefile_demo/sims/sim0161.csv
Ignored: sbatch/makefile_demo/sims/sim0162.csv
Ignored: sbatch/makefile_demo/sims/sim0163.csv
Ignored: sbatch/makefile_demo/sims/sim0164.csv
Ignored: sbatch/makefile_demo/sims/sim0165.csv
Ignored: sbatch/makefile_demo/sims/sim0166.csv
Ignored: sbatch/makefile_demo/sims/sim0167.csv
Ignored: sbatch/makefile_demo/sims/sim0168.csv
Ignored: sbatch/makefile_demo/sims/sim0169.csv
Ignored: sbatch/makefile_demo/sims/sim0170.csv
Ignored: sbatch/makefile_demo/sims/sim0171.csv
Ignored: sbatch/makefile_demo/sims/sim0172.csv
Ignored: sbatch/makefile_demo/sims/sim0173.csv
Ignored: sbatch/makefile_demo/sims/sim0174.csv
Ignored: sbatch/makefile_demo/sims/sim0175.csv
Ignored: sbatch/makefile_demo/sims/sim0176.csv
Ignored: sbatch/makefile_demo/sims/sim0177.csv
Ignored: sbatch/makefile_demo/sims/sim0178.csv
Ignored: sbatch/makefile_demo/sims/sim0179.csv
Ignored: sbatch/makefile_demo/sims/sim0180.csv
Ignored: sbatch/makefile_demo/sims/sim0181.csv
Ignored: sbatch/makefile_demo/sims/sim0182.csv
Ignored: sbatch/makefile_demo/sims/sim0183.csv
Ignored: sbatch/makefile_demo/sims/sim0184.csv
Ignored: sbatch/makefile_demo/sims/sim0185.csv
Ignored: sbatch/makefile_demo/sims/sim0186.csv
Ignored: sbatch/makefile_demo/sims/sim0187.csv
Ignored: sbatch/makefile_demo/sims/sim0188.csv
Ignored: sbatch/makefile_demo/sims/sim0189.csv
Ignored: sbatch/makefile_demo/sims/sim0190.csv
Ignored: sbatch/makefile_demo/sims/sim0191.csv
Ignored: sbatch/makefile_demo/sims/sim0192.csv
Ignored: sbatch/makefile_demo/sims/sim0193.csv
Ignored: sbatch/makefile_demo/sims/sim0194.csv
Ignored: sbatch/makefile_demo/sims/sim0195.csv
Ignored: sbatch/makefile_demo/sims/sim0196.csv
Ignored: sbatch/makefile_demo/sims/sim0197.csv
Ignored: sbatch/makefile_demo/sims/sim0198.csv
Ignored: sbatch/makefile_demo/sims/sim0199.csv
Ignored: sbatch/makefile_demo/sims/sim0200.csv
Ignored: sbatch/makefile_demo/sims/sim0201.csv
Ignored: sbatch/makefile_demo/sims/sim0202.csv
Ignored: sbatch/makefile_demo/sims/sim0203.csv
Ignored: sbatch/makefile_demo/sims/sim0204.csv
Ignored: sbatch/makefile_demo/sims/sim0205.csv
Ignored: sbatch/makefile_demo/sims/sim0206.csv
Ignored: sbatch/makefile_demo/sims/sim0207.csv
Ignored: sbatch/makefile_demo/sims/sim0208.csv
Ignored: sbatch/makefile_demo/sims/sim0209.csv
Ignored: sbatch/makefile_demo/sims/sim0210.csv
Ignored: sbatch/makefile_demo/sims/sim0211.csv
Ignored: sbatch/makefile_demo/sims/sim0212.csv
Ignored: sbatch/makefile_demo/sims/sim0213.csv
Ignored: sbatch/makefile_demo/sims/sim0214.csv
Ignored: sbatch/makefile_demo/sims/sim0215.csv
Ignored: sbatch/makefile_demo/sims/sim0216.csv
Ignored: sbatch/makefile_demo/sims/sim0217.csv
Ignored: sbatch/makefile_demo/sims/sim0218.csv
Ignored: sbatch/makefile_demo/sims/sim0219.csv
Ignored: sbatch/makefile_demo/sims/sim0220.csv
Ignored: sbatch/makefile_demo/sims/sim0221.csv
Ignored: sbatch/makefile_demo/sims/sim0222.csv
Ignored: sbatch/makefile_demo/sims/sim0223.csv
Ignored: sbatch/makefile_demo/sims/sim0224.csv
Ignored: sbatch/makefile_demo/sims/sim0225.csv
Ignored: sbatch/makefile_demo/sims/sim0226.csv
Ignored: sbatch/makefile_demo/sims/sim0227.csv
Ignored: sbatch/makefile_demo/sims/sim0228.csv
Ignored: sbatch/makefile_demo/sims/sim0229.csv
Ignored: sbatch/makefile_demo/sims/sim0230.csv
Ignored: sbatch/makefile_demo/sims/sim0231.csv
Ignored: sbatch/makefile_demo/sims/sim0232.csv
Ignored: sbatch/makefile_demo/sims/sim0233.csv
Ignored: sbatch/makefile_demo/sims/sim0234.csv
Ignored: sbatch/makefile_demo/sims/sim0235.csv
Ignored: sbatch/makefile_demo/sims/sim0236.csv
Ignored: sbatch/makefile_demo/sims/sim0237.csv
Ignored: sbatch/makefile_demo/sims/sim0238.csv
Ignored: sbatch/makefile_demo/sims/sim0239.csv
Ignored: sbatch/makefile_demo/sims/sim0240.csv
Ignored: sbatch/makefile_demo/sims/sim0241.csv
Ignored: sbatch/makefile_demo/sims/sim0242.csv
Ignored: sbatch/makefile_demo/sims/sim0243.csv
Ignored: sbatch/makefile_demo/sims/sim0244.csv
Ignored: sbatch/makefile_demo/sims/sim0245.csv
Ignored: sbatch/makefile_demo/sims/sim0246.csv
Ignored: sbatch/makefile_demo/sims/sim0247.csv
Ignored: sbatch/makefile_demo/sims/sim0248.csv
Ignored: sbatch/makefile_demo/sims/sim0249.csv
Ignored: sbatch/makefile_demo/sims/sim0250.csv
Untracked files:
Untracked: R/demo.nmf.R
Untracked: R/nmfmu.R
Untracked: R/scd.R
Untracked: analysis/linreg_methods_demo_cache/
Unstaged changes:
Deleted: R/flashier_sparse.R
Modified: R/linreg_methods_demo_functions.R
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 repository in which changes were
made to the R Markdown (analysis/log1p_vs_pois_nmf.Rmd) and
HTML (docs/log1p_vs_pois_nmf.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 0465dcc | Peter Carbonetto | 2026-07-21 | wflow_publish("log1p_vs_pois_nmf.Rmd", verbose = T, view = F) |
| Rmd | 8bf2b0a | Peter Carbonetto | 2026-07-20 | Added a note to the log1p_vs_pois_nmf analysis. |
| html | 30f5123 | Peter Carbonetto | 2026-07-10 | Build site. |
| Rmd | b8e30a0 | Peter Carbonetto | 2026-07-10 | wflow_publish("log1p_vs_pois_nmf.Rmd", view = F, verbose = T) |
| html | 0ef40da | Peter Carbonetto | 2026-07-10 | Build site. |
| Rmd | 24d1c93 | Peter Carbonetto | 2026-07-10 | wflow_publish("log1p_vs_pois_nmf.Rmd", verbose = T, view = F) |
| html | 23da987 | Peter Carbonetto | 2026-07-10 | Modified the scatterplots in the log1p_vs_pois_nmf analysis. |
| Rmd | 407c9d4 | Peter Carbonetto | 2026-07-10 | wflow_publish("log1p_vs_pois_nmf.Rmd", verbose = T, view = F) |
| html | 2b2d8be | Peter Carbonetto | 2026-07-09 | Some improvements to the log1p_vs_pois_nmf analysis. |
| Rmd | 5216309 | Peter Carbonetto | 2026-07-09 | wflow_publish("log1p_vs_pois_nmf.Rmd", view = F, verbose = T) |
| Rmd | 4bbe955 | Peter Carbonetto | 2026-07-08 | Plots are done for the log1p_vs_pois_nmf analysis; need to add some notes and explanations. |
| Rmd | 7698c4c | Peter Carbonetto | 2026-07-08 | Added some plots to the log1p_vs_pois_nmf analysis. |
| Rmd | 3dc8b96 | Peter Carbonetto | 2026-07-07 | Added ‘right’ and ‘wrong’ model fits to log1p_vs_pois_nmf analysis. |
| html | da53752 | Peter Carbonetto | 2026-07-07 | Small tweak to the log1p_vs_pois_nmf analysis. |
| Rmd | 9d9e08f | Peter Carbonetto | 2026-07-07 | wflow_publish("log1p_vs_pois_nmf.Rmd") |
| html | 92eac40 | Peter Carbonetto | 2026-07-07 | Build site. |
| Rmd | ab435b7 | Peter Carbonetto | 2026-07-07 | wflow_publish("log1p_vs_pois_nmf.Rmd") |
Here we will examine in some detail the behaviour of Poisson NMF in data simulated from the Poisson log1p NMF model.
The idea is that the topic model and Poisson NMF are specific models that make specific modeling assumptions, and when these modeling assumptions are not correct, the model behaves in unexpected ways.
Traditional Poisson NMF (as well as the topic model, which is equivalent to Poisson NMF) assumes a model of count data \(x_{ij}\) in which the Poisson rates are a linear combination of the \(K\) non-negative factors (“topics”): \[ \begin{aligned} x_{ij} &\sim \mathrm{Pois}(\lambda_{ij}) \\ \lambda_{ij} &= \sum_{k=1}^K l_{ik} f_{jk}, \quad l_{ik}, f_{jk} \geq 0. \end{aligned} \] Here we consider the case when the Poisson rates are also controlled by a combination of \(K\) non-negative factors, but in nonlinear way: \[ \begin{aligned} x_{ij} &\sim \mathrm{Pois}(\lambda_{ij}), \\ \quad \lambda_{ij} &= e^{\mu_{ij} - 1}, \\ \quad \mu_{ij} &= \sum_{k=1}^K l_{ik} f_{jk}, \quad l_{ik}, f_{jk} \geq 0. \end{aligned} \] Then we study the behaviour of traditional Poisson NMF for data simulated under this model. What we will see is that traditional Poisson NMF produces a qualitatively different result than what we would get if we fit the true model to these data.
First load the packages needed for this analysis:
library(MCMCpack)
library(fastTopics)
library(log1pNMF)
library(ebnm)
library(flashier)
library(singlecelljamboreeR)Set the seed for reproducibility:
set.seed(1)Small function used to summarize some of the results below:
scale_cols <- function (A, b)
t(t(A) / b)Simulate data
Simulate a data set from the Poisson log1p NMF model with \(\alpha = 1\). Let’s imagine for the sake of concreteness that these count data are gene expression data from \(n\) samples or cells so that the columns of the data matrix correspond to genes. These are some parameters controlling the simulation:
n <- 100
m <- 500
m1 <- 25
k <- 4First generate the \({\bf L}\) matrix:
L <- rdirichlet(n,rep(1/(2*k),k))
colnames(L) <- paste0("k",1:k)Next generate the \({\bf F}\) matrix:
F <- matrix(0,m,k)
colnames(F) <- paste0("k",1:k)
F[,1] <- abs(rnorm(m))
changed <- rep(0,m)
for (i in 2:k) {
x <- rep(0,m)
j <- sample(m,m1)
x[j] <- abs(2*rnorm(m1))
F[,i] <- F[,1] + x
changed[j] <- TRUE
}This is a “sparse” model in the sense only a few genes have different expression levels compared to the first factor:
par(mar = c(4,4,2,1))
hist(F[,-1] - F[,1],breaks = 64,xlab = "f_jk - f_j1",ylab = "",main = "")
| Version | Author | Date |
|---|---|---|
| 2b2d8be | Peter Carbonetto | 2026-07-09 |
Now simulate the counts from the Poisson log1p NMF model with \(\alpha = 1\):
M <- tcrossprod(L,F)
X <- matrix(rpois(n*m,exp(as.vector(M)) - 1),n,m)Remove genes with no expression in the data set:
j <- which(colSums(X) > 0)
X <- X[,j]
F <- F[j,]
changed <- changed[j]Fit the “right model”:
fit_log1p <- fit_poisson_log1p_nmf(X,K = k,s = FALSE,cc = 1,loglik = "exact",
init_LL = L,init_FF = F,
control = list(verbose = FALSE))Now fit the “wrong model”:
fit0 <- init_poisson_nmf(X,F = F,L = L)
fit_pnmf <- fit_poisson_nmf(X,fit0 = fit0,method = "scd",
numiter = 40,verbose = "none",
control = list(extrapolate = FALSE))
fit_pnmf <- fit_poisson_nmf(X,fit0 = fit_pnmf,method = "scd",
numiter = 40,verbose = "none",
control = list(extrapolate = TRUE))Also fit an approximation to the “right model”, a Gaussian NMF model fit to the shifted log counts.
First compute the shifted log counts:
s <- rowSums(X)
s <- s/mean(s)
Y <- log1p(X/s)Set a lower bound on the variances:
n <- nrow(X)
x <- rpois(1e7,1/n)
s1 <- sd(log(x + 1))Now fit an NMF to the shifted log counts using flashier:
set.seed(1)
fl_nmf <- flashier_nmf(Y,k = 4,greedy_init = FALSE,var_type = 2,S = s1,
verbose = 0,maxiter = 100)
# Warning in c_nnmf(A, as.integer(k), init.mask$Wi, init.mask$Hi, init.mask$Wm, :
# Target tolerance not reached. Try a larger max.iter.
# Warning in report.maxiter.reached(verbose.lvl): Maximum number of iterations
# reached.First let’s quickly check that the estimates from the “right model” are reasonably accurate:
par(mfrow = c(1,2),mar = c(4,4,2,1))
L_log1p <- fit_log1p$LL
F_log1p <- fit_log1p$FF
d <- apply(L_log1p,2,max)
L_log1p <- scale_cols(L_log1p,d)
F_log1p <- scale_cols(F_log1p,1/d)
plot(L,L_log1p,pch = 20,cex = 0.65,xlab = "true membership",
ylab = "log1p NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
plot(F[,-1] - F[,1],F_log1p[,-1] - F_log1p[,1],pch = 20,cex = 0.65,
xlab = "true gene expression",ylab = "log1p NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
| Version | Author | Date |
|---|---|---|
| 2b2d8be | Peter Carbonetto | 2026-07-09 |
Indeed they are.
Interestingly Gaussian NMF fit to the shifted log counts is also quite accurate in the memberships, although a little less accurate in the expression levels:
par(mfrow = c(1,2),mar = c(4,4,2,1))
out <- ldf(fl_nmf,type = "i")
ks <- apply(cor(L,out$L),1,which.max)
L_fl <- out$L[,ks]
F_fl <- with(out,F[,ks] %*% diag(D[ks]))
plot(L[,-1],L_fl[,-1],pch = 20,cex = 0.65,xlab = "true membership",
ylab = "Gaussian NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
plot(F[,-1],F_fl[,-1],pch = 20,cex = 0.65,
xlab = "true gene expression",ylab = "Gaussian NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
Now let’s compare the estimates of F from the right model (log1p NMF) with the wrong model (Poisson NMF):
par(mfrow = c(1,2),mar = c(4,4,2,1))
L_pnmf <- fit_pnmf$L
F_pnmf <- fit_pnmf$F
d <- apply(L_pnmf,2,max)
L_pnmf <- scale_cols(L_pnmf,d)
F_pnmf <- scale_cols(F_pnmf,1/d)
plot(F[,-1] - F[,1],F_log1p[,-1] - F_log1p[,1],pch = 20,cex = 0.65,
xlab = "true gene expression",ylab = "log1p NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
plot(F[,-1] - F[,1],log(F_pnmf[,-1]/F_pnmf[,1]),pch = 20,cex = 0.65,
xlab = "true gene expression",ylab = "Poisson NMF estimate",
ylim = c(-1,5))
abline(a = 0,b = 1,lty = "dotted",col = "red")
| Version | Author | Date |
|---|---|---|
| 2b2d8be | Peter Carbonetto | 2026-07-09 |
It is interesting that the estimates from the wrong model are still quite accurate, although overall appear to be a little “noisier”—notice the wide spread of nonzero estimates when the true gene expression is zero.
However it is perhaps more interesting to compare the membership estimates from the two models:
par(mfrow = c(1,2),mar = c(4,4,2,1))
plot(L,L_log1p,pch = 20,cex = 0.65,xlab = "true membership",
ylab = "log1p NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
plot(L,L_pnmf,pch = 20,cex = 0.65,xlab = "true membership",
ylab = "Poisson NMF estimate")
abline(a = 0,b = 1,lty = "dotted",col = "red")
| Version | Author | Date |
|---|---|---|
| 2b2d8be | Peter Carbonetto | 2026-07-09 |
Although the Poisson NMF estimates are broadly correlated with the true memberships, the Poisson NMF estimates do not capture the “admixed” samples well (in the plots I have highlighted the “admixed” samples in which the memberships to two or more factors are greater than 0.1):
par(mfrow = c(1,2),mar = c(4,4,2,1))
plot(L[,c(2,2,3)],L[,c(3,4,4)],pch = 20,cex = 0.65,
xlab = "factor k1",ylab = "factor k2",
main = "true memberships",
cex.main = 1,font.main = 1)
i1 <- which(rowSums(L[,c(2,3)] > 0.1) > 1)
i2 <- which(rowSums(L[,c(2,4)] > 0.1) > 1)
i3 <- which(rowSums(L[,c(3,4)] > 0.1) > 1)
points(L[i1,2],L[i1,3],pch = 20,cex = 0.65,col = "red")
points(L[i2,2],L[i2,4],pch = 20,cex = 0.65,col = "red")
points(L[i3,3],L[i3,4],pch = 20,cex = 0.65,col = "red")
plot(L_pnmf[,c(2,2,3)],L_pnmf[,c(3,4,4)],pch = 20,cex = 0.65,
xlab = "factor k1",ylab = "factor k2",
main = "Poisson NMF estimates",
cex.main = 1,font.main = 1)
points(L_pnmf[i1,2],L_pnmf[i1,3],pch = 20,cex = 0.65,col = "red")
points(L_pnmf[i2,2],L_pnmf[i2,4],pch = 20,cex = 0.65,col = "red")
points(L_pnmf[i3,3],L_pnmf[i3,4],pch = 20,cex = 0.65,col = "red")
The result is that the Poisson NMF tends to act more like a clustering compared to the log1p NMF model.
sessionInfo()
# R version 4.3.3 (2024-02-29)
# Platform: aarch64-apple-darwin20 (64-bit)
# Running under: macOS 15.7.4
#
# Matrix products: default
# BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
#
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#
# time zone: America/Chicago
# tzcode source: internal
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] singlecelljamboreeR_0.1-39 flashier_1.0.59
# [3] ebnm_1.1-42 log1pNMF_0.1-6
# [5] fastTopics_0.7-46 MCMCpack_1.7-1
# [7] MASS_7.3-60.0.1 coda_0.19-4.1
#
# loaded via a namespace (and not attached):
# [1] distr_2.9.3 pbapply_1.7-2 rlang_1.1.6
# [4] magrittr_2.0.3 git2r_0.33.0 matrixStats_1.2.0
# [7] susieR_0.14.4 compiler_4.3.3 vctrs_0.6.5
# [10] reshape2_1.4.4 quantreg_5.98 quadprog_1.5-8
# [13] stringr_1.5.1 pkgconfig_2.0.3 crayon_1.5.3
# [16] fastmap_1.2.0 mcmc_0.9-8 promises_1.3.3
# [19] rmarkdown_2.29 MatrixModels_0.5-3 purrr_1.0.4
# [22] xfun_0.52 cachem_1.1.0 trust_0.1-8
# [25] jsonlite_2.0.0 progress_1.2.3 later_1.4.2
# [28] reshape_0.8.9 irlba_2.3.5.1 parallel_4.3.3
# [31] prettyunits_1.2.0 R6_2.6.1 bslib_0.9.0
# [34] stringi_1.8.7 RColorBrewer_1.1-3 SQUAREM_2021.1
# [37] jquerylib_0.1.4 Rcpp_1.1.0 knitr_1.50
# [40] httpuv_1.6.14 Matrix_1.6-5 splines_4.3.3
# [43] tidyselect_1.2.1 dichromat_2.0-0.1 yaml_2.3.10
# [46] lattice_0.22-5 tibble_3.3.0 plyr_1.8.9
# [49] evaluate_1.0.4 Rtsne_0.17 survival_3.5-8
# [52] RcppParallel_5.1.10 startupmsg_0.9.6.1 pillar_1.11.0
# [55] whisker_0.4.1 plotly_4.11.0 softImpute_1.4-3
# [58] generics_0.1.4 rprojroot_2.0.4 invgamma_1.2
# [61] truncnorm_1.0-9 hms_1.1.3 ggplot2_3.5.2
# [64] scales_1.4.0 ashr_2.2-66 gtools_3.9.5
# [67] RhpcBLASctl_0.23-42 glue_1.8.0 scatterplot3d_0.3-44
# [70] lazyeval_0.2.2 tools_4.3.3 data.table_1.17.6
# [73] SparseM_1.84-2 fs_1.6.6 cowplot_1.1.3
# [76] grid_4.3.3 tidyr_1.3.1 colorspace_2.1-0
# [79] sfsmisc_1.1-18 NNLM_0.4.4 deconvolveR_1.2-1
# [82] cli_3.6.5 Polychrome_1.5.1 workflowr_1.7.1
# [85] mixsqp_0.3-54 viridisLite_0.4.2 dplyr_1.1.4
# [88] uwot_0.2.3 gtable_0.3.6 sass_0.4.10
# [91] digest_0.6.37 ggrepel_0.9.6 htmlwidgets_1.6.4
# [94] farver_2.1.2 htmltools_0.5.8.1 lifecycle_1.0.4
# [97] httr_1.4.7