7 votos

Resultados completamente diferentes de lme () y lmer ()

He estado jugando con ambos nlme::lme y lme4::lmer. He instalado una simple aleatorio intercepta modelo de uso de la lme() y lmer(). Como se puede ver a continuación, llegué a resultados completamente diferentes de lmer() y lme(). Incluso los signos de los coeficientes son diferentes! Estoy haciendo algo mal? También he ajustado un modelo vacío con dos paquetes. En este caso, los resultados fueron prácticamente los mismos (resultados no mostrados). Habría que educar a mí entender este problema? A menos que cometí un error, creo que hay algo mal con el lme4 paquete.

     multi<-structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95, 
4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12, 
2.07, 1.05, 1, 1.07, 1.16, 1.11), y = c(3.2, 3.21, 3.256, 3.25, 
3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628, 
3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954, 
4, 3.994, 4), pid = 1:25, gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
5L, 5L)), class = "data.frame", row.names = c(NA, -25L))

#lme
> lme(y~x, random=~1|gid,data=multi,method="REML")
Linear mixed-effects model fit by REML
  Data: multi 
  Log-restricted-likelihood: 41.76745
  Fixed: y ~ x 
(Intercept)           x 
  4.1846756  -0.1928357 

#lmer

 lmer(y~x+(1|(gid)), data=multi, REML=T)
    Linear mixed model fit by REML ['lmerMod']
    Formula: y ~ x + (1 | (gid))
       Data: multi
    REML criterion at convergence: -78.4862
    Random effects:
     Groups   Name        Std.Dev.
     (gid)    (Intercept) 0.70325 
     Residual             0.02031 
    Number of obs: 25, groups:  (gid), 5
    Fixed Effects:
    (Intercept)            x  
         2.8152       0.2638 

13voto

user219012 Puntos 1

Como se señaló en esta respuesta, y también se menciona en uno de los comentarios, el problema parece ser de un máximo local. Para ver esto más claramente, he escrito a continuación un código simple para calcular la negativa de la log-verosimilitud de este modelo y hacer la optimización del uso de optim(). A partir de diferentes valores iniciales conduce a las dos soluciones diferentes:

# data
multi <- structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95, 
                              4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12, 
                              2.07, 1.05, 1, 1.07, 1.16, 1.11), 
                        y = c(3.2, 3.21, 3.256, 3.25, 
                              3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628, 
                              3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954, 
                              4, 3.994, 4), 
                        pid = 1:25, 
                        gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 
                                2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 
                                5L, 5L)), class = "data.frame", row.names = c(NA, -25L))

# function to calculate the negative log-likelihood of the random intercepts model
library("mvtnorm")
logLik <- function (thetas, y, X, id) {
    ncX <- ncol(X)
    betas <- thetas[seq_len(ncX)]
    sigma_b <- exp(thetas[ncX + 1])
    sigma <- exp(thetas[ncX + 2])
    eta <- c(X %*% betas)
    unq_id <- unique(id)
    n <- length(unq_id)
    lL <- numeric(n)
    for (i in seq_len(n)) {
        id_i <- id == unq_id[i]
        n_i <- sum(id_i)
        V_i <- matrix(sigma_b^2, n_i, n_i)
        diag(V_i) <- diag(V_i) + sigma^2
        lL[i] <- dmvnorm(y[id_i], mean = eta[id_i], sigma = V_i, log = TRUE)
    }
    - sum(lL, na.rm = TRUE)
}

# optimization using as initial values 0 for the fixed effects, 
# and 1 for the variance components 
opt <- optim(rep(0, 4), logLik, method = "BFGS", 
             y = multi

8voto

Ben Bolker Puntos 8729

Estoy de acuerdo con @DimitrisRizopoulos la respuesta, y tienen un par de puntos más para hacer.

  • Voy a empezar diciendo que soy infeliz que lmer no encontrar la mejor respuesta, aunque sospecho que esta situación es probablemente limitada a pequeñas, inusual (ver más abajo) conjuntos de datos. Una de las razones por las que lme puede hacer mejor es que encaja en el registro de-desviación de la escala, lo que puede hacer el mínimo de cerca de cero "más amplia".
  • Usted puede obtener lmer a replicar el lme de resultados mediante el establecimiento de un explícito, menor valor de inicio para la ampliación de la desviación estándar (start=...); sobre la base de las exploraciones a continuación, start=8 o menor valor debería funcionar bien. Para lo que vale, esto conducirá a una estimación de efectos aleatorios varianza de 0 (y un "singular ajuste" del mensaje, y una respuesta que es el equivalente a dejar fuera la de efectos aleatorios componente completamente y el uso de lm() ...)
  • En este caso en particular el uso de la "nloptwrap" optimizador no ayuda; de hecho, todos los optimizadores de que lmer puede utilizar, a partir de los valores iniciales predeterminados ($\theta$ (en escala de desviación estándar) = 1.0), encontrar el más mínimo local lejos de cero.
  • aquí está el código equivalente para el planteamiento lmer se utiliza para encontrar el valor inicial por defecto:
v0 <- with(multi,var(ave(y,gid)))  ## variance among group values  
v.e <- var(multi$y)-v0             ## residual var ~ total var - group variance
sqrt(v0/v.e)                       ## convert to scaled standard deviation

Esto conduce a un valor inicial de $\theta=10.8$.

  • Podemos ver de forma sistemática cómo los diferentes valores de partida dan resultados diferentes:
m0 <- lmer(y~x+(1|(gid)), data=multi, REML=TRUE)
tvec2 <- seq(0,20,length=51)
ff <- function(t0) getME(update(m0,start=t0),"theta")
v <- sapply(tvec2,ff)
plot(tvec2,v)
abline(v=10.8,col="red")

enter image description here

  • También podemos explícitamente visualizar el (logaritmo negativo-)la probabilidad de la superficie:
f <- as.function(m0)
tvec <- seq(0,100,length=101)
dvec <- sapply(tvec,f)
m3 <- update(m0,REML=FALSE)
f2 <- as.function(m3)
dvec2 <- sapply(tvec,f2)
par(las=1,bty="l")
matplot(tvec,cbind(dvec,dvec2),type="l",
        ylab="deviance/REMLcrit",
        xlab="scaled standard dev")
legend("bottomright",c("REML","ML"),
       col=1:2,lty=1:2)

enter image description here

Para lo que vale, el ML ajuste da $\theta=0$ más que el valor más alto.

  • Son estos datos artificiales? La izquierda gráfico siguiente muestra los datos por grupo; el derecho del gráfico muestra los valores con su grupo de medios resta. Casi no hay variación entre los 5 valores dentro de cada grupo ...

enter image description here

  • Si queremos simular los datos con las mismas propiedades (a partir de los coeficientes estimados), pero donde la variación es en realidad de Gauss, que no reciben el mismo tipo de multimodal de la superficie del todo:
multi_sim <- transform(multi,y=simulate(m0,seed=101)[[1]])
f3 <- as.function(update(m0,data=multi_sim))
dvec3 <- sapply(tvec,f3)
plot(tvec,dvec3,type="l")

enter image description here

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X