J'ai travaillé avec certaines données qui ont des problèmes avec les mesures répétées. Ce faisant, j'ai remarqué un comportement très différent entre lme()
et en lmer()
utilisant mes données de test et je veux savoir pourquoi.
Le faux ensemble de données que j'ai créé contient des mesures de taille et de poids pour 10 sujets, prises deux fois chacune. J'ai configuré les données de sorte qu'entre les sujets il y aurait une relation positive entre la taille et le poids, mais une relation négative entre les mesures répétées au sein de chaque individu.
set.seed(21)
Height=1:10; Height=Height+runif(10,min=0,max=3) #First height measurement
Weight=1:10; Weight=Weight+runif(10,min=0,max=3) #First weight measurement
Height2=Height+runif(10,min=0,max=1) #second height measurement
Weight2=Weight-runif(10,min=0,max=1) #second weight measurement
Height=c(Height,Height2) #combine height and wight measurements
Weight=c(Weight,Weight2)
DF=data.frame(Height,Weight) #generate data frame
DF$ID=as.factor(rep(1:10,2)) #add subject ID
DF$Number=as.factor(c(rep(1,10),rep(2,10))) #differentiate between first and second measurement
Voici un tracé des données, avec des lignes reliant les deux mesures de chaque individu.
J'ai donc exécuté deux modèles, un avec lme()
le nlme
package et un avec lmer()
de lme4
. Dans les deux cas, j'ai effectué une régression du poids en fonction de la taille avec un effet aléatoire d'ID pour contrôler les mesures répétées de chaque individu.
library(nlme)
Mlme=lme(Height~Weight,random=~1|ID,data=DF)
library(lme4)
Mlmer=lmer(Height~Weight+(1|ID),data=DF)
Ces deux modèles ont souvent (mais pas toujours dépendant de la graine) généré des résultats complètement différents. J'ai vu où ils génèrent des estimations de variance légèrement différentes, calculent différents degrés de liberté, etc., mais ici les coefficients sont dans des directions opposées.
coef(Mlme)
# (Intercept) Weight
#1 1.57102183 0.7477639
#2 -0.08765784 0.7477639
#3 3.33128509 0.7477639
#4 1.09639883 0.7477639
#5 4.08969282 0.7477639
#6 4.48649982 0.7477639
#7 1.37824171 0.7477639
#8 2.54690995 0.7477639
#9 4.43051687 0.7477639
#10 4.04812243 0.7477639
coef(Mlmer)
# (Intercept) Weight
#1 4.689264 -0.516824
#2 5.427231 -0.516824
#3 6.943274 -0.516824
#4 7.832617 -0.516824
#5 10.656164 -0.516824
#6 12.256954 -0.516824
#7 11.963619 -0.516824
#8 13.304242 -0.516824
#9 17.637284 -0.516824
#10 18.883624 -0.516824
Pour illustrer visuellement, modélisez avec lme()
Et modèle avec lmer()
Pourquoi ces modèles divergent-ils autant?
la source
Réponses:
tl; dr si vous changez l'optimiseur en "nloptwrap", je pense que cela évitera ces problèmes (probablement).
Félicitations, vous avez trouvé l'un des exemples les plus simples d'optima multiples dans un problème d'estimation statistique! Le paramètre qui
lme4
utilise en interne (donc pratique pour l'illustration) est l' écart-type échelonné des effets aléatoires, c'est-à-dire l'écart std parmi les groupes divisé par l'écart std résiduel.Extraire ces valeurs pour l'original
lme
etlmer
s'adapte:Réinitialisez avec un autre optimiseur (ce sera probablement la valeur par défaut dans la prochaine version de
lme4
):Matchs
lme
... voyons ce qui se passe. La fonction de déviance (-2 * log vraisemblance), ou dans ce cas la fonction de critère REML analogue, pour les LMM avec un seul effet aléatoire ne prend qu'un seul argument, car les paramètres à effet fixe sont profilés ; ils peuvent être calculés automatiquement pour une valeur donnée de l'écart type RE.Je continuais à obsèdent encore en cours et couru les ajustements pour les graines au hasard de 1 à 1000, montage
lme
,lmer
etlmer
+ nloptwrap pour chaque cas. Voici les nombres sur 1000 où une méthode donnée obtient des réponses qui sont au moins 0,001 unités de déviance pire qu'une autre ...En d'autres termes, (1) il n'y a pas de méthode qui fonctionne toujours mieux; (2)
lmer
avec l'optimiseur par défaut est le pire (échoue environ 1/3 du temps); (3)lmer
avec "nloptwrap" est le meilleur (pire quelme
4% du temps, rarement pire quelmer
).Pour être un peu rassurant, je pense que cette situation risque d'être pire pour les petits cas mal spécifiés (c'est-à-dire que l'erreur résiduelle ici est uniforme plutôt que normale). Il serait intéressant d'explorer cela plus systématiquement cependant ...
la source