Discussion :

[Margotte02]

Bonjour,

j'aimerais faire une analyse split plot à 3 facteurs sur R.

Description de mon fichier: j'ai 3 variables mesurées:
- variable 1: la taille des racines (mm)
-variable 2: le poids des racines (g)
-variable 3: le poids des plantes (Kg)

et les variables qui peuvent influencer/expliquer ces mesures:
-variété (N=5) (= unité expérimentale)
-traitement (N=4) (3 traitements et un témoin non traité) (= grande parcelle)
-Field treatment (oui/non) (= petite parcelle)
-année de l'essai (N=3 pour variable 1 et 2 et N=2 pour variable 3) = répétitions / random factor

Mes 3 facteurs sont donc: variété, traitement et field-treatment et ma répétition est l'année.

J'aimerais faire un split pot pour chacune 3 variables mesurées, j'ai essayé pour la première (poids_racines) en faisant le code ci-dessous:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
library(lme4)
 
splitplot4.mod <- lmer(poids_racines ~ Traitement * field.treatment * Variete +(1|annee)+(Traitement|annee), data=fi)

et j'ai eu le message d'erreur suivant:

boundary (singular) fit: see ?isSingular ??

Pouvez-vous m'aider s'il vous plait?

J'ai également une autre question: comment utiliser les tests de Shapiro et Levene ici pour vérifier la normalité et l'égalité des variances?

Je mets mon fichier en PJ.

Merci par avance.
Salutations.

[faubry]

Tout simplement, il n'y a pas assez d'individus puisqu'on n'a qu'une occurrence du triplet (Traitement, field.treatment, Variete) par annee. Ce modèle avec l'année comme effet aléatoire ne peut donc pas être analysé. Si on se restreint à un modèle à effets fixes excluant l'année, il n'y a que 3 occurrences par triplet, ce qui en fait un modèle non robuste et sans intérêt. Il faut donc revoir totalement tes ambitions d'analyse. Le problème est que même dans ce cas, l'analyse est peu fiable puisqu'il est évident que le modèle à analyser ne repose pas sur un plan d'analyse défini a priori mais sur une analyse a posteriori.

[Margotte02]

Tout simplement, il n'y a pas assez d'individus puisqu'on n'a qu'une occurrence du triplet (Traitement, field.treatment, Variete) par annee. Ce modèle avec l'année comme effet aléatoire ne peut donc pas être analysé. Si on se restreint à un modèle à effets fixes excluant l'année, il n'y a que 3 occurrences par triplet, ce qui en fait un modèle non robuste et sans intérêt. Il faut donc revoir totalement tes ambitions d'analyse. Le problème est que même dans ce cas, l'analyse est peu fiable puisqu'il est évident que le modèle à analyser ne repose pas sur un plan d'analyse défini a priori mais sur une analyse a posteriori.

Merci pour le retour. Oui je suis bien embêtée par le fait que ces essais été déjà en place à mon arrivée et le plan d'analyse n'a pas été défini à priori.

Dès le départ l'année a été prévu pour servir de répétition "random factor", l'idée est de moyenner les 3 années par la suite dans mes graphiques si cela est possible mais je dois vérifier qu'il n'y a pas d’interaction entre mon facteur année et mes autres facteurs.

Je suis partie sur une analyse split plot car c'est ce qui semble se rapprocher le plus de mon modèle avec mes donnés. Ce que je n'ai pas précisé aussi c'est que:

j'ai aussi réalisé les mesures (destructives) à 3 dates d'observations différentes mais pour plus de simplicité je pensais faire une analyse séparée pour chacune des dates d'observation et aussi une analyse pour chacune des variables 1, 2 et 3. Donc en commençant par la variable 1 j'aimerais voir l'effet de la variété, du traitement; du field.treatment et aussi de leurs interactions éventuelles sur la variable "poids des racines" et vérifier l'impact de l'année.

Est-ce que vous pensez que c'est possible et fiable si je fais mes analyses "variable par variable" comme ceci (avec ce code, cela fonctionne):

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

splitplot4.mod <- lmer(poids_racines ~ Traitement * field.treatment * Variete +(1|annee)+(1|annee:Traitement), data=fi)

[faubry]

1) Utiliser l'année comme facteur aléatoire signifie de poser comme hypothèse qu'il n'y a aucune différence d'une année sur l'autre.
2) Traitement ne peut pas être à la fois utilisé comme facteur fixe (poids_racines ~ Traitement * field.treatment * Variete) et comme composante d'un facteur aléatoire (1|annee:Traitement). Il faut remplacer cette dernière formule par (Traitement|annee) comme tu l'avais fait.
3) Quel est le lien entre les années et les 3 mesures destructives. Si c'est une mesure par année, alors ton modèle peut se tenir à la condition d'utiliser toutes les mesures ensemble et non analyser mesure par mesure. Mais dans ce cas, est-ce l'année qui est le facteur aléatoire. Dans ton fichier, cette notion de duplication de mesure n'apparait pas vraiment.

Pour l'instant, d'après ton fichier, mon conseil est i) d'exclure l'année de l'analyse et ii) de bien tester la robustesse de ton modèle fixe, 3 valeurs par valeur du triplet (Traitement, field.treatment, Variete) fait largement douter de la robustesse de l'analyse. Un minimum raisonnable serait autour d'une dizaine.

[Margotte02]

1) Utiliser l'année comme facteur aléatoire signifie de poser comme hypothèse qu'il n'y a aucune différence d'une année sur l'autre.
2) Traitement ne peut pas être à la fois utilisé comme facteur fixe (poids_racines ~ Traitement * field.treatment * Variete) et comme composante d'un facteur aléatoire (1|annee:Traitement). Il faut remplacer cette dernière formule par (Traitement|annee) comme tu l'avais fait.
3) Quel est le lien entre les années et les 3 mesures destructives. Si c'est une mesure par année, alors ton modèle peut se tenir à la condition d'utiliser toutes les mesures ensemble et non analyser mesure par mesure. Mais dans ce cas, est-ce l'année qui est le facteur aléatoire. Dans ton fichier, cette notion de duplication de mesure n'apparait pas vraiment.

Pour l'instant, d'après ton fichier, mon conseil est i) d'exclure l'année de l'analyse et ii) de bien tester la robustesse de ton modèle fixe, 3 valeurs par valeur du triplet (Traitement, field.treatment, Variete) fait largement douter de la robustesse de l'analyse. Un minimum raisonnable serait autour d'une dizaine.

Ok je comprends un peu mieux; alors si je pose comme hypothèse que l'année ne va pas influencer l'efficacité des traitement et du field.treatment et que je la prends comme facteur aléatoire je dois écrire le code comme ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

splitplot4.mod <- lmer(poids_racines ~ Traitement * field.treatment * Variete +(1|annee)+(Traitement|annee), data=fi)

par contre j'ai le message d'erreur

boundary (singular) fit: see ?isSingular

En effet je n'ai pas expliqué le lien entre les années et les 3 mesures destructives car j'ai fait 4 fichiers différents par date d'observation en pensant réaliser les analyses pour chacune des dates.
Pour chacune des 3 années j'ai réalisé les mesures à 0 mois (sur un control avant de faire des "Traitements", un control a été traité au champs et l'autre non et ce pour les 5 variétés), ensuite j'ai effectué des autres mesures à 3, 5 et 7 mois sur les 5 variétés traitées ou non au champs (petites parcelles) et aussi traités avec 4 "Traitements différents" dans les "grandes parcelles".

Les plantes sont jetées après les mesures.
Voici le fichier complet en PJ. Est-ce que si je fais un schéma de mon dispositif expérimental cela serait plus clair? j'ai du mal à m'expliquer.
En fait je n'ai que 3 valeurs par valeur du triplet car j'ai moyenné les mesures dans mon fichier pour certaines variables comme ceci:
- variable 1: la taille des racines (mm) --> la mesure est réalisé sur 25 plantes donc j'ai 25 mesures de racines que j'ai moyenné dans mon fichier pour n'avoir qu'une mesure
-variable 2: le poids des racines (g) --> ici je n'ai qu'une mesure car le poids des racines des 25 plantes à été mesuré en une fois pour l'année 1 et 2 et pour l'année 3 cela a été mesuré par lot de 5 et aussi par 25 (je n'ai donc gardé que la mesure par 25 pour avoir comme les 2 autres années)
-variable 3: le poids des plantes (Kg) --> ici pareil, je n'ai qu'une mesure car le poids des 25 plantes a été mesuré en une fois pour l'année 1 et 2 et pour l'année 3 cela a été mesuré par lot de 5 + 25

Je comprends bien que ce n'est pas le schéma idéal du tout. Ce qui serait envisageable c'est de faire un test lmer uniquement pour la taille des racines pour laquelle je peux avoir 25 mesures par variété, traitement, field.treatment en gardant la formule de code ci-dessus (mes années resteront des réplicats c'est juste?) et de faire un test plus simple pour les 2 autres paramètres?

Je peux aussi envisager de "laisser tomber" le paramètre "poids des plantes" qui contient trop de NA et la date d'observation "o mois" qui en contient beaucoup également afin d'avoir un fichier plus "propre" et de pouvoir faire une analyse correcte.

Merci pour votre temps.

[Margotte02]

je complète mon message précédent, je viens d'essayer autre chose:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
 
#fichier entier:
fi<-read.csv2("complete_file.csv", na.strings = "NA",header = TRUE)
#fichier uniquement pour l'observation à  mois
test<-subset(fi, periode=="5months", select= c (Variete, annee,Traitement, field.treatment, poids_racines, taille_racines, poids_plantes))
head(test)
fi<-test
library(nlme)
splitplot4.mod2 <- lme(poids_racines ~ Traitement* field.treatment * Variete, random=~1|annee, data=fi)
anova(splitplot4.mod2 )

en faisant comme ceci j'obtiens une pvalue pour chacun des paramètre et pour les interactions, je peux faire ça pour chacune des dates d'observation et pour chacune des variables mesurées afin de simplifier mon modèle?

Par contre comment puis-je avoir le détails de la significativité, c à d que je vois que le traitement à une influence significative sur le poids des racines mais je ne vois pas l'influence des traitements un par un ou peut-être que je pourrais aussi voir l'influence des traitements par rapport au control?

[faubry]

Pour moi l'analyse n'a de sens que sur le fichier complet. En effet, quand il y a des mesures répétées (c'est le cas des mesures année par année puis dans chaque année, période par période), faire des Anova sur chacune des mesures ou faire une Anova sur l'ensemble des mesures par une Anova pour mesures répétées ou en utilisant une approche par modèle mixte ne donne pas la même chose. La raison est le phénomène dit de comparaisons multiples : les différentes répétition peuvent être corrélée ce que prend en compte (au moins en partie) une approche mesures répétées alors que les Anova individuelles font l’hypothèse de mesures non corrélée.

Si je comprends ton plan à partir du fichier complet, chaque variétés est vu 3 ans de suite à quatre moments différents (facteur période) et subissent 3 traitements + un controle au temps 0. A chaque fois la taille et le poids de la racine sont mesurées. Comment intervient le facteur fiels.treatment est moins clair pour mois. Si je comprends bien ton fichier, à chaque moment, il y a une mesure du poids avec une option filed.treatment vrai ou fausse. Est-ce bien cela ? Donc ton plan sous forme de formule R, indépendamment que ce soit des facteurs fixes ou aléatoires me semble être (j'expliquerais ultérieurement, si nécessaire, comment passer de cette formule d'un modèle à effets fixes [+ mesures répétées et/ou analyse conjointe de variables] à un modèle à effets mixtes) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

cbind( poids_racines, tailles_racine) ~ Variété * annee * periode * Traitement * filed.treatment

Parmi toutes les interactions, les quelles sont réellement pertinentes en sachant que :
1) des interactions d'ordre supérieurs à 2 sont généralement très difficiles à interpréter ; dans certains cas, s'il y a assez de données ou des hypothèses suffisamment fortes sur les données, on va jusqu'à l'ordre 3 mais très rarement au-delà ;
2) plus il y a de facteurs ou d'interactions à analyser, plus le nombre de degrés de liberté du modèle est important et donc moins il est robuste.

Tu ne considères que comme variable mesurée la variable poids_racine. Fais-tu une analyse séparée pour taille_racine ? Il serait intéressant d'estimer leur corrélation et alors de les analyses conjointement.

Il me semble que le seul facteur éventuellement aléatoire est l'année. De plus, considères-tu qu'il y a des interaction entre les facteurs que tu considères fixes et les facteurs aléatoires. Comment les justifies-tu ta réponse à cette question ?

Mon conseil est donc que tu reprécises tes hypothèses sur les données et l'hypothèse que tu veux tester sinon je ne pense pas que tu puisses réellement t'en sortir.

[Margotte02]

Merci beaucoup pour votre temps et toutes ces explications, cela m'est vraiment utile.

Pour moi l'analyse n'a de sens que sur le fichier complet. En effet, quand il y a des mesures répétées (c'est le cas des mesures année par année puis dans chaque année, période par période), faire des Anova sur chacune des mesures ou faire une Anova sur l'ensemble des mesures par une Anova pour mesures répétées ou en utilisant une approche par modèle mixte ne donne pas la même chose. La raison est le phénomène dit de comparaisons multiples : les différentes répétition peuvent être corrélée ce que prend en compte (au moins en partie) une approche mesures répétées alors que les Anova individuelles font l’hypothèse de mesures non corrélée.

--> vous avez raison, je commence à y voir plus clair. Je vais donc "sacrifier" la période "0month" car elle contient beaucoup trop de NA et n'est pas indispensable, ensuite je vais travailler sur 2 fichiers séparés dans le but d'avoir des fichiers sans données manquantes comme ceci dans le nouveau fichier en pièce jointe:

• feuillet "file_2var_3an" : dans lequel j'ai mes 2 variables taille et poids des racines pour les 5 variétés traités ou non au champs (field.treatment) puis traités avec 4 produits (P1,P2,P3 et un control non traité) et les mesures taille et poids des racines ont été effectuées à 3 périodes (3,5 et 7 mois) par mesure destructive et ceci 3 années de suite
  
• puis le feuillet "file_3var_2an": dans lequel j'ai mes 3 variables taille et poids des racines + le poids des plantes pour les 5 variétés traités ou non au champs (field.treatment) puis traités avec 4 produits (P1,P2,P3 et un control non traité) et les mesures ont été effectuées à 3 périodes (3,5 et 7 mois) par mesure destructive et ceci uniquement 2 années de suite

Si je comprends ton plan à partir du fichier complet, chaque variétés est vu 3 ans de suite à quatre moments différents (facteur période) et subissent 3 traitements + un controle au temps 0. A chaque fois la taille et le poids de la racine sont mesurées. Comment intervient le facteur fiels.treatment est moins clair pour mois. Si je comprends bien ton fichier, à chaque moment, il y a une mesure du poids avec une option filed.treatment vrai ou fausse. Est-ce bien cela ?

Désolée mes explications n'étaient pas très claires. J'ai ajouté un feuillet "design" dans le fichier excel pour illustrer. Mes plantes de différentes variétés sont traitées ou pas au champs (field-treatment N=2, oui ou non) (ce que je considère comme petite parcelle ou sub-unit), ensuite les plantes des 5 variétés sont mises dans des enceintes et sont traités avec 4 modalités (produits p1, p2, p3 et un control non traité) (ce que je considère comme les grandes parcelles ou whole unit) et tout cela est répété 3 années. Au sein d'une année les mesures (taille des racines et poids des racines et poids des plantes) sont réalisées à 3 périodes d'observations: 3, 5 et 7 mois sur chaque variété traitée ou non au champs puis traitées avec P1 ou P2 ou P3 ou non traitée (sur 2 années pour la taille et le poids des racines et sur 1 année pour le poids des plantes).

Parmi toutes les interactions, les quelles sont réellement pertinentes en sachant que :
1) des interactions d'ordre supérieurs à 2 sont généralement très difficiles à interpréter ; dans certains cas, s'il y a assez de données ou des hypothèses suffisamment fortes sur les données, on va jusqu'à l'ordre 3 mais très rarement au-delà ;
2) plus il y a de facteurs ou d'interactions à analyser, plus le nombre de degrés de liberté du modèle est important et donc moins il est robuste.

Si je formule des hypothèses fixes ici:
hypothèse fixe 1: je vais considérer que l'année n'a pas d'effet sur les facteurs (field.treatment, traitement et variété) mais l'année peut influence les variables mesurées (poids,taille racines et poids plantes)
hypothèse fixe 2: je vais considérer que l'ordre des variétés dans les parcelles n'a pas d'influence


Hypothèses à tester (si on prend par exemple le fichier avec uniquement les 2 variables poids et taille des racines):

• Effet du field.treatment sur le poids et la taille des racines
• Effet du traitement sur le poids et la taille des racines
• Effet de la variété sur le poids et la taille des racines
• Effet de l'interaction field.treatment *traitement sur le poids et la taille des racines
• Effet de l'interaction field.treatment *variété sur le poids et la taille des racines
• Effet de l'interaction traitement *variété sur le poids et la taille des racines
• Effet de l'interaction field.treatment *traitement *variété sur le poids et la taille des racines

Par contre je comptais tester ces hypothèses à 3,5 et 7 mois séparément car cela fait beaucoup de facteurs et tout simplement parce que je ne savais pas comment faire pour intégrer le facteur "période" mais si cela est possible et que vous pensez que c'est mieux je suis à l'écoute

Il me semble que le seul facteur éventuellement aléatoire est l'année. De plus, considères-tu qu'il y a des interaction entre les facteurs que tu considères fixes et les facteurs aléatoires. Comment les justifies-tu ta réponse à cette question ?

Donc ici comme j'ai fixé que l'année n'avait pas d'influence c'est mon facteur alétoire et ma répétition et je ne dois pas tester les interactions avec l'année c'est juste?

Tu ne considères que comme variable mesurée la variable poids_racine. Fais-tu une analyse séparée pour taille_racine ? Il serait intéressant d'estimer leur corrélation et alors de les analyses conjointement.

Oui en effet ça serait intéressant mais comment faire cela? je peux analyser les 2 sur mon premier feuillet puis même les 3 (avec le poids des plantes) sur mon deuxième feuillet?

Donc ton plan sous forme de formule R, indépendamment que ce soit des facteurs fixes ou aléatoires me semble être (j'expliquerais ultérieurement, si nécessaire, comment passer de cette formule d'un modèle à effets fixes [+ mesures répétées et/ou analyse conjointe de variables] à un modèle à effets mixtes) :

ici je ne suis pas sûr de comprendre.

Encore merci pour votre temps.

[faubry]

hypothèse fixe 1: je vais considérer que l'année n'a pas d'effet sur les facteurs (field.treatment, traitement et variété)

L'année étant une variable indépendante au même titre que field.treatment, traitement et variété, cette variable ne peut pas avoir d'effet sur les autres facteurs. La question est donc de savoir s'il y a indépendance entre l'année et les autres facteurs et donc pas d'interaction, c'est-à-dire que l'effet du traitement, du field.traetmentet de la variété est indépendant de l'année ou s'il y a une interaction, c'est-à-dire que ces effets dépendent de l'année.
La première question est de savoir si les variations en fonction de l'année relève de fluctuations statistiques de mesure (effet aléatoire) ou d'une variation en fonction de l'année, ici de type série temporelle. Dans le premier cas, je me poserais la question de la pertinence de conserver ce facteur.
La seconde concerne le facteur période. Faire des analyses séparées revient à considérer que les résultats entre les périodes sont indépendants et ne permet pas d'analyser les variations d'une période sur l'autre. SI on considère que chaque échantillon est référencé par un identificateur ID, c'est-à-dire qu'un échantillon à le même identificateur durant toutes les périodes, on peut alors s'en servir comme facteur de regroupement, ce qui pourrait donner un modèle du type :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

~ (field.treatment + traitement + variété + période)^2 + (1 |ID)

Cependant, vu le nombre de ddl, il faut se poser la question de la robustesse de l'analyse et de la simplification du modèle.

Remarque annexe : un plan d'expérience est construit pour répondre à une et une seule question. Ensuite, il peut y avoir des sous questions mais ce sont des analyses post-hoc. En fait ta liste d'hypothèses à tester ne sont que des tests post-hoc de l'analyse plus générale. Par exemple, dans un modèle à effets fixe Y ~ A * B, si au sens de la somme des carrés de type II, au moins l'un des facteurs n'est pas significatif, alors l'interaction A:B ne peut pas l'être.

[Margotte02]

Donc si je considère que l'effet du traitement, du field.traetmentet de la variété est indépendant de l'année (c'est le cas ici), je peux mettre un ID identique pour une variété donnée * un traitement * un field treatment et j'aurais donc 3 fois l'ID dans mon fichier c'est juste?

Si il y a des différences de mesures entre années je considère que c'est dû à l’échantillonnage et au facteur humain qui mesure uniquement dans mon cas et donc dû aux fluctuations statistiques de mesure (effet aléatoire).

Concernant le facteur période je ne pense pas que je peux m'en servir comme facteur de regroupement car par exemple le poids de mes racines évolue dans le temps, je vais avoir 0g de racine à 3 mois pour certaines variétés et 10 g à 5 mois. Aussi certains traitements peuvent avoir un effet à 3 et 5 mois puis ne plus avoir d'effet à 7 mois donc peut-être que c'est mieux de séparer les périodes?

[faubry]

Donc si je considère que l'effet du traitement, du field.traetmentet de la variété est indépendant de l'année (c'est le cas ici), je peux mettre un ID identique pour une variété donnée * un traitement * un field treatment et j'aurais donc 3 fois l'ID dans mon fichier c'est juste?

Oui

Si il y a des différences de mesures entre années je considère que c'est dû à l’échantillonnage et au facteur humain qui mesure uniquement dans mon cas et donc dû aux fluctuations statistiques de mesure (effet aléatoire).

Mais cet effet, est-il suffisamment important pour être corrigé ou peut-on considérer qu'il biaise si peu les mesures qu'on peut le négliger ? Dans le premier cas, il faut donc ajouter à la formule cet effet aléatoire sous la forme + (1|annee). Cependant, vu déjà la complexité du modèle et son nombre de ddl, je suis sceptique sur son utilité et pense que d'introduire cet effet aléatoire met en péril la robustesse de l'analyse. Il est toujours préférable d'utiliser un modèle légèrement biaisé plutôt qu'un modèle peu robuste, c'est le principe de la régularisation et des méthodes dites robustes.

Concernant le facteur période je ne pense pas que je peux m'en servir comme facteur de regroupement car par exemple le poids de mes racines évolue dans le temps, je vais avoir 0g de racine à 3 mois pour certaines variétés et 10 g à 5 mois.

ARgument justifiant de l'utiliser comme effet fixe.

Aussi certains traitements peuvent avoir un effet à 3 et 5 mois puis ne plus avoir d'effet à 7 mois

Ce qui justifie pleinement d'analyser toutes les périodes ensemble pour pouvoir comparer leurs effets.

[Margotte02]

Mais cet effet, est-il suffisamment important pour être corrigé ou peut-on considérer qu'il biaise si peu les mesures qu'on peut le négliger ? Dans le premier cas, il faut donc ajouter à la formule cet effet aléatoire sous la forme + (1|annee). Cependant, vu déjà la complexité du modèle et son nombre de ddl, je suis sceptique sur son utilité et pense que d'introduire cet effet aléatoire met en péril la robustesse de l'analyse. Il est toujours préférable d'utiliser un modèle légèrement biaisé plutôt qu'un modèle peu robuste, c'est le principe de la régularisation et des méthodes dites robustes.

Non je pense que l'on peut considérer qu'il biaise si peu les mesures qu'on peut le négliger, je suis d'accord avec vous sur le fait qu'il faut simplifier.
L'année peut en revanche influence la taille ou poids des racines (par ex. l'an 1 on peut avoir des racines beaucoup plus importantes sur toutes les plantes mais je considère que ça sera le même effet sur toutes les variétés, field.treatment et traitements.

Peut-être que je peux tester les deux, c à d tester l'année comme effet aléatoire sous la forme + (1|annee) et tester sans l'année en utilisant les ID comme vous le proposez.

Je vais ajouter les ID sur mon fichier et je reviens vers vous.

[faubry]

L'année peut en revanche influence la taille ou poids des racines (par ex. l'an 1 on peut avoir des racines beaucoup plus importantes sur toutes les plantes

Dans ce cas, l'année ne peut pas être un facteur aléatoire car qui dit facteur aléatoire dit que l'effet de ce facteur est constant quelque soit l'année et que les variations constatées ne sont dues qu'à des fluctuations aléatoires. Si son effet est différent d'une année sur l'autre, c'est obligatoirement un facteur fixe.

mais je considère que ça sera le même effet sur toutes les variétés, field.treatment et traitements.

Cela signifie simplement qu'il n'y a pas d'interaction avec ces autres facteurs donc la formule serait du type (en ne considérant que des interactions d'ordre 2) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

~ (variétés + field.treatment + traitements)^2 + annee

Si tu testes toutes les périodes, la présence de ID est obligatoire puisque cette technique revient à remplacer une Anova à mesures répétées par un modèle mixte. La littérature montre que cette approche est pertinente, surtout en absence de sphéricité des variables dépendantes.

[Margotte02]

OK merci, juste une petite question intermédiaire, comment faire pour attribuer un ID identique aux répétitions des 3 années (à la main c'est un peu long ), faut il faire une boucle du style:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
 
 
if fichier$Variete == fichier$field.treatment == fichier$Traitement == fichier== periode,
 
alors --> attribute ID ......

[faubry]

Le plus simple à condition que les individus soient dans le même ordre quelque soit la période.
Un fichier par période :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
don.1 <- lecture de la première période
n.ind <- nrow( don.1)
id <- paste0( "ID", 1:n.ind)
don.1$ID <- id
don.1$periode <- identificateur période
#
don.2 <- lecture de la deuxième période
don.2$ID <- id
don.2$periode <- identificateur période
#...
 
# Concaténation des fichiers
don <- rbind( don.1, don.2, ...)

Il faut modifier le code en fonction de la logique qui a servi à construire les fichiers de données.

[Margotte02]

Je me perds un peu; pour être sûre de bien comprendre:

1) si on prend un fichier sans ID et que l'on considère l'année comme facteur fixe:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
Cela signifie simplement qu'il n'y a pas d'interaction avec ces autres facteurs donc la formule serait du type (en ne considérant que des interactions d'ordre 2) :
 ~ (variétés + field.treatment + traitements)^2 + annee

quand vous parlez d'ordre 2 ça veut dire que avec cette formule nous testons l'effet de la variété puis du field.treatment puis des traitemetns sur la taille des racines et aussi leurs interractions? Puis l'année est aussi testé pour voir si il y a un effet uniquement sur la taille des racines?
comme ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
test_poids<-poids_racines ~ (variétés + field.treatment + traitements)^2 + annee
test_taille<-poids_taille ~ (variétés + field.treatment + traitements)^2 + annee

Dans cette formule on ne considère pas la période? je dois faire cette formule pour des sous-fichiers par période c'est juste?

2)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
#il y a indépendance entre l'année et les autres facteurs et donc pas d'interaction, 
#c'est-à-dire que l'effet du traitement, du field.traetmentet de la variété est indépendant de l'année:
 
test_poids<-poids_racines~ (field.treatment + traitement + variété + période)^2 + (1 |ID)
test_taille<-poids_taille~ (field.treatment + traitement + variété + période)^2 + (1 |ID)

Ici je dois utiliser un ficher avec un ID pour 1 variété * 1 traitement * 1 field.treatment * 1 période et donc j'aurais 3 fois le même ID puisque j'ai 3 années c'est bien ça? On utilise l'ID pour avoir 3 réplicats par échantillon et on ne s'occupe plus de l'année?
Et avec cette formule on test l'effet du field.treatment, du traitement, de la variété et de la période et leurs interraction sur le poids des racines puis sur la taille.

[faubry]

Mes différentes formules étaient des illustrations. La bonne approche c'est le 2). Pour tester l'influence de l'année, il suffit de comparer les résultats des deux formules :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
~ (field.treatment + traitement + variété + période)^2 +annee +  (1 |ID)
~ (field.treatment + traitement + variété + période)^2 + (1 |ID)

Pour effectuer cette comparaison, il faut utiliser la fonction anova. Si la différence n'est pas significative, alors l'année peut être supprimé. Tu peux continuer à chercher des simplifications en utilisant la fonction stepAIC.

[Margotte02]

Cela fonctionne, merci beaucoup!!

j'ai mis les ID dans mon fichier que je viens de mettre en pièce-jointe, juste pour être sûre que c'est juste j'ai bien placé les ID? (exemple: pour (var 1, traitement P1, field.treatment=No et période = 3 mois) je lui met l'ID = 1 pour l'année 1 mais aussi pour l'année 2 et 3 puis pour (var 1, traitement P1, field.treatment=No et période = 5 mois) je met un ID 2, etc..

Ensuite je "nettoie mon fichier" pour ne pas avoir de NA donc je travail sur un subset pour les différentes variables mesurées:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
fi<-read.csv2("file.csv.csv", na.strings = "NA",header = TRUE)
dim(fi)
head(fi)
fi_2<-subset(fi, periode != "0month", select= c (ID, Variete, annee, periode , Traitement, field.treatment, poids_racines, taille_racines))
dim(fi_2)
head(fi_2)
str(fi_2)
fi<-fi_2
fi$poids_racines<-as.numeric(fi$poids_racines)
fi$taille_racines<-as.numeric(fi$taille_racines)

et je test les 2 formules ANOVA

ici je considère l'année uniquement comme une répétition sans la tester, formule1:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
aov_id_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^2 + (1 |fi$ID))
summary(aov_id_poids)

Ici l'année est mon facteur fixe, formule 2:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
aov_id_an_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^2 + fi$annee + (1 |fi$ID))
summary(aov_id_an_poids)

On obtient les résultats ci-dessous qui sont légèrement différents donc est-ce que je peux quand même ne pas tenir compte de l'année et l'utiliser comme répétition uniquement et je peux alors appliquer la formule 1 pour mon fichier sur 3 années sur la taille des racines puis sur le poids des racines et avec mon fichier avec 2 années sur le poids des plantes?
j'ai mis en rouge les différences:

> summary(aov_id_poids)
Df Sum Sq Mean Sq F value Pr(>F)
fi$field.treatment 1 133480 133480 67.246 6.20e-15 ***
fi$Traitement 3 472538 157513 79.354 < 2e-16 ***
fi$Variete 4 22593 5648 2.846 0.02424 *
fi$periode 2 79375 39688 19.994 6.71e-09 ***
fi$field.treatment:fi$Traitement 3 54575 18192 9.165 7.87e-06 ***
fi$field.treatment:fi$Variete 4 18499 4625 2.330 0.05598 .
fi$field.treatment:fi$periode 2 1525 763 0.384 0.68135
fi$Traitement:fi$Variete 12 37635 3136 1.580 0.09601 .
fi$Traitement:fi$periode 6 39565 6594 3.322 0.00347 **
fi$Variete:fi$periode 8 29769 3721 1.875 0.06337 .
Residuals 314 623273 1985
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> aov_id_an_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^2 + fi$annee + (1 |fi$ID))
> summary(aov_id_an_poids)
Df Sum Sq Mean Sq F value Pr(>F)
fi$field.treatment 1 133480 133480 69.967 2.04e-15 ***
fi$Traitement 3 472538 157513 82.564 < 2e-16 ***
fi$Variete 4 22593 5648 2.961 0.020062 *
fi$periode 2 79375 39688 20.803 3.30e-09 ***
fi$annee 2 28052 14026 7.352 0.000759 ***
fi$field.treatment:fi$Traitement 3 54575 18192 9.536 4.82e-06 ***
fi$field.treatment:fi$Variete 4 18499 4625 2.424 0.048135 *
fi$field.treatment:fi$periode 2 1525 763 0.400 0.670874
fi$Traitement:fi$Variete 12 37635 3136 1.644 0.078605 .
fi$Traitement:fi$periode 6 39565 6594 3.456 0.002545 **
fi$Variete:fi$periode 8 29769 3721 1.951 0.052301 .
Residuals 312 595222 1908

Voici les résultats si je veux avoir les interactions d'ordre 3

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
aov_id_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^3 + (1 |fi$ID))
summary(aov_id_poids)
 
aov_id_an_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^3 + fi$annee + (1 |fi$ID))
summary(aov_id_an_poids)

car j'aimerais aussi voir l'effet var*field.treatment.Traitement, ici quand je test les 2 formules, cela change légèrement et j'ai des différences aux 3 endroits en rouge:

> aov_id_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^3 + (1 |fi$ID))
> summary(aov_id_poids)
Df Sum Sq Mean Sq F value Pr(>F)
fi$field.treatment 1 133480 133480 66.608 1.36e-14 ***
fi$Traitement 3 472538 157513 78.601 < 2e-16 ***
fi$Variete 4 22593 5648 2.819 0.02566 *
fi$periode 2 79375 39688 19.805 9.69e-09 ***
fi$field.treatment:fi$Traitement 3 54575 18192 9.078 9.69e-06 ***
fi$field.treatment:fi$Variete 4 18499 4625 2.308 0.05846 .
fi$field.treatment:fi$periode 2 1525 763 0.381 0.68389
fi$Traitement:fi$Variete 12 37635 3136 1.565 0.10180
fi$Traitement:fi$periode 6 39565 6594 3.291 0.00387 **
fi$Variete:fi$periode 8 29769 3721 1.857 0.06707 .
fi$field.treatment:fi$Traitement:fi$Variete 12 30737 2561 1.278 0.23125
fi$field.treatment:fi$Traitement:fi$periode 6 7701 1283 0.640 0.69779
fi$field.treatment:fi$Variete:fi$periode 8 11865 1483 0.740 0.65602
fi$Traitement:fi$Variete:fi$periode 24 43928 1830 0.913 0.58388
Residuals 264 529042 2004
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> aov_id_an_poids<-aov(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^3 + fi$annee + (1 |fi$ID))
> summary(aov_id_an_poids)
Df Sum Sq Mean Sq F value Pr(>F)
fi$field.treatment 1 133480 133480 69.805 3.86e-15 ***
fi$Traitement 3 472538 157513 82.373 < 2e-16 ***
fi$Variete 4 22593 5648 2.954 0.020582 *
fi$periode 2 79375 39688 20.755 4.29e-09 ***
fi$annee 2 28052 14026 7.335 0.000795 ***
fi$field.treatment:fi$Traitement 3 54575 18192 9.514 5.50e-06 ***
fi$field.treatment:fi$Variete 4 18499 4625 2.419 0.049008 *
fi$field.treatment:fi$periode 2 1525 763 0.399 0.671557
fi$Traitement:fi$Variete 12 37635 3136 1.640 0.080733 .
fi$Traitement:fi$periode 6 39565 6594 3.448 0.002696 **
fi$Variete:fi$periode 8 29769 3721 1.946 0.053654 .
fi$field.treatment:fi$Traitement:fi$Variete 12 30737 2561 1.340 0.196061
fi$field.treatment:fi$Traitement:fi$periode 6 7701 1283 0.671 0.673024
fi$field.treatment:fi$Variete:fi$periode 8 11865 1483 0.776 0.624571
fi$Traitement:fi$Variete:fi$periode 24 43928 1830 0.957 0.523842
Residuals 262 500991 1912

Je ne suis pas sûre d'avoir compris quel est l’inconvénient à utiliser la formule 2 avec l'année pour être "tranquille" plutôt que la formule 1?
Je me demandais aussi ce que cela change si on écrit la formule avec des "*" au lieu des "+"?

Merci par avance.

(j'ai mis le fichier avec ID en PJ)

[faubry]

1)

J'ai mis les ID dans mon fichier que je viens de mettre en pièce-jointe, juste pour être sûre que c'est juste j'ai bien placé les ID? (exemple: pour (var 1, traitement P1, field.treatment=No et période = 3 mois) je lui met l'ID = 1 pour l'année 1 mais aussi pour l'année 2 et 3 puis pour (var 1, traitement P1, field.treatment=No et période = 5 mois) je met un ID 2, etc..

Non la répétition étant la période, il faut regrouper selon les périodes. Maintenant si tu considères que c'est le même échantillon qui est vu d'une année sur l'autre, il faut alors donner le même ID pour chaque année. Les questions à résoudre avant de coder les ID sont donc en dehors de celle de l'année :
- Est-ce que filed.treatment yes/no conduit à deux individus différentes ? Si non, il faut le même ID, si oui, il faut deux ID différents.
- Utilise-t-on le même échantillon pour tous les traitements (peu probable) ou fait-on subir des traitements différents à des échantillons différents de la même varieté. Dans le premier cas, c'est le même ID, dans le second, des ID différents.

Si l'étude a été correctement défini, le cahier des charges la décrivant doit répondre à ces questions.

2) C'est quoi cette utilisation de la fonction aov. On est dans des modèles mixtes, c'est donc lme4::lmer à utiliser.

3) Les deux formules que je t'ai données sont des formules emboitées, les résultats se comparent par la fonction anova comme suit :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
model.1 <- lmer( variable.analysee ~ (field.treatment + traitement + variété + période)^2 + annee +  (1 |ID), donnees)
model.2 <- lmer( variable.analysee ~ (field.treatment + traitement + variété + période)^2 +  (1 |ID), donnees)
anova( model.1, model.2)

4) Malgré cette utilisation bizarre de aov, on peut constater qu'avec des interactions d'ordre 2, on a un modèle qui tourne autour de 50 ddl et un résidus autour de 300. C'est donc déjà très limite question de robustesse. Avec les interactions d'ordre 3, le modèle a plus de 100 ddl et le résidus autour de 250, donc les résultats ne sont pas fiables. Il faut donc vraiment penser à simplifier le modèle de manière pertinente pour arriver à entre 10 et 15 (maximum) ddl, la discussion des résultats étant toujours le plus important dans une étude.

5) ~ A + B signifie que les variables n'interagissent pas ensemble, c'est-à-dire pour des facteurs que la valeurs des niveaux de A est indépendantes des niveaux de B et vice-versa,
~ A * B est équivalent à ~ A + B + A:B, c'est-à-dire qu'il existe des interactions entre les deux facteurs. Exemple illustratif : supposons qu'on analyse la taille d'un poisson en fonction de son sexe (F ou M) et du milieu où il vit (par exemple, lac ou rivière), alors
taille ~ sexe + milieu signifie que la taille peut dépendre du sexe, indépendamment du milieu ou du milieu indépendamment du sexe,
taille ~ sexe * milieu signifie qu'il existe un effet sexe dans un milieu donné ou un effet milieu pour un sexe donné, contrairement au modèle précédent.

Je ne sais pas dans quel cadre tu fais cette analyse mais, ne le prend surtout pas mal, j'ai l'impression qu'on t'a mis devant un problème qui dépasse largement tes connaissances en statistiques (et peut-être, celles de tes encadrants ?)

[Margotte02]

1)

Non la répétition étant la période, il faut regrouper selon les périodes. Maintenant si tu considères que c'est le même échantillon qui est vu d'une année sur l'autre, il faut alors donner le même ID pour chaque année. Les questions à résoudre avant de coder les ID sont donc en dehors de celle de l'année :
- Est-ce que filed.treatment yes/no conduit à deux individus différentes ? Si non, il faut le même ID, si oui, il faut deux ID différents.

Oui cela conduit à des individus différents

- Utilise-t-on le même échantillon pour tous les traitements (peu probable) ou fait-on subir des traitements différents à des échantillons différents de la même varieté. Dans le premier cas, c'est le même ID, dans le second, des ID différents.

je prends des plantes de la même variété que je sépare dans des unités différentes et que je traite avec des produits différents donc je pense que nous sommes plutôt dans le cas 2


Ok merci, je comprends mieux, alors selon ce raisonnement je dois mettre un ID pour un var*1traitement*1field.treatment*annee1 identique pour mes 3 périodes d'observations qui ne se font pas exactement sur le même échantillons puisque on échantillonnage est destructif de toute façon.

par ex.:
var1 * Traitement = control * field.treatment1 = no * année = 1 * période = 3 mois--> ID = 1
var1 * Traitement = control * field.treatment1 = no * année = 1 * période = 5 mois--> ID = 1
var1 * Traitement = control * field.treatment1 = no * année = 1 * période = 7 mois--> ID = 1

var1 * Traitement = control * field.treatment1 = no * année = 2 * période = 3 mois--> ID = 2
var1 * Traitement = control * field.treatment1 = no * année = 2 * période = 5 mois--> ID = 2
var1 * Traitement = control * field.treatment1 = no * année = 2 * période = 7 mois--> ID = 2

c'est juste comme ça?


2) Ok désolée j'ai refait avec lmer

3) j'ai refait tourner avec les bonnes formules et le bon fichier avec les ID placés comme expliqué dessus (le bon fichier est mtn "file2 en pj)

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
fi<-read.csv2("file2csv.csv", na.strings = "NA",header = TRUE)
dim(fi)
head(fi)
fi_2<-subset(fi, periode != "0month", select= c (ID, Variete, annee, periode , Traitement, field.treatment, poids_racines, taille_racines))
fi<-fi_2
fi$poids_racines<-as.numeric(fi$poids_racines)
fi$taille_racines<-as.numeric(fi$taille_racines)
library(lme4)
#avec année-poids
model.1<-lmer(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^2 + fi$annee + (1 |fi$ID))
 
##SANS année-poids
model.2<-lmer(fi$poids_racines~ (fi$field.treatment + fi$Traitement + fi$Variete + fi$periode)^2 + (1 |fi$ID))
#comparaison des 2 modèles pour voir si onn garde l'année
anova( model.1, model.2)

ça donne le résultat ci-dessous avec le message d'erreur en rouge, comment lire ce résultat?

> anova( model.1, model.2)
refitting model(s) with ML (instead of REML)
Data: NULL
Models:
model.2: fi$poids_racines ~ (fi$field.treatment + fi$Traitement + fi$Variete +
model.2: fi$periode)^2 + (1 | fi$ID)
model.1: fi$poids_racines ~ (fi$field.treatment + fi$Traitement + fi$Variete +
model.1: fi$periode)^2 + fi$annee + (1 | fi$ID)
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
model.2 48 3756.0 3942.5 -1830.0 3660.0
model.1 50 3750.3 3944.6 -1825.2 3650.3 9.6156 2 0.008166 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1