Linear mixed models for data matrixSource:
Fits many linear mixed effects models for analysis of gaussian data with random effects, with parallelisation and optimisation for speed. It is suitable for longitudinal analysis of high dimensional data. Wald type 2 Chi-squared test is used to calculate p-values.
lmmSeq( modelFormula, maindata, metadata, id = NULL, offset = NULL, test.stat = c("Wald", "F", "LRT"), reduced = NULL, modelData = NULL, designMatrix = NULL, control = lmerControl(), cores = 1, removeSingles = FALSE, verbose = TRUE, returnList = FALSE, progress = FALSE, ... )
the model formula. This must be of the form
"~ ..."where the structure is assumed to be
"gene ~ ...". The formula must include a random effects term. See formula structure for random effects in
data matrix with genes in rows and samples in columns
a dataframe of sample information with variables in columns and samples in rows
Optional. Used to specify the column in metadata which contains the sample IDs to be used in repeated samples for random effects. If not specified, the function defaults to using the variable after the "|" in the random effects term in the formula.
Vector containing model offsets (default = NULL). If provided the
lmer()offset is set to
Character value specifying test statistic. Current options are "Wald" for type 2 Wald Chi square test using code derived and modified from car::Anova to improve speed for matrix tests. Or "F" for conditional F tests using Saiterthwaite's method of approximated Df. This uses lmerTest::lmer and is somewhat slower.
Optional reduced model formula. If this is chosen, a likelihood ratio test is used to calculate p-values instead of the default Wald type 2 Chi-squared test.
Optional dataframe. Default is generated by call to
expand.gridusing levels of variables in the formula. Used to calculate model predictions (estimated means & 95% CI) for plotting via modelPlot. It can therefore be used to add/remove points in modelPlot.
Optional custom design matrix generated by call to
FEformula. Used to calculate model predictions for plotting.
lmeroptimizer control (default =
number of cores to use for parallelisation. Default = 1.
whether to remove individuals with no repeated measures (default = FALSE)
Logical whether to display messaging (default = TRUE)
Logical whether to return results as a list or lmmSeq object (default = FALSE). Helpful for debugging.
Logical whether to display a progress bar
Other parameters passed to
lmerTest::lmer(). Only available if
test.stat = "F".
Returns an S4 class
lmmSeq object with results for gene-wise linear
mixed models; or a list of results if
TRUE, which is
useful for debugging. If all genes return errors from
lmer, then an error
message is shown and a character vector containing error messages for all
genes is returned.
By default, p-values for each model term are computed using Wald type 2
Chi-squared test as per
car::Anova(). The underlying code for this has been
optimised for speed. However, if a reduced model formula is specified by
reduced, then a likelihood ratio test (LRT) is performed instead
anova. This will double computation
time since two LMM have to be fitted for each gene. For LRT, models being
compared are optimised by maximum likelihood and not REML (
Two key methods are used to speed up computation above and beyond simple
parallelisation. The first is to speed up
lme4::lmer() by calling
lme4::lFormula() once at the start and then updating the
with new data. The 2nd speed up is through optimised code for repeated type 2
Wald Chi-squared tests (original code was derived from car::Anova). For
example, elements such as the hypothesis matrices are generated only once to
reduce unnecessarily repetitive computation, and the generation of p-values
from Chi-squared values is vectorised and performed at the end. F-tests using
lmerTest package have not been optimised and are therefore slower.
Parallelisation is performed using parallel::mclapply on unix/mac and parallel::parLapply on windows. Progress bars use pbmcapply::pbmclapply on unix/mac and pbapply::pblapply on windows.
id argument is usually optional. By default the
id column in the
metadata is determined as the term after the bar in the random effects term
of the model. Note that
id is not passed to
lmer. It is only really used
to remove singletons from the
maindata matrix and
metadata dataframe. The
id is also stored in the output from
lmmSeq and used by plotting function
modelPlot(). However, due to its flexible nature, in theory
allow for more than one random effect term, although this has not been tested
formally. In this case, it is probably prudent to specify a value for
data(PEAC_minimal_load) logtpm <- log2(tpm +1) lmmtest <- lmmSeq(~ Timepoint * EULAR_6m + (1 | PATID), maindata = logtpm[1:2, ], metadata = metadata, verbose = FALSE) names(attributes(lmmtest)) #>  "info" "formula" "stats" "predict" "reduced" "maindata" #>  "metadata" "modelData" "optInfo" "errors" "vars" "class"