glmmTMB建模后，emmeans分析（Tukey失效）能否获取p值？

Hey there! This is such a common pain point when working with big datasets—Tukey's honestly significant difference test can grind to a halt because it calculates the joint distribution of all pairwise comparisons, which gets computationally expensive fast. Let's break down your questions clearly:

Can you still get p-values using emmeans? #

Absolutely! You just need to swap out Tukey's correction for a more efficient method, or narrow down your comparisons. Here are your best options:

• Use lighter multiple-testing corrections
  Instead of adjust="tukey", try adjust="bonferroni" or adjust="fdr" (False Discovery Rate). These are way faster for large datasets because they don't rely on computing the complex multivariate distribution Tukey uses. Bonferroni applies a strict per-test correction, while FDR is more lenient and ideal for exploratory analyses. Example code:

  # Pairwise comparisons with Bonferroni correction
  emmeans(model, pairwise ~ your_grouping_factor, adjust = "bonferroni")
  
  # Or use FDR correction for better power
  emmeans(model, pairwise ~ your_grouping_factor, adjust = "fdr")
  

• Focus on pre-planned contrasts (not all pairwise comparisons)
  If you don't need to compare every group to every other one (most studies don't!), specify only the contrasts you care about (e.g., treatment vs. control, specific subgroup comparisons). This cuts computation time drastically. Example:

  # First generate the emmeans object
  emm <- emmeans(model, ~ your_grouping_factor)
  
  # Define a custom contrast (e.g., Group 1 vs. Group 2)
  custom_contrast <- contrast(emm, list("Group1 vs Group2" = c(1, -1, 0, 0)))
  
  # Get p-values for this specific contrast
  summary(custom_contrast, infer = c(TRUE, TRUE))
  

• Use uncorrected p-values (with extreme caution)
  If you're just exploring trends or have a small number of pre-specified tests, you can skip correction entirely. But always explicitly report that you didn't apply multiple-testing correction, as this increases the risk of Type I errors. Example:

  emmeans(model, pairwise ~ your_grouping_factor, adjust = "none")
  

Can you extract p-values directly from the glmmTMB emmeans results? #

It depends on what p-values you're targeting:

• Single mean p-values (vs. 0)：If you want to test whether each group's estimated mean differs significantly from 0, use infer = c(TRUE, TRUE) in the summary() call. This will return p-values alongside confidence intervals for each mean:

  emm <- emmeans(model, ~ your_grouping_factor)
  summary(emm, infer = c(TRUE, TRUE))
  

• Group comparison p-values：These aren't included in the basic emmeans output—you need to run a pairwise comparison or custom contrast first (like the examples above). Once you do that, you can extract p-values by converting the results to a data frame:

  # Run pairwise comparisons with FDR correction
  pairwise_results <- emmeans(model, pairwise ~ your_grouping_factor, adjust = "fdr")
  
  # Convert contrasts to a data frame and extract p-values
  p_values_df <- as.data.frame(pairwise_results$contrasts)
  # Access the p-value column directly
  p_values_df$p.value
  

A quick pro tip: If your dataset is extremely large, double-check if your model can be simplified (e.g., removing unnecessary random effects) before running post-hoc tests—this speeds up both model fitting and emmeans computations.

内容的提问来源于stack exchange，提问作者Klara Albert