Тестирование предположения нормальности для повторных измерений anova? (в R)

Таким образом, предполагая, что есть смысл проверить допущение нормальности для ановы (см. 1 и 2 )

Как это можно проверить в R?

Я ожидал бы сделать что-то вроде:

## From Venables and Ripley (2002) p.165.
utils::data(npk, package="MASS")
npk.aovE <- aov(yield ~  N*P*K + Error(block), npk)
residuals(npk.aovE)
qqnorm(residuals(npk.aov))

Что не работает, так как «остатки» не имеют метода (и не предсказывают, в этом отношении) для случая повторных измерений anova.

Так что же делать в этом случае?

Могут ли остатки просто быть извлечены из одной и той же подходящей модели без условия ошибки? Я недостаточно знаком с литературой, чтобы знать, действительно ли это так или нет, заранее спасибо за любые предложения.

Вы можете не получить простой ответ, residuals(npk.aovE)но это не означает, что в этом объекте нет остатков. Сделайте strи убедитесь, что на уровнях все еще есть остатки. Я полагаю, что вас больше всего интересует уровень "Внутри"

> residuals(npk.aovE$Within)
          7           8           9          10          11          12 
 4.68058815  2.84725482  1.56432584 -5.46900749 -1.16900749 -3.90234083 
         13          14          15          16          17          18 
 5.08903669  1.28903669  0.35570336 -3.27762998 -4.19422371  1.80577629 
         19          20          21          22          23          24 
-3.12755705  0.03910962  2.60396981  1.13730314  2.77063648  4.63730314 

Мои собственные тренировки и практика не состояли в том, чтобы использовать тестирование нормальности, вместо этого использовать графики QQ и параллельное тестирование с надежными методами.

Другой вариант - использовать lmeфункцию nlmeпакета (а затем передать полученную модель anova). Вы можете использовать residualsна его выходе.

Я думаю, что допущение нормальности можно оценить для каждого из повторных измерений, прежде чем проводить анализ. Я бы изменил структуру данных таким образом, чтобы каждый столбец соответствовал повторяемой мере, а затем выполнил shapiro.test для каждого из этих столбцов.

apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)

Венейблс и Рипли объясняют, как проводить остаточную диагностику для схемы повторных измерений, позже в своей книге (стр. 284), в разделе о случайных и смешанных эффектах.

Функция остатков (или остаток) реализована для результатов aov для каждого слоя:

из их примера: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

Чтобы получить установленные значения или остатки:

«Таким образом, fitted(oats.aov[[4]])и resid(oats.aov[[4]])являются векторами длины 54, представляющими установленные значения и остатки от последней страты».

Важно отметить, что они добавляют:

«Невозможно связать их однозначно с сюжетами оригинального эксперимента».

Для диагностики они используют проекцию:

plot(fitted(oats.aov[[4]]), studres(oats.aov[[4]]))
abline(h=0, lty=2)
oats.pr <- proj(oats.aov)
qqnorm(oats.pr[[4]][, "Residuals"], ylab = "Stratum 4 residuals")
qqline(oats.pr[[4]][, "Residuals"])

Они также показывают, что модель может быть выполнена с использованием lme, как написал другой пользователь.

Вот два варианта, с aov и с lme (я думаю, что второй предпочтительнее):

require(MASS) ## for oats data set
require(nlme) ## for lme()
require(multcomp) ## for multiple comparison stuff

Aov.mod <- aov(Y ~ N * V + Error(B/V), data = oats)
the_residuals <- aov.out.pr[[3]][, "Residuals"]

Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)
the_residuals <- residuals(Lme.mod)

Исходный пример был получен без взаимодействия ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)), но, похоже, он работает с ним (и дает разные результаты, поэтому он что-то делает).

Вот и все...

но для полноты:

1 - Краткое изложение модели

summary(Aov.mod)
anova(Lme.mod)

2 - тест Тьюки с повторными измерениями anova (3 часа в поисках этого !!).

summary(Lme.mod)
summary(glht(Lme.mod, linfct=mcp(V="Tukey")))

3 - Графики нормальности и гомоскедастичности

par(mfrow=c(1,2)) #add room for the rotated labels
aov.out.pr <- proj(aov.mod)                                            
#oats$resi <- aov.out.pr[[3]][, "Residuals"]
oats$resi <- residuals(Lme.mod)
qqnorm(oats$resi, main="Normal Q-Q") # A quantile normal plot - good for checking normality
qqline(oats$resi)
boxplot(resi ~ interaction(N,V), main="Homoscedasticity", 
        xlab = "Code Categories", ylab = "Residuals", border = "white", 
        data=oats)
points(resi ~ interaction(N,V), pch = 1, 
       main="Homoscedasticity",  data=oats)