使用lme&glht做单因素混合效应ANOVA的Tukey检验时p值异常

I’ve run into similar headaches when trying to replicate classic textbook ANOVA results with modern mixed-effects tools, so I totally get where you’re coming from! The fact that your F-statistic matches Hays’ 43.41 means your core mixed-effects model is correctly specified for the fixed and random effects—so the issue almost certainly lies in how the Tukey test is calculating standard errors or degrees of freedom compared to the traditional repeated-measures ANOVA approach in the textbook.

Here’s what to check and fix:

1. Match the Traditional ANOVA’s Variance Assumption #

Traditional repeated-measures ANOVA assumes compound symmetry (equal variances across treatments and equal correlations between all pairs of repeated measures). By default, lme doesn’t enforce this structure, which leads to mismatched standard errors compared to the textbook’s approach.

Align your model with Hays’ assumptions by adding a compound symmetry correlation structure:

library(nlme)
# Replace placeholders with your actual variable names
model <- lme(dv ~ treatment, 
             random = ~1 | subject,
             correlation = corCompSymm(form = ~1 | subject), # Enforces compound symmetry
             data = your_data_frame)

This ensures the variance-covariance structure matches what’s used in the classic ANOVA, which directly impacts the standard errors driving the Tukey test.

2. Use emmeans Instead of glht for Mixed-Model Friendly Tukey Tests #

glht relies on Wald tests and often uses automated degrees of freedom corrections (like Satterthwaite) that don’t match the fixed error df from traditional ANOVA. The emmeans package is built for post-hoc comparisons in mixed models and lets you explicitly set parameters to match textbook results.

Try this workflow:

library(emmeans)
# Get estimated marginal means for each treatment
emm <- emmeans(model, ~ treatment)
# Run Tukey-adjusted comparisons using the textbook's error df (likely 27 here: (10-1)*(4-1))
tukey_comparisons <- pairs(emm, adjust = "tukey", df = 27)
print(tukey_comparisons)

Forcing the df to 27 ensures the test uses the same error degrees of freedom as Hays’ example, which is likely why your original p-values felt off.

3. Cross-Check the HSD Value Manually #

Since Hays reports HSD (honestly significant difference) instead of p-values, calculate it yourself to verify alignment:

# Extract MSE from your model (should match Hays' reported MSE)
mse <- anova(model)$"Mean Sq"[2] # Second row is the error term
# Calculate HSD using the textbook's parameters
hsd <- qtukey(0.95, nmeans = 4, df = 27) * sqrt(mse / 10) # n=10 subjects per treatment

Compare this calculated HSD to the textbook’s value. Any pairwise mean difference larger than this threshold should be significant—cross-reference this with the pairwise differences from emmeans to confirm results, even if p-values aren’t explicitly listed in the book.

Wrap-Up #

The core mismatch comes down to modern mixed-effects packages defaulting to flexible (but non-traditional) variance structures and df corrections, while Hays’ example uses strict classic ANOVA assumptions. By enforcing compound symmetry in lme and using emmeans with the correct fixed df, your post-hoc results should align perfectly with the textbook’s HSD values.

内容的提问来源于stack exchange，提问作者Sarah Stanley