How to bootstrap non-negative matrix factorization results?

Question

I have RNA-seq data from 9 samples and around 15,000 genes. I know that these 9 samples consist of varying proportions of two cell types, each with their own expression profile. I am using non-negative matrix factorization with a rank of 2 to determine the cell-type ratios in each sample. In other words, I use NMF to decompose my 9 x 15,000 matrix into a 2 x 15,000 matrix and a 2 x 9 matrix.

I want to obtain confidence intervals for my estimates of the 2 x 9 matrix. From what I understand about bootstrapping, I could do this by resampling from my 9 observations a bunch of times and keep track of each estimate for each observation.

My question is, since NMF doesn't have any notion of observations vs. features, could I also estimate the variance by resampling genes (features) rather than observations?

Answer 1

It turns out the answer is yes.
I read Efron and Tibshirani's textbook on the bootstrap and they gave an example of bootstrapping PCA applied to standardized test data where a few students answered many questions. In this example they resampled questions (equivalent to genes in my example). Since NMF is a special case of the more general multinomial PCA, it seems to follow that this approach should also work with NMF.

Answer 2

NNMF is unique to within a permutation and a scaling, per Laurberg, 2007, 14th IEEE/SP Workshop on Statistical Signal Processing, August 2007. Related, Donoho, Stodden "When does NNMF give a correct decomposition into parts?", Stanford, 2003.

Answer 3

Executive summary: "yes".

IIUC, your observed data is 9 rows by 15000 columns, and you have considered sampling-with-replacement from the rows to generate 9x15000 for variance estimation. As you probably suspected, this has the problem that a large number of the rows are likely to be duplicates of each other in each sample.

IIUC, you are proposing instead sampling columns. This sounds reasonable to me. Since there are 15000 columns, you will expect relatively few sampled columns to have a duplicate. Since you are only estimating a 2-rank matrix, I would sample less than 15000 columns, say only 1000 columns, and I would do the 1000 column sampling without replacement.