Linear model strategy, post hoc tests, and modeling specific taxa in MicrobiomeStat

[lauraflows]

Hello,

Thank you for developing MicrobiomeStat — it's a very helpful and well-designed tool!
I'm new to applying linear models to microbiome data and would appreciate your advice on a few analysis questions.

My experimental design is relatively simple:

4 groups (n = 6 per group): 1 control and 3 treatment conditions

2 time points: T0 and 20 weeks later (T20)

Here are my main questions:

• When using generate_alpha_test_pair, the complex model often fails with the message:

Complex model failed. Trying a simpler model...
boundary (singular) fit: see help('isSingular')

In such cases, is it statistically valid to rely on the simpler linear model fallback?

• I'm interested in between-group differences at T0, within-group changes over time, and differences between treatments at T20 — ideally including multiple comparisons.
  What kind of post hoc comparisons would you recommend in this context? For example, would you use something like Tukey HSD, or is there a more appropriate approach for microbiome data?

• Finally, I’d like to analyze specific families (e.g., changes in the relative abundance of a single taxon).
  What modeling strategy would you suggest? Should I extract those features manually and analyze them separately (e.g., using CLR-transformed or raw relative abundance values)?

Thank you again for your time — I really appreciate both the tool and any guidance you can share.

Best regards,
Laura

[cafferychen777]

Hello Laura,

Thank you for using MicrobiomeStat and for your thoughtful questions! I'm happy to address your questions about linear model strategies, post hoc tests, and modeling specific taxa.

Regarding the Complex Model Failure

When you see the message "Complex model failed. Trying a simpler model..." and "boundary (singular) fit", this typically indicates that your data structure doesn't support the complex random effects structure. In the generate_alpha_test_pair function, we first attempt a complex model with random slopes and intercepts:

# Complex model (with random slopes and intercepts)
formula_part <- paste(index, "~", group.var, "*", time.var, " + (1 +", time.var, "|", subject.var, ")")

When this model fails, we fall back to a simpler model with only random intercepts:

# Simple model (with only random intercepts)
formula_part <- paste(index, "~", group.var, "*", time.var, " + (1|", subject.var, ")")

From a statistical perspective, using this simpler model is completely valid. Singular fits usually indicate that the data is insufficient to estimate all random effect parameters, and the simplified model addresses this by reducing the number of parameters to estimate. This is a common strategy when dealing with longitudinal microbiome data where sample sizes are often limited.

Recommendations for Post Hoc Comparisons

For your experimental design (4 groups, 2 time points), I recommend the following strategy:

1. Between-group differences at T0: Use generate_alpha_test_single or generate_taxa_test_single with a subset of data containing only the T0 time point.

2. Within-group changes over time: Use generate_alpha_change_test_pair or generate_taxa_change_test_pair, running it separately for each group.

3. Between-group differences at T20: Similar to point 1, but with the T20 data subset.

4. Multiple comparisons: Functions in MicrobiomeStat already implement appropriate multiple comparison corrections. For between-group comparisons, when the number of groups is >2, we use ANOVA to test for overall differences, which provides built-in protection against multiple comparisons. If you need more detailed pairwise comparisons, consider using the emmeans package for Tukey HSD or other post hoc tests:

# Example code (to be run manually after MicrobiomeStat analysis)
library(emmeans)
# Assuming lme.model is the model you obtained from MicrobiomeStat
emm <- emmeans(lme.model, specs = ~ group.var | time.var)
pairs(emm, adjust = "tukey")

Recommendations for Analyzing Specific Taxa

For analyzing specific taxa (e.g., changes in the relative abundance of a particular family), I recommend:

1. Use MicrobiomeStat's taxa functions: Functions like generate_taxa_test_pair and generate_taxa_change_test_pair are designed for analyzing specific taxa. These functions automatically handle multiple comparison issues and provide appropriate visualizations.

2. Data transformation: For relative abundance data, we default to using TSS (Total Sum Scaling) normalization for count data, which is a common method for handling compositional data. You don't need to manually extract and transform the data.

3. Example code:

# Analyzing specific taxa
taxa_results <- generate_taxa_test_pair(
  data.obj = your_data_obj,
  subject.var = "subject_id",
  time.var = "time_point",
  change.base = "T0",  # Specify the reference time point (optional)
  group.var = "treatment_group",
  adj.vars = NULL,  # Add covariates if needed
  feature.level = c("Family"),  # Or other taxonomic level
  prev.filter = 0.1,  # Filter low prevalence taxa
  abund.filter = 0.001,  # Filter low abundance taxa
  feature.dat.type = "count"  # Data type: count, proportion, or other
)

# Visualize results
volcano_plot <- generate_taxa_volcano_single(
  data.obj = your_data_obj,
  group.var = "treatment_group",
  test.list = taxa_results,
  feature.sig.level = 0.05,
  feature.mt.method = "fdr"  # Multiple testing correction method
)

If you're interested in analyzing a single specific taxon, you can first use filtering functions to create a data object containing only the taxa of interest, then apply the analyses above.

I hope this information helps! If you have any follow-up questions or need more specific code examples, please don't hesitate to ask.

Best regards,

Chen

[lauraflows]

Hi Chen,

Thanks again for your detailed response and the helpful guidance!

I just wanted to double-check that I’ve correctly understood your suggested strategy.

• I tried fitting a mixed model using lmer with only a random intercept for subject ((1 | subject)), but I’m still getting a singular fit. From your explanation, I gather that this often reflects limitations in the data structure, and that simplifying the model is an appropriate response. That said, to be sure:
  In cases like mine, where even the simpler model is singular, would you recommend splitting the data (e.g., analyzing groups separately at T0 and T20, and changes over time within each group using lm)? Or would you still go ahead with the full dataset in a single model, despite the singularity, and use emmeans for post hoc contrasts?

And just to clarify: is the reason for splitting mainly due to limited sample size and lack of covariates, making it difficult to support a valid random effects structure?

• Also, I’d appreciate your thoughts on how you would approach beta diversity in this same context:
  Would you prefer testing all time points and groups together (e.g., with PERMANOVA-style tests), or do you think splitting by time or group is better, as with alpha?
  And if global effects are found, how would you recommend handling post hoc comparisons for beta diversity?

I just want to make sure I’m interpreting your suggestions correctly and applying the most reliable approach given the paired design and limited sample size.

Thanks again for your support!

Best regards,

Laura

[cafferychen777]

Hello Laura,

Thank you for your thoughtful questions about MicrobiomeStat and for sharing your experimental design details. I'm happy to provide some guidance on linear model strategies, post hoc tests, and approaches for analyzing specific taxa in your microbiome data.

Regarding Singular Fit in Mixed Models

When you encounter messages like "Complex model failed. Trying a simpler model..." and "boundary (singular) fit", this indicates that your data structure cannot support the complex random effects structure initially attempted. This is quite common in microbiome studies with limited sample sizes.

In generate_alpha_test_pair, we first try a complex model with random slopes and intercepts:

# Complex model (with random slopes and intercepts)
formula_part <- paste(index, "~", group.var, "*", time.var, " + (1 +", time.var, "|", subject.var, ")")

When this fails, we fall back to a simpler model with only random intercepts:

# Simple model (with only random intercepts)
formula_part <- paste(index, "~", group.var, "*", time.var, " + (1|", subject.var, ")")

If even the simpler model gives a singular fit warning, it suggests that your data might be too limited to support even this reduced random effects structure. In such cases, you have several options:

1. Continue with the singular fit model: Despite the warning, the model can still provide valid fixed effects estimates, though you should interpret random effects with caution.

2. Split your analysis by time point: For your specific design (4 groups, 2 time points), you could:

   • Use generate_alpha_test_single or generate_taxa_test_single for between-group comparisons at T0 and T20 separately
   • Use generate_alpha_change_test_pair or generate_taxa_change_test_pair for within-group changes over time

3. Use a completely fixed effects approach: If you're primarily interested in group comparisons rather than subject-specific effects, you could consider analyzing your data with standard linear models rather than mixed models.

Recommended Post Hoc Testing Strategy

For your experimental design with 4 groups and 2 time points, I recommend:

1. For between-group differences at T0:

   t0_results <- generate_alpha_test_single(
     data.obj = your_data_obj,
     alpha.name = c("shannon", "observed_species", "chao1"), # Specify the required alpha diversity indices
     time.var = "time_point",
     t.level = "T0",
     group.var = "treatment_group",
     adj.vars = NULL # Specify adjustment variables here if needed
   )

2. For between-group differences at T20:

   t20_results <- generate_alpha_test_single(
     data.obj = your_data_obj,
     alpha.name = c("shannon", "observed_species", "chao1"), # Specify the required alpha diversity indices
     time.var = "time_point",
     t.level = "T20",
     group.var = "treatment_group",
     adj.vars = NULL # Specify adjustment variables here if needed
   )

3. For within-group changes over time:

   # For each group separately
   for (group in unique(your_data_obj$meta.dat$treatment_group)) {
     # Create a subset for this group
     group_condition <- paste("treatment_group == '", group, "'", sep="")
     group_data <- mStat_subset_data(your_data_obj, condition = group_condition)
     
     # Analyze changes over time within this group
     change_results <- generate_alpha_change_test_pair(
       data.obj = group_data,
       subject.var = "subject_id",
       time.var = "time_point",
       change.base = "T0"
     )
   }

For detailed pairwise comparisons and post hoc tests, you can use the emmeans package with the models generated by MicrobiomeStat:

library(emmeans)
# For between-group comparisons at T20
t20_model <- lm(alpha_diversity ~ treatment_group, data = your_alpha_data_t20)
emm_t20 <- emmeans(t20_model, specs = ~ treatment_group)
pairs(emm_t20, adjust = "tukey")  # Tukey's HSD adjustment

# For longitudinal comparisons
pair_model <- lmer(alpha_diversity ~ treatment_group * time_point + (1|subject_id), data = your_alpha_data)
emm_pair <- emmeans(pair_model, specs = ~ treatment_group | time_point)
pairs(emm_pair, adjust = "tukey")  # Between groups at each time
emm_time <- emmeans(pair_model, specs = ~ time_point | treatment_group)
pairs(emm_time)  # Within groups across time

Beta Diversity Analysis Strategy

For beta diversity in your experimental design:

1. Global test first: For longitudinal beta diversity analysis, we can use two approaches:

   a. Analyzing changes between time points:

   beta_change_results <- generate_beta_change_test_pair(
     data.obj = your_data_obj,
     dist.obj = NULL,  # Will be calculated if NULL
     time.var = "time_point",
     subject.var = "subject_id",
     group.var = "treatment_group",
     adj.vars = NULL,
     change.base = "T0",
     dist.name = c('BC', 'Jaccard') # Specify the required distance metrics
   )

   b. Analyzing each time point separately:

   # Analysis for the T0 time point
   beta_t0_results <- generate_beta_test_single(
     data.obj = your_data_obj,
     dist.obj = NULL,
     time.var = "time_point",
     t.level = "T0",
     group.var = "treatment_group",
     adj.vars = NULL,
     dist.name = c('BC', 'Jaccard')
   )
   
   # Analysis for the T20 time point
   beta_t20_results <- generate_beta_test_single(
     data.obj = your_data_obj,
     dist.obj = NULL,
     time.var = "time_point",
     t.level = "T20",
     group.var = "treatment_group",
     adj.vars = NULL,
     dist.name = c('BC', 'Jaccard')
   )

2. Stratified analyses: If significant effects are found, follow up with:

   • Between-group tests at each time point using generate_beta_test_single
   • Within-group tests across time using generate_beta_change_test_pair

3. Pairwise comparisons: For post-hoc testing in beta diversity, you can use the pairwise parameter of the adonis2 function:

   # First calculate the distance matrix
   dist_obj <- mStat_calculate_beta_diversity(
     data.obj = your_data_obj,
     dist.name = c('BC', 'Jaccard')
   )
   
   # Extract data for the T20 time point
   t20_data <- subset(your_data_obj$meta.dat, time_point == "T20")
   
   # Use the adonis2 function from the vegan package for pairwise comparisons
   library(vegan)
   adonis2_result <- adonis2(dist_obj$BC ~ treatment_group, data = t20_data, pairwise = TRUE)
   print(adonis2_result)

Analyzing Specific Taxa

For analyzing specific families or taxa:

1. Using MicrobiomeStat's built-in functions:

   # For differential abundance testing of all taxa
   taxa_results <- generate_taxa_test_pair(
     data.obj = your_data_obj,
     subject.var = "subject_id",
     time.var = "time_point",
     group.var = "treatment_group",
     feature.level = c("Family"),  # Or other taxonomic level
     prev.filter = 0.1,
     abund.filter = 0.001
   )
   
   # For visualizing results
   volcano_plot <- generate_taxa_volcano_single(
     data.obj = your_data_obj,
     group.var = "treatment_group",
     test.list = taxa_results,
     feature.sig.level = 0.05,
     feature.mt.method = "fdr"
   )

2. For specific taxa of interest:
   If you want to focus on just one or a few specific taxa, you can:

   a. Extract the specific taxon data:

   # Assuming you're interested in a specific family
   family_of_interest <- "Lachnospiraceae"
   
   # Extract the family-level data
   if (is.null(data.obj$feature.agg.list[["Family"]])) {
     data.obj <- mStat_aggregate_by_taxonomy(data.obj, feature.level = "Family")
   }
   family_data <- data.obj$feature.agg.list[["Family"]]
   
   # Extract the specific family
   if (family_of_interest %in% rownames(family_data)) {
     specific_family_abundance <- family_data[family_of_interest, ]
     
     # Create a data frame for analysis
     analysis_df <- data.frame(
       sample = names(specific_family_abundance),
       abundance = as.numeric(specific_family_abundance)
     ) %>%
     inner_join(data.obj$meta.dat %>% rownames_to_column("sample"), by = "sample")
   }

   b. Analyze using standard statistical methods:

   # For longitudinal analysis
   library(lme4)
   model <- lmer(abundance ~ treatment_group * time_point + (1|subject_id), data = analysis_df)
   anova(model)
   
   # Post-hoc tests
   library(emmeans)
   emm <- emmeans(model, specs = ~ treatment_group | time_point)
   pairs(emm, adjust = "tukey")

Conclusion

Given your experimental design (4 groups, 2 time points, n=6 per group), I recommend:

1. For singular fit issues: Split your analysis by time point and/or group to reduce model complexity, which is appropriate given your limited sample size.

2. For post hoc comparisons: Use the emmeans package with Tukey's HSD adjustment for multiple comparisons after identifying significant overall effects.

3. For specific taxa: Use MicrobiomeStat's functions for comprehensive analysis, or extract specific taxa of interest for targeted analysis with standard statistical methods.

I hope this guidance helps with your analysis! Let me know if you have any further questions or need additional clarification.

Best regards,
Chen