modèle mixte

Questions sur les fonctions statistiques de R

Modérateur : Groupe des modérateurs

Claire Bissery
Messages : 11
Enregistré le : 19 Juin 2006, 06:31

modèle mixte

Messagepar Claire Bissery » 19 Juin 2006, 09:56

Bonjour,

J’ai quelques problèmes avec la fonction lme sous R et sous Splus.

J’étudie la relation allométrique entre le diamètre des branches et la biomasse soutenue par ces branches chez le manguier.
Je voudrais regarder l’influence de la variété sur la relation entre le diamètre et la biomasse.
J’ai des données pour 6 variétés. Environ 30 données par variétés pour 4 variétés et 10 par variétés pour 2 variétés.

Sous R le modèle ne converge pas
lme1<-lme(lmstoter~cldbm,random=~cldbm|variete)
Error in lme.formula(lmstoter ~ cldbm, random = ~cldbm | variete) :
iteration limit reached without convergence (9)

avec lmstoter=log(biomasse totale de l’ensemble raméal )
cldbm=log(diamètre à la base de l’ensemble raméal ), centré par la moyenne des log(dbm)

Sous S il converge mais donne une corrélation entre l’intercept et cldbm de –1.

Je ne trouve pas ce qui dans les données peut poser problème.
Est ce quelqu’un pourrait m’aider à comprendre ce qui ne va pas ?
Merci par avance
Claire

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 19 Juin 2006, 11:55

Bonjour Claire,

Les modèles mixtes sont faits pour analyser des données répétées sur les mêmes individus (plusieurs obs sur les mêmes manguiers par exemple). Ce que tu décris ressemble à un plan factoriel classique analysable par un modèle linéaire (fonction lm):

Code : Tout sélectionner

lm(lmstoter ~ cldbm * variete)

Si tu as des données répétées sur les mêmes individus, il faut spécifier ce niveau d'agrégation dans la partie aléatoire du modèle. Par exemple si chaque arbre est repéré par une variable id:

Code : Tout sélectionner

lme(lmstoter ~ cldbm * variete,
    random = ~ 1 | id)

pour un modèle avec un seul effet aléatoire lié à l'intercept.

Je te conseille de lire le bouquin de P&B qui te permettra de te fixer les idées:

Pinheiro, J.C., Bates, D.M., 2000. Mixed-effect models in S and S-Plus. Springer, New York, 598 p.

Amicalement,

Renaud

Claire Bissery
Messages : 11
Enregistré le : 19 Juin 2006, 06:31

Messagepar Claire Bissery » 19 Juin 2006, 12:54

Merci pour votre réponse
Je veux regarder l’effet aléatoire de la variété sur la relation entre le diamètre et la biomasse.
Je recherche le modèle le plus adapté aux données (pente et intercept aléatoire, pente aléatoire, pas d’effet aléatoire de la variété).

Seulement, je ne comprend pas pourquoi le modèle pente et intercept aléatoire ne converge pas sous R et donne des valeurs de corrélation aberrantes sous S.
J’ai bien regardé mes données, et je n’ai pas trouvé d’explications.
J’ai utilisé des modèles mixtes dans d’autres études et je n’ai jamais rencontré ce type de problèmes.
Si vous avez le temps de regarder, je peux vous envoyer mes données.
Merci par avance
Claire

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 19 Juin 2006, 13:41

Pour moi, la variété n'est pas un effet aléatoire mais un effet fixe. Elle pourrait être considérée comme aléatoire si les variétés avaient été tirées aux sort dans une (grande) population de variétés, et que l'on souhaite faire une inférence sur l'effet de la variété sur la moyenne de pop. Je ne pense pas que vous soyez dans cette situation. Je veux bien jeter un coup d'oeil sur vos données, mais j'ai peu de temps.

Si vous persistez dans cette voie, essayez la fonction lmer dans le package lme4 (dans Matrix, en fait): est supposée être plus robuste que lme.

Renaud

Claire Bissery
Messages : 11
Enregistré le : 19 Juin 2006, 06:31

Messagepar Claire Bissery » 20 Juin 2006, 06:29

Je peux réaliser mon étude en considérant la variété comme un effet fixe, mais j’aimerai prendre en compte un effet aléatoire de la variété pour plusieurs raisons.
Premièrement j’aimerai voir si la variété a un effet aléatoire sur la relation allometrique générale entre le diamètre et la biomasse. Ensuite, les modèles mixte me permettent de réaliser une étude avec toutes mes données malgré le déséquilibre.
Mes variétés ne sont pas choisit aléatoirement mais il me semble qu’il suffit de vérifier la normalité de la distribution des estimations des paramètres aléatoires pour pouvoir utiliser des modèles mixtes.

La pente semble être peu différente entre les variétés mais j’aimerai le vérifier en comparant des modèles pente et intercept aléatoire et intercept aléatoire par un test de ratio de vraisemblance.
J’aimerai surtout comprendre pourquoi le modèle ne converge pas sous R et donne une corrélation de –1 sous S.
J’ai bien regardé mes données, j’ai essayé différents modèles sur des sous jeux de données et je n'ai pas trouvé d’où pouvait provenir ce problème. La fonction de lme4 ne converge pas non plus.

Si vous pouvez m’aider à comprendre pourquoi ce modèle ne marche pas.
Les données se trouvent sur ftp://ftp.cirad.fr/pub/group-r/groupe-r ... metrie.txt

Er : ensemble raméal
Dbm : diamètre à la base de l’er
Mstoter : masse sèche totale de l’er
Les er sont détruits pour être étudiés. Ils sont échantillonnées aléatoirement, en suivant une gamme de diamètre, sur une dizaine de plant génétiquement identiques.
Pour 4 varietes l’étude a été réalisée en 2004, 2005 et 2006.
2 variétés ont été rajoutées en 2006.
L’année a a priori peu d’effet sur la relation allometrique, et je pense l’étudier séparément si cela est confirmé pour la troisième année.

Claire

Renaud Lancelot
Messages : 2484
Enregistré le : 16 Déc 2004, 08:01
Contact :

Messagepar Renaud Lancelot » 20 Juin 2006, 07:03

Bonjour Claire,

Pour comprendre le pb, il suffit de faire un graphe de vos données. J'ai placé le tableau dans un répertoire personnel ==> code à adapter à vos besoins:

Code : Tout sélectionner

tab <- read.table("d:/analyses/travail/data/allometrie.txt", header = TRUE)

library(lattice)
xyplot(log(mstoter) ~ I(log(dbm) - log(mean(dbm))) | variete, data = tab,
  prepanel = function(x, y) prepanel.lmline(x, y), aspect = "xy",
  panel = function(x, y){
    panel.grid()
    panel.xyplot(x, y)
    panel.lmline(x, y)
    })


Ce graphe montre qu'il il a trés, très peu de variation de l'intercept et de la pente en fonction de la variété (avec la paramétrisation que vous avez retenue). En conséquence, on sait avant tout calcul que la fonction lme va ramer pour estimer les effets aléatoires liés à l'intercept et à la pente. Dans ces conditions, il est prudent de commencer par un modèle aléatoire le plus simple possible, i.e. sans la covariance entre pente et intercept:

Code : Tout sélectionner

> library(nlme)
> fm1 <- lme(log(mstoter) ~ I(log(dbm) - log(mean(dbm))),
+           random = list(variete = pdDiag(~ I(log(dbm) - log(mean(dbm))))), data = tab)
>
> summary(fm1)
Linear mixed-effects model fit by REML
 Data: tab
       AIC      BIC    logLik
  55.88947 70.97587 -22.94473

Random effects:
 Formula: ~I(log(dbm) - log(mean(dbm))) | variete
 Structure: Diagonal
        (Intercept) I(log(dbm) - log(mean(dbm)))  Residual
StdDev:   0.1241014                    0.0638946 0.2654799

Fixed effects: log(mstoter) ~ I(log(dbm) - log(mean(dbm)))
                                Value  Std.Error  DF  t-value p-value
(Intercept)                  3.524295 0.05642922 146 62.45515       0
I(log(dbm) - log(mean(dbm))) 2.667351 0.06305499 146 42.30199       0
 Correlation:
                             (Intr)
I(log(dbm) - log(mean(dbm))) 0.073

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.87820166 -0.70277749  0.03419046  0.71424256  2.40789680

Number of Observations: 153
Number of Groups: 6
> intervals(fm1)
Approximate 95% confidence intervals

 Fixed effects:
                                lower     est.    upper
(Intercept)                  3.412771 3.524295 3.635819
I(log(dbm) - log(mean(dbm))) 2.542733 2.667351 2.791970
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: variete
                                       lower      est.     upper
sd((Intercept))                  0.060140273 0.1241014 0.2560872
sd(I(log(dbm) - log(mean(dbm)))) 0.001709319 0.0638946 2.3883893

 Within-group standard error:
    lower      est.     upper
0.2364470 0.2654799 0.2980777


Le modèle converge sans pb mais remarquez que l'intervalle de confiance de la pente est énorme par rapport à la valeur de l'effet aléatoire. Cela est en accord avec le graphe. En conséquence, on peut s'attendre à des problèmes pour estimer le modèle plus complexe avec la covariance pente intercept:

Code : Tout sélectionner

> fm2 <- lme(log(mstoter) ~ I(log(dbm) - log(mean(dbm))),
+           random = list(variete = pdSymm(~ I(log(dbm) - log(mean(dbm))))), data = tab)
>
> summary(fm2)
Linear mixed-effects model fit by REML
 Data: tab
       AIC     BIC    logLik
  57.01002 75.1137 -22.50501

Random effects:
 Formula: ~I(log(dbm) - log(mean(dbm))) | variete
 Structure: General positive-definite
                             StdDev     Corr 
(Intercept)                  0.11926127 (Intr)
I(log(dbm) - log(mean(dbm))) 0.05874384 -1   
Residual                     0.26541168       

Fixed effects: log(mstoter) ~ I(log(dbm) - log(mean(dbm)))
                                Value  Std.Error  DF  t-value p-value
(Intercept)                  3.520635 0.05441790 146 64.69626       0
I(log(dbm) - log(mean(dbm))) 2.658298 0.06102595 146 43.56013       0
 Correlation:
                             (Intr)
I(log(dbm) - log(mean(dbm))) -0.303

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.87633430 -0.66212787  0.02015185  0.70313969  2.37963793

Number of Observations: 153
Number of Groups: 6
>
> intervals(fm2)
Erreur dans intervals.lme(fm2) : Cannot get confidence intervals on var-cov components: Non-positive definite approximate variance-covariance


Le modèle converge mais la matrice de variance covariance des effets aléatoires est "malade", indiquant un pb d'estimation. La corrélation de -1 entre pente et intercept indique que ces deux effets sont parfaitement corrélés ==> il serait suffisant de spécifier un seul des deux dans le modèle alaétoire, l'intercept par exemple:

Code : Tout sélectionner

> fm3 <- lme(log(mstoter) ~ I(log(dbm) - log(mean(dbm))),
+           random = ~ 1 | variete, data = tab)
>
> summary(fm3)
Linear mixed-effects model fit by REML
 Data: tab
       AIC      BIC    logLik
  53.98664 66.05576 -22.99332

Random effects:
 Formula: ~1 | variete
        (Intercept)  Residual
StdDev:   0.1264046 0.2662686

Fixed effects: log(mstoter) ~ I(log(dbm) - log(mean(dbm)))
                                Value  Std.Error  DF  t-value p-value
(Intercept)                  3.523360 0.05719728 146 61.60013       0
I(log(dbm) - log(mean(dbm))) 2.668937 0.05660827 146 47.14747       0
 Correlation:
                             (Intr)
I(log(dbm) - log(mean(dbm))) 0.077

Standardized Within-Group Residuals:
         Min           Q1          Med           Q3          Max
-2.848505650 -0.717101092  0.009730304  0.708112670  2.415555764

Number of Observations: 153
Number of Groups: 6
> intervals(fm3)
Approximate 95% confidence intervals

 Fixed effects:
                                lower     est.    upper
(Intercept)                  3.410318 3.523360 3.636401
I(log(dbm) - log(mean(dbm))) 2.557059 2.668937 2.780814
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: variete
                   lower      est.     upper
sd((Intercept)) 0.062141 0.1264046 0.2571271

 Within-group standard error:
    lower      est.     upper
0.2374758 0.2662686 0.2985524


NB: pour comparer les modèles avec le test du rapport des vraisemblances, il faut les ajuster avec la méthode du max de vraisemblances (method = "ML"), ce qui n'est pas la valeur par défaut.

Code : Tout sélectionner

> fm1 <- update(fm1, method = "ML")
> fm3 <- update(fm3, method = "ML")
>
> anova(fm1, fm3)
    Model df      AIC      BIC    logLik   Test      L.Ratio p-value
fm1     1  5 48.09807 63.25026 -19.04904                           
fm3     2  4 46.09845 58.22020 -19.04923 1 vs 2 0.0003823962  0.9844


Cela confirme que la pente n'est pas "significative": tous les indicateurs concordent (AIC, BIC et test).

Une autre manière d'aborder le problème est de faire une exploration graphique plus soigneuse des effets aléatoires: je vous conseille de lire le bouqin de P&B qui donne de nombreux détails et exemples. Voir en particulier la fct lmList. J'ai également écrit une publi sur cette question.

Cela dit, je persiste à penser qu'il n'est pas pertinent de considérer la variété comme un effet aléatoire dans cette étude...

Bien cordialement,

Renaud


Retourner vers « Archives : Fonctions statistiques »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité