anova hiérarchique, avec facteur aléatoire et design déséquilibré

[mathilde dufay]

Bonjour à tous,

après un peu de recherches dans les docs et sur les forums, je ne parviens toujours pas à faire l'analyse suivante:
j'ai des données de nombre de grains de pollen par étamine, mesuré sur 98 plantes appartenant à 3 groupes (31 plantes G1, 30 plantes G2, 34 plantes M et 3 plantes non assignées), et je dispose de 3-4 réplicas par plante.
Je souhaite tester l'effet du groupe et de l'individu niché dans le groupe, l'effet individu étant aléatoire.

j'ai essayé plusieurs choses:

1. model1<-aov(Nb_pollen~Sexe+Error(Sexe/Indiv),data=poll_rep)

et j'obtiens le message suivant :

Warning message:
In aov(Nb_pollen ~ Sexe + Error(Sexe/Indiv), data = poll_rep) :
le modèle Error() est singulier

2. J'ai lu ensuite qqpart (sans trop comprendre pourquoi) qu'il fallait spécifier autrement le terme d'erreur:

model1<-aov(Nb_pollen~Sexe+Error(Indiv),data=poll_rep)

et j'obtiens le message suivant :

Error in qr.qty(qr.e, resp) :
'qr' and 'y' must have the same number of rows
>
est ce que c'est parce que mon design est déséquilibré ?

3. Je suis ensuite allée chercher du côté de la fonction lme qui apparemment autorise les designs déséquilibrés (??) mais je n'arrive pas à trouver de doc simple pour utiliser cette fonction donc j'ai probablement fait des erreurs:

model2<-lme(Nb_pollen~Indiv,random=~Indiv,na.action=na.omit)

mais j'obtiens:

Error in getGroups.data.frame(dataMix, groups) :
formule incorrecte pour les groupes
>

et je ne vois pas comment formuler autrement...

si qqn voit ce que je fais de mal, et surtout comment faire ce type d'analyse que je vais souvent etre amenée à réaliser régulièrement je suis preneuse de toute suggestion!
Merci d'avance,

Mathilde

[Florent Aubry]

Voilà la procédure que je conseille pour construire un modèle à effets mixtes avec lme ou lmer :
1) commencer par écrire la formule comme si tous les facteurs étaient fixes ;
2) à partir de cette formule construire la formule de l'effet aléatoire (random= pour lme, + (formule) pour lme4) de la façon suivante, en appelant 1 le facteur aléatoire :
- pour les effets principaux des facteurs aléatoires, on aura ~ 1 | A puisque l'effet aléatoire signifie qu'il existe une ligne de base (ordonnée à l'origine) aléatoire dépendant du niveau du facteur (à lire comme ordonnée aléatoire par niveau du facteur A) ;
- pour les interactions entre un prédicteur fixe X et le facteur aléatoire ; ~ X | A qui signifie qu'il y a un effet aléatoire dépendant des niveaux de A par paramètre estimé de X.

Exemple, on mesure dans le temps T la variation du variable Y et on a différents échantillons E qu'on suit dans le temps. Pour chaque échantillon e, on a y.e = a.e t + b.e. L'hypothèse est que y = a t + b. De ce fait, on aura :
- a.e = a + da.e ; da.e étant une pente propre à l'échantillon en excès ou défaut par rapport à la pente moyenne
- b.e = b + db.e, db.e étant l'excès ou défaut propre à l'échantillon par rapport à l'ordonnée à l'originee moyenne.
On a donc le modèle fixe de départ (cf. supra) :
Y ~ T * E ou Y ~ T + E + T:E ou Y ~ T + (1 + T):E
E étant un facteur aléatoire, on aura :
lme( Y ~ T, data=donnees, random=~ T | E)


Notes :
- Conformément à la syntaxe des formules sous R, ~ X | A est équivalent à (1 + X) | A donc, si on ne veut pas d'ordonnée à l'origine aléatoire, il faut écrire : (0 + X) | A.
- Si X est un facteur et que chaque niveau de A appartient à un et un seul niveau de X, alors, il n'y a pas d'interaction entre X et A ; par exemple le sexe d'un individu est unique, la formule Sexe * Individu n'a pas de sens et seule la formule Sexe + Individu a un sens. Si l'individu est un facteur aléatoire, alors la formule fixe est Y ~ Sexe et la formule aléatoire est ~ 1 | Individu.
- Un facteur ne peut pas être à la fois fixe et aléatoire, c'est-à-dire, dans les formules de R, se retrouver comme prédicteur de la formule fixe et à droite de | dans la formule aléatoire.

J'espère que cela répond à ta question pour lme que je conseille plutôt qu'aov car lme gère beaucoup mieux qu'aov les design déséquilibrés et, de plus, permet beaucoup plus facilement les tests post-hoc.

[Stéphane Adamowicz]

Bonjour Mathilde,

d'abord merci de nous proposer une ANOVA qui modélise le nombre de grains de pollen en fonction du sexe !

plus sérieusement,je vous confirme que la fonction aov utilisée avec un terme Error ne donne de résultats cohérents que pour les dispositifs équilibrés. De plus, votre modèle est mal écrit. Au lieu de

Code : Tout sélectionner

model1<-aov(Nb_pollen~Sexe+Error(Sexe/Indiv),data=poll_rep)

il aurait fallu écrire :

Code : Tout sélectionner

model1 <- aov(Nb_pollen ~ Sexe + Error(Indiv), data= poll_rep)

Mais mon petit doigt me dit qu'une étude statistique n'est peut être pas nécessaire dans ce cas ...

Les règles d'écriture du modèle sont les suivantes pour aov : y ~ F1*F2*F3*B1*B2*B3 + Error(sujet / (F1*F2*F3))
     Fi et Bi sont des facteurs fixes dont on veut connaître la signification ;
     sujet est un facteur aléatoire, typiquement un bloc ;
     Fi sont les facteurs "within" sujet, i.e. nichés dans sujet
     Bi sont les facteurs "between" sujet
Le symbole de division "/" dans Error sert à décrire la hiérarchie des facteurs fixes et aléatoires, et on peut en mettre plusieurs, par exemple :
     y ~ F1*F2*F3*B1*B2*B3 + Error(sujet / F1 / (F2+F3))
     Cette écriture est typique d'un split-plot où F1 a été distribué aléatoirement dans le sujet (un bloc, par ex.), et où F1 a été subdivisé pour y distribuer les diverses combinaisons de F2 & F3

La fonction lme de la bibliothèque nlme convient à un dispositif déséquilibré. Mais votre modèle n'est pas cohérent puisqu'il met Indiv à la fois comme facteur fixe et comme facteur aléatoire. Au lieu de :

Code : Tout sélectionner

model2<-lme(Nb_pollen~Indiv,random=~Indiv, na.action=na.omit)

il faudrait écrire, en s'inspirant du model1 (et en n'oubliant pas le champ data !)

Code : Tout sélectionner

model2 <- lme(fixed= Nb_pollen ~ Sexe, random= ~ 1|Indiv, , data= poll_rep, na.action=na.omit)

[mathilde dufay]

merci à tous les deux!!
je vais suivre vos conseils et je reviendrai si j'ai des pb supplémentaires (ou juste pour vous dire que ca a marché!)

Mathilde

ps: oui le pollen en fonction du sexe peut paraitre bizarre, mais ici pour la petite histoire les deux niveaux de facteurs sont "male" et "hermaphrodite"...
merci encore!