Analyse de variance mixte: pas d'alternative à aov ?

[Stéphane Adamowicz]

Lorsqu'une analyse de variance comporte plusieurs termes d'erreurs, il est nécessaire de les modéliser et R semble proposer diverses solutions. Malheureusement, elles sont très pauvrement documentées pour le simple utilisateur que je suis. Si j'ai bien compris, la fonction aov() propose deux syntaxes pour les termes d'erreur :
1) Error(Facteur.aleatoire / Facteur.fixe) : typique du split-plot, cette syntaxe s'emploie lorsque le facteur aléatoire est au sommet de la hiérarchie. Elle semble aisée à transcrire dans les fonctions alternatives lme() ou lmer(), et on en trouve de nombreux exemples sur le net.
2) Error(Facteur.aleatoire %in% Facteur.fixe) : À l'inverse de la précédente, cette syntaxe exprime que le facteur aléatoire est niché dans le facteur fixe. Je n'ai pas compris comment cette disposition peut être spécifiée dans les fonctions lme(), lmer() ou ezANOVA(), ni trouvé d'exemples semblables traités avec elles dans la littérature ou sur le net.

Pourtant, le problème semble banal si l'on considère l'exemple suivant qui motive la présente requête : Nous avons cultivé des sujets (des plantes) sous divers régimes nutritionnels (Facteur fixe Nutrition) dans un dispositif aléatoire équilibré. À la récolte, chaque plante a été séparée en trois classes d'organes (Facteur fixe Organe: Feuilles, Tiges et Racines) qui ont été analysées pour leur composition minérale et organique. Nous attendons de l'analyse statistique de mesurer la sensibilité des organes à la nutrition, et si cette sensibilité est la même pour tous les organes (interaction Organe x Nutrition). En somme, les facteurs fixes "Nutrition" et "Organe" sont croisés, mais les Organes sont appariés sur les mêmes plantes, elle-mêmes nichées dans la nutrition. Voici, à titre d'exemple, l'analyse de variance de la teneur en saccharose :

Code : Tout sélectionner

> summary(aov(Saccharose ~ Nutrition*Organe + Error(Plant %in% Nutrition), data=Data))

Error: Plant:Nutrition
          Df Sum Sq Mean Sq F value   Pr(>F)   
Nutrition  3 0.7754 0.25847   43.25 2.73e-05 ***
Residuals  8 0.0478 0.00598                     
---
Error: Within
                 Df  Sum Sq  Mean Sq F value Pr(>F)
Organe            2 0.01450 0.007252   1.007  0.387
Nutrition:Organe  6 0.06621 0.011036   1.532  0.231
Residuals        16 0.11526 0.007204 

Cette procédure me satisferait pleinement si elle ne souffrait pas de deux limitations :
1) les fonctions residuals() et fitted() ne sont pas disponibles lorsqu'aov() comporte un terme Error. Avec quelques lignes de programmation, on arrive néanmoins à les extraire, mais ce n'est pas très générique ...
2) Suite à un accident, l'un des composés analysés comporte 2 données manquantes. Dans ce cas, en plus des problèmes posés par les calculs de SCE de type I, aov() "disjoncte" en testant le facteur Organe deux fois, une fois par terme d'erreur ...

Quelqu'un pourrait-il m'aider à trouver une alternative à aov() pour ce type de modèle statistique ?

[Stéphane Laurent]

Hello,
Avec lme4 ça devrait ressembler à ça :

Code : Tout sélectionner

lmer(Saccharose ~ Nutrition*Organe + (Plant | Nutrition), data=Data)

[Eric Wajnberg]

Si j'ai bien compris, la fonction aov() propose deux syntaxes pour les termes d'erreur :
1) Error(Facteur.aleatoire / Facteur.fixe) : typique du split-plot, cette syntaxe s'emploie lorsque le facteur aléatoire est au sommet de la hiérarchie. Elle semble aisée à transcrire dans les fonctions alternatives lme() ou lmer(), et on en trouve de nombreux exemples sur le net.
2) Error(Facteur.aleatoire %in% Facteur.fixe) : À l'inverse de la précédente, cette syntaxe exprime que le facteur aléatoire est niché dans le facteur fixe. Je n'ai pas compris comment cette disposition peut être spécifiée dans les fonctions lme(), lmer() ou ezANOVA(), ni trouvé d'exemples semblables traités avec elles dans la littérature ou sur le net.

Salut Stephane, j'espère que tu vas bien!

D'un point de vue syntaxique, le slash "/" ou le "%in%" sont strictement synonymes. Il est strictement équivalent d'écrire:

Code : Tout sélectionner

Error(Facteur.aleatoire / Facteur.fixe)

ou bien

Code : Tout sélectionner

Error(Facteur.fixe %in% Facteur.aleatoire)

C'est inhérent à la syntaxe de l'écriture d'un modèle sous R, et ne dépend pas de la fonction d'ajustement que tu utilises (aov, glm, lm, etc.).

Hope this helps!

Eric.

[Stéphane Adamowicz]

Désolé de n'être d'accord avec aucun de vous: aucune de ces écritures n'est équivalente à:

Code : Tout sélectionner

( summary( aov(Saccharose.~Nutrition*Organe + Error(Plant%in%Nutrition), data=Data)))
Error: Plant:Nutrition
          Df Sum Sq Mean Sq F value   Pr(>F)   
Nutrition  3 0.7754 0.25847   43.25 2.73e-05 ***
Residuals  8 0.0478 0.00598                     

Error: Within
                 Df  Sum Sq  Mean Sq F value Pr(>F)
Organe            2 0.01450 0.007252   1.007  0.387
Nutrition:Organe  6 0.06621 0.011036   1.532  0.231
Residuals        16 0.11526 0.007204

Cette écriture, qui me paraît correcte donne 2 termes d'erreurs. Le facteur nutrition y est testé par rapport à la variabilité entre plantes sachant la nutrition.

Si on utilise la syntaxe alternative avec barre de division, en mettant le facteur aléatoire en premier, ce qui doit toujours être le cas avec aov(), on obtient 3 termes d'erreur et le F du facteur Nutrition est très surestimé !

Code : Tout sélectionner

> (summary(aov(Saccharose.pc~Nutrition*Organe + Error(Plant/Nutrition), data=Data)))

Error: Plant
          Df  Sum Sq Mean Sq F value Pr(>F)
Residuals  2 0.03389 0.01695               

Error: Plant:Nutrition
          Df Sum Sq Mean Sq F value   Pr(>F)   
Nutrition  3 0.7754 0.25847   111.5 1.19e-05 ***
Residuals  6 0.0139 0.00232                     

Error: Within
                 Df  Sum Sq  Mean Sq F value Pr(>F)
Organe            2 0.01450 0.007252   1.007  0.387
Nutrition:Organe  6 0.06621 0.011036   1.532  0.231
Residuals        16 0.11526 0.007204

Enfin, si on inverse le facteur aléatoire et le facteur fixe dans le terme d'erreur, la Nutrition prend le statut de facteur aléatoire, ce qui conduit à l'absurdité suivante :

Code : Tout sélectionner

> (summary(aov(Saccharose.pc~Nutrition*Organe + Error(Nutrition/Plant), data=Data)))
Error: Nutrition
          Df Sum Sq Mean Sq
Nutrition  3 0.7754  0.2585

Error: Nutrition:Plant
          Df  Sum Sq  Mean Sq F value Pr(>F)
Residuals  8 0.04781 0.005976               

Error: Within
                 Df  Sum Sq  Mean Sq F value Pr(>F)
Organe            2 0.01450 0.007252   1.007  0.387
Nutrition:Organe  6 0.06621 0.011036   1.532  0.231
Residuals        16 0.11526 0.007204

Il me semble donc que les deux écritures ne sont absolument pas équivalentes.

Quant à la syntaxe proposée avec lmer, elle me paraît avoir des problèmes lorsqu'un facteur aléatoire est niché dans un facteur fixe:

Code : Tout sélectionner

> anova(lmer(Saccharose ~ Nutrition*Organe + (Plant | Nutrition), data=Data))
Analysis of Variance Table
                 Df  Sum Sq  Mean Sq F value
Nutrition         3 0.36314 0.121045 19.2443
Organe            2 0.01450 0.007252  1.1530
Nutrition:Organe  6 0.06621 0.011036  1.7545

En tout cas, merci d'avoir essayé de m'aider, et si d'autres idées vous viennent, je suis preneur ...

PS pour Eric : Oui merci, je vais bien. En tout cas mieux que mon code R !

[Eric Wajnberg]

Stephane,

Regarde bien en détail tes sorties pour:

Code : Tout sélectionner

( summary( aov(Saccharose.pc~Nutrition*Organe + Error(Plant%in%Nutrition), data=Data)))

et

Code : Tout sélectionner

(summary(aov(Saccharose.pc~Nutrition*Organe + Error(Nutrition/Plant), data=Data)))

Tu retrouves bien les mêmes résultats. Les deux syntaxes sont équivalentes.

Amitiés, Eric.

[Maxime Hervé]

Bonjour,
pour insister sur un point déjà évoqué, Plant%in%Nutrition équivaut exactement à Nutrition/Plant, et non pas à Plant/Nutrition.
Pour la syntaxe lmer(), je dirai plutôt que l'équivalent de

Code : Tout sélectionner

aov(Saccharose.pc~Nutrition*Organe + Error(Plant%in%Nutrition), data=Data)

est

Code : Tout sélectionner

lmer(Saccharose.pc~Nutrition*Organe + (1|Nutrition/Plant), data=Data)

Maxime

[Stéphane Adamowicz]

Certes, Éric,

mais tout dépend de ce qu'on entend par syntaxe équivalente. Il semble que tu considères qu'elles le sont parce qu'elles font la même décomposition des SCE. Mais on n'a pas besoin d'un terme Error() pour cela. Ensuite, on a toujours la ressource de terminer une analyse à la main ...
L'unique intérêt de ce terme Error() est de de faire les bons appariements de facteurs et de résidus. En ce sens les 2 syntaxes diffèrent réellement.

Mais, pour intéressante que soit notre discussion, ces considérations ne nous apportent pas d'alternative à aov() en cas de données manquantes ...

[Eric Wajnberg]

Certes, Éric,

mais tout dépend de ce qu'on entend par syntaxe équivalente. Il semble que tu considères qu'elles le sont parce qu'elles font la même décomposition des SCE. Mais on n'a pas besoin d'un terme Error() pour cela. Ensuite, on a toujours la ressource de terminer une analyse à la main ...
L'unique intérêt de ce terme Error() est de de faire les bons appariements de facteurs et de résidus. En ce sens les 2 syntaxes diffèrent réellement.

Mais, pour intéressante que soit notre discussion, ces considérations ne nous apportent pas d'alternative à aov() en cas de données manquantes ...

Les deux syntaxes sont équivalentes dans le sens où l'utilisateur demande avec ces deux syntaxes la même chose au logiciel (le calcul derrière est strictement le même).

Eric.

[Stéphane Adamowicz]

Merci pour cette suggestion dont voici le résultat:

Code : Tout sélectionner

> anova(lmer(Saccharose.pc~Nutrition*Organe + (1|Nutrition/Plant), data=Data))
Analysis of Variance Table
                 Df   Sum Sq  Mean Sq F value
Nutrition         3 0.211487 0.070496 10.3755
Organe            2 0.014504 0.007252  1.0673
Nutrition:Organe  6 0.066213 0.011036  1.6242

Comme vous voyez, c'est encore une analyse différente. Le test du facteur Nutrition est faux. Je n'en suis pas vraiment étonné parce que j'ai cru comprendre que le premier facteur intervenant après la barre verticale est interprété comme aléatoire.

Et c'est bien mon problème. Je ne suis pas parvenu à comprendre comment nicher un facteur aléatoire dans un facteur fixe avec lme() et lmer(). Mais d'un autre côté, il n'est écrit nul part que ce n'est pas possible ...

Ceci me paraît symptomatique de la faiblesse de la documentation en ligne et des ouvrages disponibles : j'ai déjà eu deux propositions de syntaxes différentes pour un même problème, a priori très simple puisqu'il ne s'agit que d'une analyse de variance à 2 facteurs ! Nous avons déjà commencé à faire de la combinatoire et à procéder par essais-erreurs ...

[Stéphane Laurent]

En fait je ne comprends pas la syntaxe aov() avec Error(). Pourquoi y'a 2 résiduelles dans ton exemple ?

[Stéphane Adamowicz]

La nécessité des deux termes d'erreur tient à ceci :
1) Les nutritions sont aléatoirement distribuées sur les plants. Le facteur Nutrition se teste naturellement par rapport à la variabilité entre plants sachant la nutrition ;
2) En revanche, l'interaction Organe x Nutrition n'est pas aléatoirement distribuée : les Organes étant appariés sur les mêmes plants, ils ne sont pas indépendants les uns des autres.

Cette disposition n'est pas si différente d'un split-plot, où la Nutrition tiendrait le rôle de Grand Bloc et le Plant, celui de sous-bloc. C'est sur la nature aléatoire ou fixe des éléments de la hiérarchie que se fait la différence avec un vrai split-plot.

[Nicolas Péru]

Je ne suis pas un spécialiste de la question mais je me suis demandé si les effets sont bien nichés au sens strict du terme à savoir que les catégories du facteur niché ne sont bien présentes qu'une seule fois dans celles du facteur au dessus hiérarchiquement.

La gestion des effets nichés est particulière dans lme4 ou encore nlme. Souvent, il n'y a pas besoin de préciser. Mais généralement, il est vrai que la question est traitée pour les effets aléatoires et pas pour le cas d'un effets aléatoire niché dans un effet fixe : ce n'est pas si courant que ça car il y a quand même pas mal de bouquin qui abordent ces questions et dans ceux que je possèdent, il n'y a rien là dessus (en tout cas en passant rapidement...)

Nicolas

[Stéphane Adamowicz]

Merci beaucoup de prendre le temps de réfléchir à ce type de problème.

N'étant pas non plus ce qu'on peut appeler un spécialiste, je me suis aussi posé ce genre de question, et j'ai aussi des doutes sur le vocabulaire qu'il convient d'utiliser. Je ne suis sûr que de deux choses:
1) du statut des facteurs (fixe ou aléatoire) dans cette étude, ce qui n'est pas toujours aussi tranché ;
2) de la nécessité de 2 termes d'erreur.

J'ai été surpris de ne pas trouver d'exemple comparable traité avec R. Pourtant, la littérature agronomique et celle de physiologie végétale regorgent d'exemples semblables. Pensez, par exemple, aux études en fonction de l'ambiance sous serres: effets de la brumisation, de l'enrichissement de l'air en CO2, de la température etc. On y compare des serres subissant chacune un certain traitement d'ambiance. Au sein de chaque serre, on compare, par exemple, plusieurs variétés végétales dans un dispositif en blocs complets équilibrés. Si on considère les serres en tant que facteur aléatoire, on a un split-plot typique, traitable avec aov et lmer qui donnent les mêmes analyses. En revanche, si on considère la serre en tant que facteur fixe, parce qu'on l'identifie à un traitement climatique expérimental, seul aov semble permettre de résoudre le problème. Je me dois d'ajouter que, dans les publications, les "matériels et méthodes" se contentent souvent de dire que les données ont été traitées par analyse de variance et de préciser le nom du logiciel utilisé, ce qui ne nous aide pas beaucoup ...

Stéphane

[Nicolas Péru]

Voilà ce que j'ai récupéré de mon côté quand je cherchais des définition aussi précises que possible du terme "niché".

Tant qu'à faire, on va citer un grand nom :)

We have seen an example of multiple random-effects expressions for the same grouping factor. It is more common to create random-effects expressions associated with different grouping factors. For example in a multicenter trial we may have observations over time on each of several patients at each of several institutions and it could be appropriate to use random effects for both patient and institution. If each patient is observed at only one institution
we say that the levels of patient are nested within the levels of institution; otherwise the factors are nonnested. Notice that it only takes one patient migrating between institutions to cause the grouping factors to lose the nesting property.

En comparant ce design et le tien, tu devrais pouvoir répondre à ta question de définition au moins du point de vue R.

Nicolas

[Stéphane Adamowicz]

Oui, d'accord ...

Enfin, je ne sais plus si nous sommes d'accord ou non !

Dans mon message initial, j'écrivais que les Plantes sont nichées dans les Nutritions, puisque chaque plante n'est affectée qu'à une seule Nutrition. Ce me semble correspondre à la définition du grand homme.

Mais il commence par se placer dans le cadre habituel (common) de relations entre facteurs aléatoires. Est-ce que cela exclut obligatoirement le cadre présent (qui est common dans mon domaine de recherches) ?

J'avais simplement utilisé le terme niché par opposition au terme croisé, ne connaissant pas de troisième mot:

Pour résumer:
1) Les facteurs Organe et Nutrition sont croisés, parce que chaque catégorie d'organes est exposée à chaque Nutrition ;
2) Les Organes sont croisés aux plantes, parce que chaque plante possède toutes les catégories d'organes ;
3) Les plantes ne sont pas croisées aux Nutritions, parce que chaque plante n'est exposée qu'à une seule nutrition.

Stéphane