Я хотел бы провести апостериорные тесты Tukey HSD для повторного измерения ANOVA. Введенная формула "TukeyHSD" возвращает мне ошибку. Не могу найти ответ на форуме. Могу ли я попросить о помощи?
«лечить» — фактор повторных измерений, «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)
TukeyHSD()
не может работать с aovlist
результатом повторных измерений ANOVA. В качестве альтернативы вы можете подобрать эквивалентную модель смешанных эффектов, например. 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
В качестве альтернативы вы также можете сделать это, как в репрексе ниже. Обратите внимание, что часть 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.
Created on 2022-12-20 with reprex v2.0.2