lme4 message d'erreur

**Margotte02** · 17/07/2019, 15h07

j'ai aussi encore une autre question si vous permettez:

avec quelles commandes R / où voyez vous dans les résultats la robustesse et fiabilité du modèle? par exemple pour le modèle avec l'année vous disiez que la variance était nulle alors qu'on ne pouvait pas l'utiliser.

**faubry** · 17/07/2019, 15h31

il n'y a pas eu de plan statistique fait à priori (d'où mes problèmes), l'idée de mon superviseur était vraiment simplement de voir l'éffet des traitements, field.treat. et l'effet variété sur le poids et la taille des racines afin de trouver le ou les traitements les plus efficaces.

Existe-t-il un intérêt à comparer les variétés entre elles ? Si non, pourquoi ne pas faire une analyse variété par variété ?

il n'y a pas eu de plan statistique fait à priori (d'où mes problèmes), l'idée de mon superviseur était vraiment simplement de voir l'éffet des traitements, field.treat. et l'effet variété sur le poids et la taille des racines

Il avait bien une idée de comment il allait procéder.

sur le poids et la taille des racines afin de trouver le ou les traitements les plus efficaces

Et si le traitement qui se révèle le plus efficace pour l'une des caractéristiques est le moins efficace pour l'autre ? Comment définir un compromis ?

- Différence entre les formules anova : les deux premières utilisent une correction des ddl du modèle mixte de type Satterthwaite, la dernière utilise les formules du package lme4. Type="III", calcule les significativités marginales, c'est-à-dire pour explique rapidement, les significativités des interactions puis de l'apport du facteur quand on a retiré l'effet de l'interaction, le type="II", c'est l'apport du facteur en entier.

- Le problème robustesse vs. nombre de ddl est un problème général des statistiques. Quand à la variance nulle pour un facteur aléatoire, c'est summary qui le donne.

**Margotte02** · 17/07/2019, 15h47

Par exemple ici, dans la formule que vous aviez donné avec l'année et la période, quand je travaille sur le file3 (un ID par unité donc un ID pour les 3 périodes et les 3 années), cette fois la formule fonctionne et la variance ne semble pas nulle, quand je fais "summary" j'obtiens les résultats ci-dessous:

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("file3csv.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)
fi <- droplevels( subset( fi, periode != "0month")) 
library(lme4)
> model.1<-lmer(fi$poids_racines~ (field.treatment + Traitement + Variete + periode)^2 + annee + (1 |ID), fi)
> summary(model.1)

> summary(model.1)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: fi$poids_racines ~ (field.treatment + Traitement + Variete + periode)^2 + annee + (1 | ID)
Data: fi

REML criterion at convergence: 3376.6

Scaled residuals:
Min 1Q Median 3Q Max
-3.12008 -0.52707 -0.02464 0.49347 2.98721

Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 75.54 8.691
Residual 1881.61 43.378
Number of obs: 360, groups: ID, 40

Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 84.0194 17.4287 27.5929 4.821 4.68e-05 ***
field.treatmentYes -60.1389 16.4164 16.7896 -3.663 0.001959 **
TraitementP1 -15.8778 20.6193 18.5517 -0.770 0.450968
TraitementP2 -91.8444 20.6193 18.5517 -4.454 0.000287 ***
TraitementP3 -36.4444 20.6193 18.5517 -1.767 0.093590 .
Varietevar2 48.6528 21.4544 20.0190 2.268 0.034568 *
Varietevar3 64.7778 21.4544 20.0190 3.019 0.006769 **
Varietevar4 3.6944 21.4544 20.0190 0.172 0.865009
Varietevar5 33.5694 21.4544 20.0190 1.565 0.133327
periode5months 53.8500 16.8001 300.0000 3.205 0.001494 **
periode7months 78.6583 16.8001 300.0000 4.682 4.31e-06 ***
anneean2 0.8083 5.6000 300.0000 0.144 0.885325
anneean3 -18.3083 5.6000 300.0000 -3.269 0.001203 **
field.treatmentYes:TraitementP1 31.0889 15.0892 12.0000 2.060 0.061734 .
field.treatmentYes:TraitementP2 69.0000 15.0892 12.0000 4.573 0.000640 ***
field.treatmentYes:TraitementP3 26.3111 15.0892 12.0000 1.744 0.106748
field.treatmentYes:Varietevar2 -18.0278 16.8703 12.0000 -1.069 0.306269
field.treatmentYes:Varietevar3 -31.6389 16.8703 12.0000 -1.875 0.085268 .
field.treatmentYes:Varietevar4 6.6944 16.8703 12.0000 0.397 0.698467
field.treatmentYes:Varietevar5 -22.9167 16.8703 12.0000 -1.358 0.199324
field.treatmentYes:periode5months 0.6000 11.2000 300.0000 0.054 0.957312
field.treatmentYes:periode7months 9.0167 11.2000 300.0000 0.805 0.421425
TraitementP1:Varietevar2 -42.8889 23.8581 12.0000 -1.798 0.097420 .
TraitementP2:Varietevar2 -19.5000 23.8581 12.0000 -0.817 0.429669
TraitementP3:Varietevar2 -1.6667 23.8581 12.0000 -0.070 0.945458
TraitementP1:Varietevar3 -7.3333 23.8581 12.0000 -0.307 0.763829
TraitementP2:Varietevar3 -18.7222 23.8581 12.0000 -0.785 0.447827
TraitementP3:Varietevar3 -15.4444 23.8581 12.0000 -0.647 0.529598
TraitementP1:Varietevar4 7.3333 23.8581 12.0000 0.307 0.763829
TraitementP2:Varietevar4 3.1667 23.8581 12.0000 0.133 0.896607
TraitementP3:Varietevar4 -8.3333 23.8581 12.0000 -0.349 0.732927
TraitementP1:Varietevar5 -47.1111 23.8581 12.0000 -1.975 0.071774 .
TraitementP2:Varietevar5 -5.8889 23.8581 12.0000 -0.247 0.809212
TraitementP3:Varietevar5 -8.9444 23.8581 12.0000 -0.375 0.714276
TraitementP1:periode5months 2.3667 15.8393 300.0000 0.149 0.881324
TraitementP2:periode5months -35.5333 15.8393 300.0000 -2.243 0.025602 *
TraitementP3:periode5months -24.1000 15.8393 300.0000 -1.522 0.129179
TraitementP1:periode7months 11.8000 15.8393 300.0000 0.745 0.456864
TraitementP2:periode7months -51.7667 15.8393 300.0000 -3.268 0.001208 **
TraitementP3:periode7months -7.2000 15.8393 300.0000 -0.455 0.649750
Varietevar2:periode5months -19.1667 17.7088 300.0000 -1.082 0.279978
Varietevar3:periode5months -34.7083 17.7088 300.0000 -1.960 0.050928 .
Varietevar4:periode5months -17.1667 17.7088 300.0000 -0.969 0.333134
Varietevar5:periode5months -28.0417 17.7088 300.0000 -1.583 0.114364
Varietevar2:periode7months -58.0833 17.7088 300.0000 -3.280 0.001161 **
Varietevar3:periode7months -53.5000 17.7088 300.0000 -3.021 0.002736 **
Varietevar4:periode7months -27.9583 17.7088 300.0000 -1.579 0.115440
Varietevar5:periode7months -35.7917 17.7088 300.0000 -2.021 0.044155 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation matrix not shown by default, as p = 48 > 12.
Use print(x, correlation=TRUE) or
vcov(x) if you need it

et quand j'appelle l'anova:

Code :

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

anova(model.1)

ça donnne ceci:

> anova(model.1)
Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
field.treatment 98053 98053 1 12 52.1110 1.059e-05 ***
Traitement 347122 115707 3 12 61.4936 1.482e-07 ***
Variete 16597 4149 4 12 2.2051 0.1298296
periode 79375 39688 2 300 21.0923 2.686e-09 ***
annee 28052 14026 2 300 7.4542 0.0006927 ***
field.treatment:Traitement 40090 13363 3 12 7.1021 0.0053304 **
field.treatment:Variete 13589 3397 4 12 1.8056 0.1926324
field.treatment:periode 1525 763 2 300 0.4052 0.6671821
Traitement:Variete 27647 2304 12 12 1.2244 0.3657338
Traitement:periode 39565 6594 6 300 3.5045 0.0022952 **
Variete:periode 29769 3721 8 300 1.9777 0.0489445 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Existe-t-il un intérêt à comparer les variétés entre elles ? Si non, pourquoi ne pas faire une analyse variété par variété ?

le but premier est de comparer l’efficacité des traitements et field.treatment et de s'assurer que ces traitements soient efficaces sur toutes les variétés. Comparer les variétés entre elles n'est pas le but premier de l'étude même si c'est intéressant de voir les différences entre variétés, nous pourrions en effet faire une analyse séparée pour chaque variété, à ce moment là je dois creer un subset variété par variété et rappeler la formule ci-dessus sans le "var"?

**faubry** · 17/07/2019, 17h29

le but premier est de comparer l’efficacité des traitements et field.treatment et de s'assurer que ces traitements soient efficaces sur toutes les variétés. Comparer les variétés entre elles n'est pas le but premier de l'étude même si c'est intéressant de voir les différences entre variétés, nous pourrions en effet faire une analyse séparée pour chaque variété, à ce moment là je dois creer un subset variété par variété et rappeler la formule ci-dessus sans le "var"?

Je reviens toujours à ce que j'ai dit sur le nombre de ddl, moins il y en a, plus l'analyse est robuste. C'est un premier point. Ensuite, il faut d'abord passer par anova avant d'essayer d'interpréter les résultats de summary. Que constate-t-on ? La très faible influence de la variété. Il semble d'après anova que la seule interaction significative soit avec période, ce qui signifie que l'effet des périodes est différent d'une variété à l'autre mais pas le reste.

L'utilisation du facteur année comme effet fixe me pose problème. Je m'en suis déjà expliqué. Comment interpréter les variations d'une année sur l'autre alors que les années ne sont pas caractérisées : il a fait plus chaud, moins chaud, il a plus plu ? Ou ne sont-ce que des variations naturelles aléatoires et cet effet fixe de l'année un artefact ? Sauf justification solide qui invalide mon opinion, l'année est juste utilisée pour répéter la mesure en s'affranchissant des variations possibles d'une année sur l'autre et en espérant ainsi stabiliser les mesures. On est proche de ce qu'on appelle des moyennes mobiles, faire une moyenne sur plusieurs mesures sachant que chacune des mesures est affectée d'un bruit aléatoire.

Ma suggestion est de simplifier l'étude statistique en :
1) faisant une étude variété par variété,
2) utilisant comme valeurs mesurées (longueur et poids), les moyennes des trois années,
3) de ce fait en utilisant exclusivement des modèles fixes (fonction lm) avec comme formules mesures.moyennes ~ (traitement + field.traitement + periode)^2, avec possibilité de les simplifier (voir fonctions step et stepAIC),
4) éventuellement, utilisant comme approche de type Manova pour analyser conjointement la taille et le poids.

Dernière chose, les fonctions read.table, read.csv et ead.csv2 ont deux argument dec et sep qui permettent de spécifier le codage décimale (dec="." pour le codage anglo-saxon du point décimal et dec="," pour le codage par virgule) et sep pour le séparateur des colonnes. Il faut les spécifier dans ton appel. Par exemple si tes données utilisent la virgule décimale, alors il faut écrire dec=",", read.csv2 renvoie alors automatiquement des nombres.

**Margotte02** · 18/07/2019, 05h18

Merci beaucoup pour votre suggestion. Je vais simplifier le modèle comme vous le suggérez.

ou voyez-vous le nombre de dll?

Pour "analyser" la fonction lm, je peux aussi faire comme ceci?

Code :

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

1
2
lm1 <- lm(mesures.moyennes ~ (traitement + field.traitement + periode)^2)
anova(lm1)

Que considérez vous comme modèle avec une bonne variance ou bonne ddl? j'ai bien compris que le nombre de ddl doit être bas mais à partir de quand le considère t on élevé?

Merci encore.

**faubry** · 18/07/2019, 08h52

ou voyez-vous le nombre de dll?

Le nombre de ddl d'une variable dépendante est donné par la colonne Df ou NumDF (selon les cas) de la fonction anova ou Anova. Le nombre de ddl de l'analyse, c'est-à-dire le nombre de coefficients estimés, est la somme pour toutes les variables dépendantes, excepté le nombre de ddl du résidus quand il est donné par anova (dans les modèle fixe, ligne Residuals). C'est aussi le nombre de lignes de la partie fixed effects de summary moins la ligne (intercept).

Pour "analyser" la fonction lm, je peux aussi faire comme ceci?

Question de vocabulaire. On n'analyse pas la fonction lm, on utilise la fonction lm qui est la fonction par défaut de R pour analyser un modèle linéaire à effets fixes. Son équivalent pour un modèle linéaire mixte est lmer.

Code :

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

lm1 <- lm(mesures.moyennes ~ (traitement + field.traitement + periode)^2)

Ne pas oublier l'argument des données. Deux solutions, soit utiliser toutes les données, c'est-à-dire faire une analyse toutes variétés confondues, soit faire une analyse par variété.

Que considérez vous comme modèle avec une bonne variance ou bonne ddl? j'ai bien compris que le nombre de ddl doit être bas mais à partir de quand le considère t on élevé?

Une bonne variance ne peut pas être définie puisque cela dépend des données mais il existe des formules qui permettent d'estimer le nombre d'individus en fonction de connaissances a priori sur la variance du phénomène analysé (en utilisant généralement la littérature sur le sujet ou à partir d'une première analyse exploratoire) et, si le type de question que l'on se pose le demande, la taille de l'effet.
Le bon nombre de ddl ? C'est assez empirique. Pour un modèle fixe, je considère que le rapport entre le nombre d'individus et le nombre de ddl doit être au minimum entre 10 et 20. Pour un modèle mixte, le nombre réel de ddl doit aussi tenir compte des facteurs aléatoires or ce nombre de ddl fait l'objet de débat, d'où les corrections proposées par lmerTest. La littérature estime que le nombre de ddl d'un modèle mixte est compris entre celui du modèle fixe plus le nombre de paramètres de la matrice de variance de covariance du modèle aléatoire et le nombre de ddl du modèle fixe + le nombre de ddl du modèle aléatoire. Empiriquement, je dirai que le nombre d'individus doit être compris entre ces deux nombres en respectant la règle empirique d'un facteur 10 à 20 minimum.

Précisions : 1) Pour un facteur à N niveaux, le nombre de ddl est N - 1 puisque, pour faire rapide, les coefficients estimés sont les différences à la moyenne. Donc puisque la somme des N coefficients doit être nulle, il suffit d'en connaître N - 1 pour calculer le dernier.
2) Pourquoi le rapport 10 à 20 ? Simplement, parce qu'on considère qu'une estimation correcte de la variance d'un phénomène gaussien nécessite au moins 12 échantillons et qu'elle devient vraiment robuste à partir de 30. Ensuite, selon les disciplines, il n'est pas toujours possible de trop multiplier le nombre d'individus, souvent pour une question de coût.

**Margotte02** · 18/07/2019, 09h42

merci pour toutes ces explications, cela m'est utile et me sera très utile pour mes autres analyses.

Question de vocabulaire. On n'analyse pas la fonction lm, on utilise la fonction lm qui est la fonction par défaut de R pour analyser un modèle linéaire à effets fixes. Son équivalent pour un modèle linéaire mixte est lmer.

Donc j'utilise lm pour un modèle à effets fixes donc ici sans la variété et lmer pour un modèle mixte donc avec la variété? Dans les 2 cas, l'année sert de répétition mais pour l'inclure dans le modèle je dois m'assurer que le ddl et la variance restent convenables, c'est juste?

Si j'utilise lm et que j'analyse variété par variété avec cette formule, est-ce que je peux ajouter l'année ci-dessous? Ici il n'y a donc plus d'utilité à avoir l'ID?

Code :

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

1
2
lm1 <- lm(mesures.moyennes ~ (traitement + field.traitement + periode)^2)
anova(lm1)

**faubry** · 18/07/2019, 10h31

Donc j'utilise lm pour un modèle à effets fixes donc ici sans la variété et lmer pour un modèle mixte donc avec la variété?

Ce n'est absolument pas ce que j'ai voulu dire. Je dis seulement que lm est la fonction qui estime un modèle linéaire à effets fixes et lmer un modèle linéaire mixte.

Dans les 2 cas, l'année sert de répétition mais pour l'inclure dans le modèle je dois m'assurer que le ddl et la variance restent convenables, c'est juste?

SI tu calcule la moyenne sur les tois années, ce que je te propose, alors il n'y a plus lieu d'inclure l'année puisque ce facteur disparaît.

Dans ton script, quelles sont tes données ?

Enfin, je ne peux que te conseiller de chercher des introductions aux analyses statistiques, ça ne manque pas, et à R sinon tu vas être en grande difficulté pour la suite de ton travail. Pour R, il y a pas mal de tutoriels sur le site de R (https://cran.r-project.org) et sur ce site.

**Margotte02** · 19/07/2019, 06h48

SI tu calcule la moyenne sur les tois années, ce que je te propose, alors il n'y a plus lieu d'inclure l'année puisque ce facteur disparaît

je dois la calculer avant de lancer la fonction lm ou la fonction lm le fait automatiquement?

Enfin, je ne peux que te conseiller de chercher des introductions aux analyses statistiques, ça ne manque pas, et à R sinon tu vas être en grande difficulté pour la suite de ton travail. Pour R, il y a pas mal de tutoriels sur le site de R (https://cran.r-project.org) et sur ce site.

Oui merci pour tout, je vais lire des tutoriels. Ici les gens utilisent statistica mais j'aimerais apprendre plutôt R.

**faubry** · 19/07/2019, 08h12

je dois la calculer avant de lancer la fonction lm ou la fonction lm le fait automatiquement?

Tu peut passer ton tableau initial à lm. Il ne calcule pas vraiment automatiquement la moyenne mais c'est tout comme.

Personnellement, je vois trois avantages principaux à R par rapport à Statitistica que j'ai aussi utilisé. Ce sont ces avantages qui m'ont fait d'ailleurs passer à R :
1) R est gratuit et open-source,
2) Statistica est presse bouton, c'est-à-dire qu'on utilise généralement les procédures définies par défaut sans se poser de question ce qui peut amener à des erreurs d'analyse tandis que R oblige à réfléchir plus ; pour moi c'est un avantage, pour certains c'est un défaut de R,
3) la base de R et certains package sont certifiés par la FDA pour les essais cliniques.

**Margotte02** · 19/07/2019, 08h53

Oui vous avez raison, je vais m'entrainer sur R pour pouvoir l'utiliser facilement pour la suite.
Encore merci pour toute votre aide et le temps passé sur mon étude, c'est vraiment gentil.

lme4 message d'erreur

R

Discussions similaires

Partager

Partager