Tester l'hypothèse de normalité pour des mesures répétées anova? (en R)

11

En supposant donc qu'il est utile de tester l'hypothèse de normalité pour anova (voir 1 et 2 )

Comment peut-il être testé en R?

Je m'attendrais à faire quelque chose comme:

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

Ce qui ne fonctionne pas, car les "résidus" n'ont pas de méthode (ni de prédiction d'ailleurs) pour le cas des mesures répétées anova.

Alors, que faut-il faire dans ce cas?

Les résidus peuvent-ils simplement être extraits du même modèle d'ajustement sans le terme d'erreur? Je ne connais pas assez la littérature pour savoir si elle est valide ou non, merci d'avance pour toute suggestion.

Tal Galili
la source

Réponses:

5

Vous ne pouvez pas obtenir une réponse simple à residuals(npk.aovE)mais cela ne signifie pas qu'il n'y a pas de résidus dans cet objet. Faites stret voyez qu'à l'intérieur des niveaux, il y a encore des résidus. J'imagine que vous étiez le plus intéressé par le niveau "Within"

> 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 

Ma propre formation et pratique n'a pas consisté à utiliser des tests de normalité, mais à utiliser des tracés QQ et des tests parallèles avec des méthodes robustes.

DWin
la source
Merci Dwin. Je me demande lequel des différents résidus devrait être exploré (en plus de l'intérieur). Cheers, Tal
Tal Galili
npk.aovE est une liste de trois éléments. Les deux premiers sont des estimations de paramètres et la normalité n'est pas supposée pour eux, il ne semble donc pas approprié de tester quoi que ce soit sauf $ Within. names(npk.aovE)renvoie `[1]" (Intercept) "" block "" Within "`
DWin
@Dwin, pourriez-vous vérifier la dernière approche publiée par trev (dernière réponse)? Il utilise une sorte de projection pour calculer les résidus. C'est l'approche la plus simple pour moi afin de tracer un graphe avec homogénéité des variances entre les groupes.
toto_tico
@Dwin, qqplot semble également n'accepter que lm, et non aov.
toto_tico du
2

Une autre option serait d'utiliser la lmefonction du nlmepackage (puis de passer le modèle obtenu à anova). Vous pouvez utiliser residualssur sa sortie.

Nico
la source
1

Je pense que l'hypothèse de normalité peut être évaluée pour chacune des mesures répétées, avant d'effectuer l'analyse. Je remodèle le cadre de données afin que chaque colonne corresponde à une mesure répétée, puis effectue un test shapiro.test pour chacune de ces colonnes.

apply(cast(melt(npk,measure.vars="yield"), ...~N+P+K)[-c(1:2)],2,function(x) shapiro.test(x)$p.value)
George Dontas
la source
Merci gd047 - la question est de savoir ce que nous faisons quand nous avons un scénario plus complexe d'aov (rendement ~ N P K + erreur (bloc / (N + K)), npk) le test que vous proposez ferait-il le travail?
Tal Galili
Seriez-vous assez aimable pour expliquer la différence entre les scénarios? Malheureusement, je ne suis pas assez familier avec l'utilisation du terme d'erreur dans le modèle (au fait, pouvez-vous suggérer un bon livre à ce sujet?). Ce que je viens de proposer est la manière SPSS de vérifier l'hypothèse de normalité, telle que je l'ai apprise. Voir ceci comme un exemple goo.gl/p45Bx
George Dontas
Salut gd047. Merci pour le lien. Les choses que je sais sur ces modèles sont toutes liées à partir d'ici: r-statistics.com/2010/04/… Si vous apprenez à connaître d'autres ressources - j'adorerais les connaître. Cheers, Tal
Tal Galili
1

Venables et Ripley expliquent comment effectuer des diagnostics résiduels pour une conception à mesures répétées plus loin dans leur livre (p. 284), dans la section sur les effets aléatoires et mixtes.

La fonction des résidus (ou resid) est implémentée pour les résultats aov pour chaque strate:

à partir de leur exemple: oats.aov <- aov(Y ~ N + V + Error(B/V), data=oats, qr=T)

Pour obtenir les valeurs ajustées ou les résidus:

"Ainsi fitted(oats.aov[[4]])et resid(oats.aov[[4]])sont des vecteurs de longueur 54 représentant les valeurs ajustées et les résidus de la dernière strate."

Surtout, ils ajoutent:

"Il n'est pas possible de les associer uniquement aux tracés de l'expérience originale."

Pour les diagnostics, ils utilisent une projection:

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

Ils montrent également que le modèle peut être fait en utilisant lme, comme l'a signalé un autre utilisateur.

trev
la source
selon cela , ce devrait être [[3]] et non [[4]]. Je l'ai testé, et cela fonctionne juste pour [[3]]
toto_tico
1

Voici deux options, avec aov et avec lme (je pense que la 2e est préférée):

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)

L'exemple original est venu sans l'interaction ( Lme.mod <- lme(Y ~ N * V, random = ~1 | B/V, data = oats)) mais il semble fonctionner avec (et produire des résultats différents, donc il fait quelque chose).

Et c'est tout...

mais pour être complet:

1 - Les résumés du modèle

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

2 - Le test de Tukey avec des mesures répétées anova (3 heures à chercher !!).

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

3 - Les tracés de normalité et d'homoscédasticité

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)

entrez la description de l'image ici

toto_tico
la source