Skip to contents

Compute linear mixed-effects models (LMM) on longitudinal microbiome data

Source: R/run_lmm_model.R
run_lmm_model.Rd

This function takes an animalcules-formatted MultiAssayExperiment and runs an independent LMM model for each taxon. The model predicts taxon log CPM abundance as a product of fixed-effects covariates with a random effect, usually the unit on which repeated measurements were taken. Note, the 'broom', 'lmerTest', and 'broom.mixed' packages are required to use this function; they can be downloaded from CRAN.

Usage

run_lmm_model(
  dat,
  taxon_level = "genus",
  unit_var,
  fixed_cov,
  p_adj_method = "fdr",
  plot_out = FALSE,
  plotsave_loc = ".",
  plot_terms = NULL,
  ...
)

Arguments

dat

A MultiAssayExperiment object specially formatted as an animalcules output.

taxon_level

Character string, default is "genus".

unit_var

Character string giving the variable containing the identifiers for the unit on which multiple measurements were conducted, e.g. subjects. Default is NULL; must be supplied if paired = FALSE.

fixed_cov

A character vector naming covariates to be tested.

p_adj_method

A character string specifying the correction method. Can be abbreviated. See details. Default is "fdr".

plot_out

Logical indicating whether plots should be output alongside the model results. Default is FALSE.

plotsave_loc

A character string giving the folder path to save plot outputs. This defaults to the current working directory.

plot_terms

Character vector. Which terms should be examined in the plot output? Can overlap with the fixed_cov inputs.

...

Further arguments passed to ggsave for plot creation.

Value

A data.frame of modeling results.

Details

P-values are adjusted for the model coefficients within each taxon. The following methods are permitted: c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none")

Examples

dat <- system.file("extdata/MAE.RDS", package = "LegATo") |>
  readRDS() |>
  filter_MAE()
#> The overall range of relative abundance counts between samples is (13218, 3016276) 
#> Number of OTUs that exhibit a relative abundance >3% in at least 5% of the total samples: 17/1690
out <- run_lmm_model(dat, taxon_level = "genus", unit_var = "Subject",
                     fixed_cov = c("HIVStatus", "timepoint"))
#> Warning: Model may not have converged with 1 eigenvalue close to zero: 8.6e-09
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
#> boundary (singular) fit: see help('isSingular')
head(out)
#>   effect Coefficient Coefficient Estimate Standard Error Statistic        df
#> 1  fixed (Intercept)         155224.18209   40395.201882  3.842639 118.81582
#> 2  fixed (Intercept)             40.10844      23.584492  1.700628  26.32407
#> 3  fixed (Intercept)             15.74961       6.753684  2.332003  34.56949
#> 4  fixed (Intercept)              9.80022       4.584216  2.137818 164.00000
#> 5  fixed (Intercept)              4.22703       2.145746  1.969958  90.57703
#> 6  fixed (Intercept)              5.52326       4.591306  1.202982 164.00000
#>   Unadj p-value               Taxon Adj p-value
#> 1  0.0001969403               Other 0.003347984
#> 2  0.1007992523            Belnapia 0.244798184
#> 3  0.0256565784         Actinotalea 0.144554929
#> 4  0.0340129244          Leptothrix 0.144554929
#> 5  0.0518978791 Limosilactobacillus 0.176452789
#> 6  0.2307169043      Staphylococcus 0.392218737