lme () y lmer () dan resultados contradictorios

He estado trabajando con algunos datos que tienen algunos problemas con las mediciones repetidas. Al hacerlo, noté un comportamiento muy diferente entre lme()y lmer()usando mis datos de prueba y quiero saber por qué.

El conjunto de datos falsos que creé tiene medidas de altura y peso para 10 sujetos, tomados dos veces cada uno. Configuré los datos para que entre los sujetos hubiera una relación positiva entre la altura y el peso, pero una relación negativa entre las medidas repetidas dentro de cada individuo.

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

Aquí hay una gráfica de los datos, con líneas que conectan las dos medidas de cada individuo.

Así que ejecuté dos modelos, uno con lme()del nlmepaquete y otro con lmer()del lme4. En ambos casos, realicé una regresión de peso contra altura con un efecto aleatorio de ID para controlar las mediciones repetidas de cada individuo.

library(nlme)
Mlme=lme(Height~Weight,random=~1|ID,data=DF)
library(lme4)
Mlmer=lmer(Height~Weight+(1|ID),data=DF)

Estos dos modelos a menudo (aunque no siempre dependiendo de la semilla) generaron resultados completamente diferentes. He visto dónde generan estimaciones de varianza ligeramente diferentes, calculan diferentes grados de libertad, etc., pero aquí los coeficientes están en direcciones opuestas.

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

Para ilustrar visualmente, modelar con lme()

Y modelar con lmer()

¿Por qué estos modelos divergen tanto?

tl; dr si cambia el optimizador a "nloptwrap", creo que evitará estos problemas (probablemente).

¡Felicitaciones, ha encontrado uno de los ejemplos más simples de óptimos múltiples en un problema de estimación estadística! El parámetro que se lme4usa internamente (por lo tanto, conveniente para la ilustración) es la desviación estándar escalada de los efectos aleatorios, es decir, el std dev entre grupos dividido por el std dev residual.

Extraiga estos valores para el original lmey los lmerajustes:

(sd1 <- sqrt(getVarCov(Mlme)[[1]])/sigma(Mlme))
## 2.332469
(sd2 <- getME(Mlmer,"theta")) ## 14.48926

Vuelva a instalar con otro optimizador (probablemente será el valor predeterminado en la próxima versión de lme4):

Mlmer2 <- update(Mlmer,
  control=lmerControl(optimizer="nloptwrap"))
sd3 <- getME(Mlmer2,"theta")   ## 2.33247

Partidos lme... veamos qué está pasando. La función de desviación (-2 * probabilidad de registro), o en este caso la función de criterio REML análoga, para LMM con un solo efecto aleatorio toma solo un argumento, porque los parámetros de efecto fijo se perfilan ; se pueden calcular automáticamente para un valor dado de la desviación estándar de RE.

ff <- as.function(Mlmer)
tvec <- seq(0,20,length=101)
Lvec <- sapply(tvec,ff)
png("CV38425.png")
par(bty="l",las=1)
plot(tvec,Lvec,type="l",
     ylab="REML criterion",
     xlab="scaled random effects standard deviation")
abline(v=1,lty=2)
points(sd1,ff(sd1),pch=16,col=1)
points(sd2,ff(sd2),pch=16,col=2)
points(sd3,ff(sd3),pch=1,col=4)
dev.off()

Continué obsesionarse más sobre esto y corriendo los ajustes para las semillas aleatorias de 1 a 1000, encajando lme, lmery lmer+ nloptwrap para cada caso. Aquí están los números de 1000 donde un método dado obtiene respuestas que son al menos 0.001 unidades de desviación peor que otro ...

          lme.dev lmer.dev lmer2.dev
lme.dev         0       64        61
lmer.dev      369        0       326
lmer2.dev      43        3         0

En otras palabras, (1) no existe un método que siempre funcione mejor; (2) lmercon el optimizador predeterminado es peor (falla aproximadamente 1/3 del tiempo); (3) lmercon "nloptwrap" es mejor (peor que el lme4% de las veces, rara vez peor que lmer).

Para ser un poco tranquilizador, creo que es probable que esta situación sea peor para casos pequeños y mal especificados (es decir, el error residual aquí es uniforme en lugar de Normal). Sin embargo, sería interesante explorar esto de manera más sistemática ...