Fit Poisson GLM-PCA Model to Count Data

Description

Fit a Poisson GLM-PCA model by maximum-likelihood.

Usage

fit_glmpca_pois(
  Y,
  K,
  fit0 = init_glmpca_pois(Y, K),
  verbose = TRUE,
  control = list()
)

fit_glmpca_pois_control_default()

init_glmpca_pois(
  Y,
  K,
  U,
  V,
  X = numeric(0),
  Z = numeric(0),
  B = numeric(0),
  W = numeric(0),
  fixed_b_cols = numeric(0),
  fixed_w_cols = numeric(0),
  col_size_factor = TRUE,
  row_intercept = TRUE
)

Arguments

Details

In generalized principal component analysis (GLM-PCA) based on a Poisson likelihood, the counts y_{ij} stored in an n \times m matrix Y are modeled as

y_{ij} \sim Pois(\lambda_{ij}),

in which the logarithm of each rate parameter \lambda_{ij} is defined as a linear combination of rank-K matrices to be estimated from the data:

\log \lambda_{ij} = (UDV')_{ij},

where U and V are orthogonal matrices of dimension n \times K and m \times K, respectively, and D is a diagonal K \times K matrix in which the entries along its diagonal are positive and decreasing. K is a tuning parameter specifying the rank of the matrix factorization. This is the same as the low-rank matrix decomposition underlying PCA (that is, the singular value decomposition), but because we are not using a linear (Gaussian) model, this is called “generalized PCA” or “GLM PCA”.

To allow for additional components that may be fixed, fit_glmpca_pois can also fit the more general model

\log \lambda_{ij} = (UDV' + XB' + WZ')_{ij},

in which X, Z are fixed matrices of dimension n \times n_x and m \times n_z, respectively, and B, W are matrices of dimension m \times n_x and n \times n_z to be estimated from the data.

fit_glmpca_pois computes maximum-likelihood estimates (MLEs) of U, V, D, B and W satistifying the orthogonality constraints for U and V and the additional constraints on D that the entries are positive and decreasing. This is accomplished by iteratively fitting a series of Poisson GLMs, where each of these individual Poissons GLMs is fitted using a fast “cyclic co-ordinate descent” (CCD) algorithm.

The control argument is a list in which any of the following named components will override the default optimization algorithm settings (as they are defined by fit_glmpca_pois_control_default). Additional control arguments not listed here can be used to control the behaviour of fpiter or daarem; see the help accompanying these functions for details.

use_daarem

If use_daarem = TRUE, the updates are accelerated using DAAREM; see daarem for details.

tol

This is the value of the “tol” control argument for fpiter or daarem that determines when to stop the optimization. In brief, the optimization stops when the change in the estimates or in the log-likelihood between two successive updates is less than “tol”.

maxiter

This is the value of the “maxiter” control argument for fpiter or daarem. In brief, it sets the upper limit on the number of CCD updates.

convtype

This is the value of the “convtype” control argument for daarem. It determines whether the stopping criterion is based on the change in the estimates or the change in the log-likelihood between two successive updates.

mon.tol

This is the value of the “mon.tol” control argument for daarem. This setting determines to what extent the monotonicity condition can be violated.

training_frac

Fraction of the columns of input data Y to fit initial model on. If set to 1 (default), the model is fit by optimizing the parameters on the entire dataset. If set between 0 and 1, the model is optimized by first fitting a model on a randomly selected fraction of the columns of Y, and then projecting the remaining columns of Y onto the solution. Setting this to a smaller value will increase speed but decrease accuracy.

num_projection_ccd_iter

Number of co-ordinate descent updates be made to elements of V if and when a subset of Y is projected onto U. Only used if training_frac is less than 1.

num_ccd_iter

Number of co-ordinate descent updates to be made to parameters at each iteration of the algorithm.

line_search

If line_search = TRUE, a backtracking line search is performed at each iteration of CCD to guarantee improvement in the objective (the log-likelihood).

ls_alpha

alpha parameter for backtracking line search. (Should be a number between 0 and 0.5, typically a number near zero.)

ls_beta

beta parameter for backtracking line search controlling the rate at which the step size is decreased. (Should be a number between 0 and 0.5.)

calc_deriv

If calc_deriv = TRUE, the maximum gradient of U and V is calculated and stored after each update. This may be useful for assessing convergence of the optimization, though increases overhead.

calc_max_diff

If calc_max_diff = TRUE, the largest change in U and V after each update is calculated and stored. This may be useful for monitoring progress of the optimization algorithm.

orthonormalize

If orthonormalize = TRUE, the matrices U and V are made to be orthogonal after each update step. This improves the speed of convergence without the DAAREM acceleration; however, should not be used when use_daarem = TRUE.

You may use function set_fastglmpca_threads to adjust the number of threads used in performing the updates.

Value

An object capturing the state of the model fit. It contains estimates of U, V and D (stored as matrices U, V and a vector of diagonal entries d, analogous to the svd return value); the other parameters (X, B, Z, W); the log-likelihood achieved (loglik); information about which columns of B and W are fixed (fixed_b_cols, fixed_w_cols); and a data frame progress storing information about the algorithm's progress after each update.

References

Townes, F. W., Hicks, S. C., Aryee, M. J. and Irizarry, R. A. (2019). Feature selection and dimension reduction for single-cell RNA-Seq based on a multinomial model. Genome Biology 20, 295. doi:10.1186/s13059-019-1861-6

Collins, M., Dasgupta, S. and Schapire, R. E. (2002). A generalization of principal components analysis to the exponential family. In Advances in Neural Information Processing Systems 14.

See Also

fit_glmpca_pois

Examples

set.seed(1)
n <- 200
p <- 100
K <- 3
dat  <- generate_glmpca_data_pois(n,p,K)
fit0 <- init_glmpca_pois(dat$Y,K)
fit  <- fit_glmpca_pois(dat$Y,fit0 = fit0)

Get Fitted Values for GLM-PCA Model Fit

Description

fitted method for the “glmpca_pois_fit” class.

Usage

## S3 method for class 'glmpca_pois_fit'
fitted(object, ...)

Arguments

Value

An n x p matrix of fitted means. Calculated as

exp(UDV')

using the fit object.

Generate Data from a GLM-PCA Model

Description

Generate data from a GLM-PCA model with a specified rank.

Usage

generate_glmpca_data_pois(n, p, K, link = c("log", "log1p"))

Arguments

Details

This function assumes that each column of the data is generated from a multinomial distribution. Let

Y_j

denote column j of the generated data matrix. First, we set

sum(Y_j)

equal to a value generated from a

Uniform(50, 5000)

distribution. Then, we generate

L

and

F

from mixture distributions, and calculate

H = exp(L'F)

. Then, we generate the individual elements of

Y_j

from a multinomial model where the probability for each individual element is just

H_j

normalized.

Value

list with the following components

• LL - loadings of underlying mean structure. A K x n matrix

• FF - factors of underlying mean structure. A K x p matrix

• Y - n x p matrix of generated data.

Examples

set.seed(1)
sim_data <- generate_glmpca_data_pois(1000, 500, 1)

Mixture of 10 FACS-purified PBMC Single-Cell RNA-seq data

Description

These data are a selection of the reference transcriptome profiles generated via single-cell RNA sequencing (RNA-seq) of 10 bead-enriched subpopulations of PBMCs (Donor A), described in Zheng et al (2017). The data are unique molecular identifier (UMI) counts for 16,791 genes in 3,774 cells. (Genes with no expression in any of the cells were removed.) Since the majority of the UMI counts are zero, they are efficiently stored as a 16,791 x 3774 sparse matrix. These data are used in the vignette illustrating how ‘fastglmpca’ can be used to analyze single-cell RNA-seq data. Data for a separate set of 1,000 cells is provided as a “test set” to evaluate out-of-sample predictions.

Format

pbmc_facs is a list with the following elements:

counts

16,791 x 3,774 sparse matrix of UMI counts, with rows corresponding to genes and columns corresponding to cells (samples). It is an object of class "dgCMatrix").

counts_test

UMI counts for an additional test set of 100 cells.

samples

Data frame containing information about the samples, including cell barcode and source FACS population (“celltype” and “facs_subpop”).

samples_test

Sample information for the additional test set of 100 cells.

genes

Data frame containing information and the genes, including gene symbol and Ensembl identifier.

fit

GLM-PCA model that was fit to the UMI count data in the vignette.

Source

https://www.10xgenomics.com/resources/datasets

References

G. X. Y. Zheng et al (2017). Massively parallel digital transcriptional profiling of single cells. Nature Communications 8, 14049. doi:10.1038/ncomms14049

Examples

library(Matrix)
data(pbmc_facs)
cat(sprintf("Number of genes: %d\n",nrow(pbmc_facs$counts)))
cat(sprintf("Number of cells: %d\n",ncol(pbmc_facs$counts)))
cat(sprintf("Proportion of counts that are non-zero: %0.1f%%.\n",
            100*mean(pbmc_facs$counts > 0)))

Set up Multithreading for fastglmpca

Description

Initialize RcppParallel multithreading using a pre-specified number of threads, or using the default number of threads when n is not specified or is NA.

Usage

set_fastglmpca_threads(n)

Arguments

Value

The number of threads to be used.

Summarize GLM-PCA Model Fit

Description

summary method for objects of class “glmpcan_fit_pois”.

Usage

## S3 method for class 'glmpca_pois_fit'
summary(object, ...)

## S3 method for class 'summary.glmpca_pois_fit'
print(x, ...)

Arguments

Value

summary returns a vector of basic statistics summarizing the model fit.