Skip to contents

Conducts longitudinal trend tests to analyze how microbial taxa abundance changes over time and across groups using linear mixed-effects models.

Usage

generate_taxa_trend_test_long(
  data.obj,
  subject.var,
  time.var = NULL,
  group.var = NULL,
  ref.level = NULL,
  adj.vars = NULL,
  feature.level,
  prev.filter = 0,
  abund.filter = 0,
  feature.dat.type = c("count", "proportion", "other"),
  ...
)

Arguments

data.obj

A MicrobiomeStat data object, which is a list containing at minimum the following components:

  • feature.tab: A matrix of feature abundances (taxa/genes as rows, samples as columns)

  • meta.dat: A data frame of sample metadata (samples as rows)

Optional components include:

  • feature.ann: A matrix/data frame of feature annotations (e.g., taxonomy)

  • tree: A phylogenetic tree object (class "phylo")

  • feature.agg.list: Pre-aggregated feature tables by taxonomy

Data objects can be created using converters like mStat_convert_phyloseq_to_data_obj or importers like mStat_import_qiime2_as_data_obj.

subject.var

Character string specifying the column name in meta.dat that uniquely identifies each subject or sample unit. Required for longitudinal and paired designs to track repeated measurements.

time.var

Character string specifying the column name in meta.dat containing the time variable. Required for longitudinal and paired analyses. Supports character/factor labels (e.g., "baseline", "week4") and numeric values. Some trend/volatility methods require numeric or coercible-to-numeric time values.

group.var

Character string specifying the column name in meta.dat containing the grouping variable (e.g., treatment, condition, phenotype). Used for between-group comparisons.

ref.level

Character string specifying the reference level for categorical group.var. If NULL, the first level alphabetically is used. Ignored for continuous variables.

adj.vars

Character vector specifying column names in meta.dat to be used as covariates for adjustment in statistical models. These variables will be included as fixed effects.

feature.level

Character vector specifying the taxonomic or annotation level(s) for analysis. Should match column names in feature.ann, such as "Phylum", "Family", "Genus", etc. Use "original" to analyze at the original feature level without aggregation.

prev.filter

Numeric value between 0 and 1. Features with prevalence (proportion of non-zero samples) below this threshold will be excluded from analysis. Default is usually 0 (no filtering).

abund.filter

Numeric value. Features with mean abundance below this threshold will be excluded from analysis. Default is usually 0 (no filtering).

feature.dat.type

Character string specifying the data type of feature.tab. One of:

  • "count": Raw count data (will be normalized)

  • "proportion": Relative abundance data (should sum to 1 per sample)

  • "other": Pre-transformed data (no transformation applied)

...

Additional arguments passed to downstream functions.

Value

A nested list structure where:

  • First level: Named by feature.level (e.g., "Phylum", "Genus")

  • Second level: Named by tested effects (e.g., "time", "group:time interaction")

    • For categorical group.var: Elements named as "Level vs Reference (Reference) [Main Effect]" and "Level vs Reference (Reference) [Interaction]"

    • For continuous group.var: Elements named by model terms (e.g., group.var, group.var:time.var)

    • time.var main effect is included when present in model output

  • Each element is a data.frame with the following columns:

    • Variable: Feature/taxon name

    • Coefficient: Effect size. For time effects, represents change per unit time. For group effects, represents log2 fold change. For interactions, represents difference in trends between groups

    • SE: Standard error of the coefficient

    • P.Value: Raw p-value from the statistical test (mixed-effects model or LinDA)

    • Adjusted.P.Value: FDR-adjusted p-value using Benjamini-Hochberg method

    • Mean.Abundance: Mean abundance of the feature across all samples

    • Prevalence: Proportion of samples where the feature is present (non-zero)

Details

This function requires time.var and always evaluates temporal trends. The fixed-effects structure is:

  • time.var when group.var = NULL

  • group.var * time.var when group.var is provided

  • optional additive adjustment terms from adj.vars

Random effects are modeled as (1 + time.var | subject.var).

Output includes time main effect (when available), group main effects, and group-time interaction terms when group.var is specified.

Examples

if (FALSE) { # \dontrun{
# Example 1
data("ecam.obj")
generate_taxa_trend_test_long(
  data.obj = ecam.obj,
  subject.var = "studyid",
  time.var = "month_num",
  group.var = "delivery",
  adj.vars = "diet",
  feature.level = c("Phylum","Class"),
  feature.dat.type = c("proportion")
)
generate_taxa_trend_test_long(
  data.obj = ecam.obj,
  subject.var = "studyid",
  time.var = "month_num",
  group.var = "delivery",
  feature.level = c("Phylum","Class"),
  feature.dat.type = c("proportion")
)
generate_taxa_trend_test_long(
  data.obj = ecam.obj,
  subject.var = "studyid",
  time.var = "month_num",
  group.var = NULL,
  feature.level = c("Phylum","Class"),
  feature.dat.type = c("proportion")
)

# Example 2
data("subset_T2D.obj")
test.list <- generate_taxa_trend_test_long(
  data.obj = subset_T2D.obj,
  subject.var = "subject_id",
  time.var = "visit_number",
  group.var = "subject_race",
  adj.vars = "sample_body_site",
  prev.filter = 0.1,
  abund.filter = 0.001,
  feature.level = c("Genus","Family"),
  feature.dat.type = c("count")
)

plot.list <- generate_taxa_trend_volcano_long(
  data.obj = subset_T2D.obj,
  group.var = "subject_race",
  time.var = "visit_number_num",
  test.list = test.list,
  feature.mt.method = "none")

} # }