Получение неправильных p-значений для теста Тьюки для одностороннего смешанного эффекта ANOVA
Я пытаюсь повторить учебный пример теста Тьюки для одностороннего смешанного эффекта ANOVA (из Статистики, Уильям Л. Хейс, стр. 581-583), но значения p, которые я получаю с помощью lme & glht, не имеют смысла
Исследование повторило измерения четырех процедур и 10 предметов
Данные
subject=c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6,
6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10)
treatment=c("a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3",
"a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2",
"a3", "a4", "a1", "a2", "a3", "a4", "a1", "a2", "a3", "a4", "a1",
"a2", "a3", "a4", "a1", "a2", "a3", "a4")
response=c(11, 9, 5, 17, 14, 12, 10, 18, 15, 7, 6, 21, 17, 10, 13, 22, 15, 7, 6,
15, 14, 8, 13, 22, 9, 6, 7, 15, 17, 11, 10, 19, 10, 13, 14, 23, 12,
8, 11, 20)
dataFrame=data.frame(subject, treatment, response)
Модель
library(nlme)
model = lme(response~ treatment,random=~1|subject,data=dataFrame)
anova(model)
numDF denDF F-value p-value
(Intercept) 1 27 375.9198 <.0001
treatment 3 27 43.4507 <.0001
Это F-значение достаточно близко к Хей (F = 43.41Я уверен, что с моей моделью все в порядке.
Тест Тьюки
library(multcomp)
glht.out =glht(model, mcp(treatment="Tukey"))
summary(glht.out)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = response ~ treatment, data = dataFrame, random = ~1 |
subject)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
a2 - a1 == 0 -4.300 1.006 -4.276 <0.001 ***
a3 - a1 == 0 -3.900 1.006 -3.879 <0.001 ***
a4 - a1 == 0 5.800 1.006 5.768 <0.001 ***
a3 - a2 == 0 0.400 1.006 0.398 0.979
a4 - a2 == 0 10.100 1.006 10.044 <0.001 ***
a4 - a3 == 0 9.700 1.006 9.647 <0.001 ***
Это не согласуется с книгой. Хей сообщил только о некоторых сравнениях и дал HSD и среднее значение, не p, но для обнаруженного контраста a3-a1 не было существенной разницы с HSD = 4.02 а также mean = 3.9 который по моим подсчетам имеет p.value = 1-ptukey(3.9*sqrt(10/8.27),4,9)=0.05719563,
Вывод R также не имеет смысла, потому что p-значение теста Тьюки, которое должно контролировать множественные сравнения, намного меньше, чем p-значение из парного t-теста (p=0,0142 с использованием t.test(c(11, 14, 15, 17, 15, 14, 9, 17, 10, 12),c(5, 10, 6, 13, 6, 13, 7, 10, 14, 11), paired=TRUE),
Любая идея, что я делаю неправильно и как я могу правильно выполнить тест Тьюки в R?
1 ответ
Я исследовал ваш вопрос, и это мои выводы. К сожалению, у меня нет книги Хэя, поэтому я просто сосредоточился на информации, которую вы разместили.
1) Относительно HSD = 4,02, сообщенного в книге Хэя. lme функция использует RELM метод подбора, который максимизирует ограниченную вероятность логарифмирования в то время как LM метод максимизирует логарифмическую вероятность, ?lme для получения дополнительной информации. Модель, соответствующая методу ML, показывает то же значение F.
model <- lme(response~treatment, random=~1|subject, data=dataFrame, method="ML")
anova.lme(model)
numDF denDF F-value p-value
(Intercept) 1 27 375.9227 <.0001
treatment 3 27 43.4506 <.0001
При выполнении нескольких сравнений с glht, z value ближе к той, о которой сообщается в книге 4.088, но контраст все же находит существенную разницу.
glht.out <- glht(model, mcp(treatment="Tukey"), alternative="two.sided")
summary(glht.out)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lme.formula(fixed = response ~ treatment, data = dataFrame, random = ~1 |
subject, method = "ML")
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
a2 - a1 == 0 -4.3000 0.9539 -4.508 < 1e-04 ***
a3 - a1 == 0 -3.9000 0.9539 -4.088 0.00027 ***
a4 - a1 == 0 5.8000 0.9539 6.080 < 1e-04 ***
a3 - a2 == 0 0.4000 0.9539 0.419 0.97520
a4 - a2 == 0 10.1000 0.9539 10.588 < 1e-04 ***
a4 - a3 == 0 9.7000 0.9539 10.168 < 1e-04 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
2) Относительно расчета р-значения. Предыдущий вывод показывает, что значения p корректируются с использованием single-step method, ?summary.glht описание объясняет, что любые методы, реализованные в adjusted а также p.adjust может быть использован. Вот пример того, как это изменить:
summary(glht.out, test = adjusted(type = "single-step"))
summary(glht.out, test = adjusted(type = "none"))
К сожалению, я не смог воспроизвести p-значения из вывода glht. Итак, я проверил TukeyHSD функция вместо Это требует aov модель выхода.
aovmodel <- aov(model)
summary(aovmodel)
Df Sum Sq Mean Sq F value Pr(>F)
treatment 3 659.0 219.67 26.95 2.57e-09 ***
Residuals 36 293.4 8.15
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
При запуске TukeyHSD мы получаем следующий вывод
TukeyHSD(aovmodel, "treatment")
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = model)
$treatment
diff lwr upr p adj
a2-a1 -4.3 -7.738482 -0.8615177 0.0093920
a3-a1 -3.9 -7.338482 -0.4615177 0.0210165
a4-a1 5.8 2.361518 9.2384823 0.0003378
a3-a2 0.4 -3.038482 3.8384823 0.9891641
a4-a2 10.1 6.661518 13.5384823 0.0000000
a4-a3 9.7 6.261518 13.1384823 0.0000000
Код легче понять, чем код от multcomp, по крайней мере, для меня. P-значения могут быть воспроизведены следующим образом:
MSE <- 8.15
center <- -3.9/sqrt(MSE/10) # -4.32002
ptukey(abs(center), 4, 36, lower.tail=FALSE)
[1] 0.02101645
Обратите внимание, что HSD теперь равен -4,32002, а степень свободы равна Nk = 40 - 4 = 36. Значение p равно 0,021, что немного выше, чем у парного t-критерия. Вот Gist для воспроизведения выходных данных.
Надеюсь это поможет