Testando a suposição de normalidade para medidas repetidas anova? (em R)

Portanto, supondo que haja um ponto em testar a suposição de normalidade para a anova (ver 1 e 2 )

Como pode ser testado em R?

Eu esperaria fazer algo como:

## 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))

O que não funciona, já que os "resíduos" não têm um método (nem prevê, nesse caso) para o caso de medidas repetidas anova.

Então, o que deve ser feito neste caso?

Os resíduos podem simplesmente ser extraídos do mesmo modelo de ajuste sem o termo Erro? Não estou familiarizado o suficiente com a literatura para saber se isso é válido ou não, desde já, obrigado por qualquer sugestão.

Você pode não receber uma resposta simples, residuals(npk.aovE)mas isso não significa que não haja resíduos nesse objeto. Faça stre veja que dentro dos níveis ainda existem resíduos. Eu imagino que você estivesse mais interessado no nível "Dentro"

> 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 

Meu próprio treinamento e prática não foi usar o teste de normalidade, em vez disso, usar gráficos de QQ e testes paralelos com métodos robustos.

Outra opção seria usar a lmefunção do nlmepacote (e depois passar o modelo obtido para anova). Você pode usar residualsem sua saída.

Penso que a suposição de normalidade pode ser avaliada para cada uma das medidas repetidas antes de realizar a análise. Remodelaria o quadro de dados para que cada coluna correspondesse a uma medida repetida e, em seguida, realizasse um shapiro.test para cada uma dessas colunas.

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

Venables e Ripley explicam como fazer diagnósticos residuais para um projeto de medidas repetidas posteriormente em seu livro (p. 284), na seção sobre efeitos aleatórios e mistos.

A função de resíduos (ou resid) é implementada para os resultados de Aov para cada estrato:

do exemplo deles: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

Para obter os valores ou resíduos ajustados:

"Deste modo fitted(oats.aov[[4]])e resid(oats.aov[[4]])são vectores de comprimento 54 que representam os valores ajustados e resíduos a partir do último estrato."

Importante, eles acrescentam:

"Não é possível associá-los exclusivamente às parcelas do experimento original".

Para diagnósticos, eles usam uma projeção:

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"])

Eles também mostram que o modelo pode ser feito usando o lme, como outro usuário postou.

Aqui estão duas opções, com aov e com lme (acho que a segunda é a preferida):

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)

O exemplo original veio sem a interação ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)), mas parece estar trabalhando com ela (e produzindo resultados diferentes, por isso está fazendo algo).

E é isso...

mas para ser completo:

1 - Os resumos do modelo

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

2 - O teste Tukey com medidas repetidas anova (3 horas procurando por isso !!).

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

3 - As parcelas de normalidade e homoscedasticidade

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)