lmmfit | R Documentation |
lmmfit, 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.
lmmfit(
Y,
X,
Z,
d = ncol(Z),
theta0 = NULL,
nBlocks = ceiling((ncol(Y) * 1e-08) * nrow(Y)),
method = c("REML", "ML"),
max.iter = 50,
epsilon = 1e-05,
output.cov = TRUE,
output.RE = FALSE
)
Y |
A features-by-samples matrix of responses (genes-by-cells matrix of gene expressions for scRNA-seq). |
X |
A design matrix for fixed effects, with rows corresponding to the columns of Y. |
Z |
A design matrix for random effects, with rows corresponding to the columns of Y. Z = [Z1, ..., Zk], and Zi, i=1,...,k, is the design matrix for the i-th type random factor. |
d |
= (d1,...,dk), where di = ncol(Zi), number of columns in Zi. sum(d) = ncol(Z), number of columns in Z. For the model with only one random factor, d = ncol(Z). |
theta0 |
A vector of initial values of the variance components, (s1, ...,sk, s_(k+1)), si = sigma_i^2, the variance component of the i-th type random effects. s_(k+1) = sigma^2, the variance component of model residual error. |
nBlocks |
Number of the blocks, which a big data is subdivided into, used for reducing the storage in computing the summary statistics that are computed from a block of data. The default value may not be adequate. If encountering the error: vector memory limit reached, you should increase the nBlocks value to avoid the issue. |
method |
The REML with Fisher scoring (FS) iterative algorithm, REML-FS. |
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 |
If TRUE, output the covariance matrices for the estimated coefficients, which are needed for testing contrasts. |
output.RE |
If TRUE, output the best linear unbiased prediction (BLUP) of the random effects. |
A list containing the following components:
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), 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. |
df |
Degrees of freedom for the t-statistics (values). |
p |
A matrix of two-sided p-values for the t-tests of the fixed effects. |
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. |
RE |
A matrix of the best linear unbiased prediction (BLUP) of random effects. |
lmm
#Generate data: X, Y, and Z.
set.seed(2024)
n <- 1e3
m <- 10
Y <- matrix(rnorm(m*n), m, n)
rownames(Y) <- paste0("Gene", 1:nrow(Y))
trt <- sample(c("A", "B"), n, replace = TRUE)
X <- model.matrix(~ 0 + trt)
q <- 20
sam <- rep(NA, n)
sam[trt == "A"] <- paste0("A", sample.int(q/2, sum(trt == "A"), replace = TRUE))
sam[trt == "B"] <- paste0("B", sample.int(q/2, sum(trt == "B"), replace = TRUE))
Z <- model.matrix(~ 0 + sam)
d <- ncol(Z)
#Fit LMM by the cell-level data
fit <- lmmfit(Y, X, Z, d = d)
str(fit)
#Fit LMM by 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)
fitss <- lmm(XX, XY, ZX, ZY, ZZ, Ynorm = Ynorm, n = n, d = d)
identical(fit, fitss)
#Hypothesis testing
lmmtest(fit)
lmmtest(fit, index = 2)
lmmtest(fit, contrast = cbind("B-A" = c(-1, 1)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.