flashmm: Fitting Linear Mixed-effects Models

In FLASHMM: Fast and Scalable Single Cell Differential Expression Analysis using Mixed-Effects Models

Fitting Linear Mixed-effects Models

Description

flashmm, a wrapper function of lmm, fits linear mixed-effects models (LMM) by sample-level data. The LMM parameters are estimated by either restricted maximum likelihood (REML) or maximum likelihood (ML) method with Fisher scoring (FS) gradient descent algorithm.

Usage

flashmm(
  Y,
  fixed,
  random,
  metadata,
  is.counts = FALSE,
  method = c("REML", "ML"),
  nBlocks = NULL,
  theta0 = NULL,
  max.iter = 50,
  epsilon = 1e-05,
  output.cov = TRUE,
  output.FIM = FALSE,
  output.RE = FALSE,
  output.SS = FALSE,
  verbose = FALSE
)

Arguments

Y	a features-by-samples matrix of responses. For single-cell RNA sequencing (scRNA-seq) data, Y is a genes-by-cells gene expression matrix, whose entries can be raw counts, log-transformed values, or normalized values. If Y is a matrix of raw counts, it will be log-transformed as log2(1 + Y) for differential expression analysis.
fixed	a one-sided model formula used to create the fixed-effects design matrix by model.matrix function.
random	a one-sided formula specifying a model with a single random-effects component or a list of one-sided formulas specifying multiple random-effects components, used to create the random-effects design matrix.
metadata	a data frame containing the variables named in fixed and random.
is.counts	a logical scalar indicating whether Y is a matrix of raw counts or processed values (e.g., log-transformed or normalized). If Y contains raw counts, it will be log-transformed as log2(1 + Y).
method	a character string, either "REML" or "ML". Defaults to "REML". If ‘REML’, the model is fit using REML; otherwise, using ML.
nBlocks	the number of blocks into which large datasets are subdivided to reduce memory usage when computing summary statistics. The default value, nBlocks = ceiling((ncol(Y) * 1e-08) * nrow(Y)), may not be adequate. If encountering the error: vector memory limit reached, you should increase the nBlocks value to avoid the issue.
theta0	a vector of initial values of the variance components, (s1, ...,sk, s_(k+1)), si = sigma_i^2, the i-th variance component. s_(k+1) = sigma^2, the variance component of model residual error.
max.iter	the maximal number of iterations for the iterative algorithm.
epsilon	positive convergence tolerance. If the absolute value of the first partial derivative of log likelihood is less than epsilon, the iterations converge.
output.cov	a logical scalar. If TRUE, output an array of the covariance matrices for the estimated coefficients, cov.beta.
output.FIM	a logical scalar. If TRUE, output an array of the Fisher information matrices (FIM) for fixed effects, FIM.beta[,,i] = X^T(Cov(Y[i,])^{-1}X for the i-th feature, and variance components, FIM.theta[,,i] = Var(\partial(logL_i)/\partial(theta)) = - E[\partial^2(logL_i)/{\partial(theta)\partial(theta^T)}] for the i-th feature.
output.RE	a logical scalar. If TRUE, output the best linear unbiased prediction (BLUP) of the random effects.
output.SS	a logical scalar. If TURE, output list(XX, XY, ZX, ZY, ZZ, Ynorm, n), a list of the summary-level data (summary statistics), defined by XX=t(X)%*%X, XY=t(Y%*%X), ZX=t(Z)%*%X, ZY=t(Z)%*%Y, ZZ=t(Z)%*%Z, Ynorm=rowSums(Y*Y), and n=nrow(X), which can also be computed using the sslmm function.
verbose	a logical scalar. If TRUE, print the number of features for which LMM fitting did not converge (abs(dlogL) > epsilon).

Value

A list containing the following components:

method	the method, either REML or ML, for estimating the LMM parameters (fixed effects and variance components).
epsilon	convergence tolerance.
dlogL	first partial derivatives of log-likelihoods for each feature.
logLik	maximum log-likelihoods for ML method or maximum log-restricted-likelihood for REML method.
niter	numbers of iterations for each feature.
coef	a matrix of estimated coefficients (fixed effects or beta), each column corresponds to a feature and each row one covariate.
se	a matrix of standard errors of the estimated coefficients.
t	a matrix of t-values for the fixed effects, equal to coef/se.
p	a matrix of two-sided p-values for the t-tests of the fixed effects.
df	degrees of freedom for the t-statistics (values).
cov	a array of covariance matrices of the estimated coefficients (fixed effects).
theta	a matrix of the estimated variance components, each column corresponds to a feature and each row one variance component. The last row is the variance component of the residual error.
se.theta	standard errors of the estimated theta.
FIM.beta	the Fisher information matrices (FIM) for fixed effects.
FIM.theta	the Fisher information matrices (FIM) for variance components.
RE	a matrix of the best linear unbiased prediction (BLUP) of random effects.
summary_data	a list of the summary-level data (summary statistics).

See Also

lmmfit, lmm

Examples

#Generate data
set.seed(2024)

#Y: gene expression matrix
n <- 1e3
m <- 10
Y <- matrix(rnorm(m*n), m, n)
rownames(Y) <- paste0("Gene", 1:nrow(Y))

#metadata: samples (sam) and treatments (trt)
metadata <- data.frame(trt = character(n), sam = character(n))
metadata$trt <- sample(c("A", "B"), n, replace = TRUE)
q <- 10
metadata$sam <- paste0("S", sample.int(q, n, replace = TRUE))

#Model formulas
fixed <- ~ 0 + trt
random <- ~ 0 + sam

#Fit LMM by flashmm with the model formulas based on the cell-level data
fit1 <- flashmm(Y, fixed, random, metadata = metadata, is.counts = FALSE)

#Fit LMM by lmmfit with the model design matrices based on the cell-level data
#design matrices
X <- model.matrix(fixed, metadata)
Z <- model.matrix(random, metadata)
d <- ncol(Z)
fit2 <- lmmfit(Y, X, Z, d = d)

identical(fit1, fit2)

#Fit LMM by lmm based on summary-level data
#Compute and store the summary-level data:
n <- nrow(X)
XX <- t(X)%*%X
XY <- t(Y%*%X)
ZX <- t(Z)%*%X
ZY <- t(Y%*%Z)
ZZ <- t(Z)%*%Z
Ynorm <- rowSums(Y*Y)
fit3 <- lmm(XX, XY, ZX, ZY, ZZ, Ynorm = Ynorm, n = n, d = d)

identical(fit2, fit3)

#Hypothesis testing
lmmtest(fit1)
lmmtest(fit1, index = 2)
lmmtest(fit1, contrast = cbind("B-A" = c(-1, 1)))

FLASHMM documentation built on April 8, 2026, 1:06 a.m.