在R中进行重复测量的ANOVA和TukeyHSD事后检验。

Question

我想对重复测量的ANOVA进行Tukey HSD事后检验。输入的公式"TukeyHSD"返回错误。我在论坛中找不到答案。我能寻求帮助吗？

"treat"是重复测量因子，"vo2"是因变量。

以下是产生此错误的脚本：

my_data <- data.frame(
  stringsAsFactors = FALSE,
  id = c(1L,2L,3L,4L, 5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L),
  treat = c("o","o","o","o","o","j","j","j","j","j","z","z","z","z","z","w","w","w","w","w"),
  vo2 = c("47.48","42.74","45.23","51.65","49.11","51.00","43.82","49.88","54.61","52.20","51.31",
          "47.56","50.69","54.88","55.01","51.89","46.10","50.98","53.62","52.77"))

summary(rm_result <- aov(vo2~factor(treat)+Error(factor(id)), data = my_data))
TukeyHSD(rm_result, "treat", ordered = TRUE)

- alunats

Answer 1

TukeyHSD() 无法处理重复测量的ANOVA结果中的 aovlist。作为替代方案，您可以使用相等混合效应模型进行拟合，例如lme4::lmer()，并使用 multcomp::glht() 进行事后检验。

my_data$vo2 <- as.numeric(my_data$vo2)
my_data$treat <- factor(my_data$treat)
m <- lme4::lmer(vo2 ~ treat + (1|id), data = my_data)
summary(multcomp::glht(m, linfct=mcp(treat="Tukey")))

# Simultaneous Tests for General Linear Hypotheses
# 
# Multiple Comparisons of Means: Tukey Contrasts
# 
# 
# Fit: lmer(formula = vo2 ~ treat + (1 | id), data = my_data)
# 
# Linear Hypotheses:
#            Estimate Std. Error z value Pr(>|z|)    
# o - j == 0   -3.060      0.583  -5.248   <0.001 ***
# w - j == 0    0.770      0.583   1.321   0.5497    
# z - j == 0    1.588      0.583   2.724   0.0327 *  
# w - o == 0    3.830      0.583   6.569   <0.001 ***
# z - o == 0    4.648      0.583   7.972   <0.001 ***
# z - w == 0    0.818      0.583   1.403   0.4974    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# (Adjusted p values reported -- single-step method)

将混合效应模型的ANOVA表与您的重复测量ANOVA结果进行比较，发现两种方法在处理“treat”变量时是等价的：

anova(m)
# Analysis of Variance Table
#       npar Sum Sq Mean Sq F value
# treat    3 61.775  20.592   24.23

summary(rm_result)
# Error: factor(id)
#           Df Sum Sq Mean Sq F value Pr(>F)
# Residuals  4  175.9   43.98               
# 
# Error: Within
#               Df Sum Sq Mean Sq F value   Pr(>F)    
# factor(treat)  3  61.78   20.59   24.23 2.22e-05 ***
# Residuals     12  10.20    0.85                     
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

- Robert Hacken

Answer 2

作为另一种选择，您也可以按照下面的示例进行操作。请注意，cld() 部分是可选的，仅通过“紧凑字母显示”（有关详细信息，请参见此处）尝试总结结果。

# data --------------------------------------------------------------------
my_data <- data.frame(
  stringsAsFactors = FALSE,
  id = c(1L,2L,3L,4L, 5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L,1L,2L,3L,4L,5L),
  treat = c("o","o","o","o","o","j","j","j","j","j","z","z","z","z","z","w","w","w","w","w"),
  vo2 = c("47.48","42.74","45.23","51.65","49.11","51.00","43.82","49.88","54.61","52.20","51.31",
          "47.56","50.69","54.88","55.01","51.89","46.10","50.98","53.62","52.77"))

my_data$vo2 <- as.numeric(my_data$vo2)
my_data$treat <- factor(my_data$treat)


# model -------------------------------------------------------------------
m <- lme4::lmer(vo2 ~ treat + (1|id), data = my_data)


# emmeans -----------------------------------------------------------------
library(emmeans)
emmeans <- emmeans(m, specs = "treat")
pairs(emmeans, adjust = "Tukey")
#>  contrast estimate    SE df t.ratio p.value
#>  j - o       3.060 0.583 12   5.248  0.0010
#>  j - w      -0.770 0.583 12  -1.321  0.5681
#>  j - z      -1.588 0.583 12  -2.724  0.0761
#>  o - w      -3.830 0.583 12  -6.569  0.0001
#>  o - z      -4.648 0.583 12  -7.972  <.0001
#>  w - z      -0.818 0.583 12  -1.403  0.5209
#> 
#> Degrees-of-freedom method: kenward-roger 
#> P value adjustment: tukey method for comparing a family of 4 estimates


# multcomp ----------------------------------------------------------------
library(multcomp)
library(multcompView)
cld(emmeans, Letters = letters)
#>  treat emmean   SE   df lower.CL upper.CL .group
#>  o       47.2 1.53 4.47     43.2     51.3  a    
#>  j       50.3 1.53 4.47     46.2     54.4   b   
#>  w       51.1 1.53 4.47     47.0     55.1   b   
#>  z       51.9 1.53 4.47     47.8     56.0   b   
#> 
#> Degrees-of-freedom method: kenward-roger 
#> Confidence level used: 0.95 
#> P value adjustment: tukey method for comparing a family of 4 estimates 
#> significance level used: alpha = 0.05 
#> NOTE: If two or more means share the same grouping symbol,
#>       then we cannot show them to be different.
#>       But we also did not show them to be the same.

^{创建于2022年12月20日，使用reprex v2.0.2}

- Paul Schmidt