2011-10-27 27 views
9

Para IGF datos de nlme biblioteca, que estoy recibiendo este mensaje de error:Error con nlme

lme(conc ~ 1, data=IGF, random=~age|Lot) 
Error in lme.formula(conc ~ 1, data = IGF, random = ~age | Lot) : 
    nlminb problem, convergence error code = 1 
    message = iteration limit reached without convergence (10) 

Pero todo está bien con este código

lme(conc ~ age, data=IGF) 
Linear mixed-effects model fit by REML 
    Data: IGF 
    Log-restricted-likelihood: -297.1831 
    Fixed: conc ~ age 
(Intercept)   age 
5.374974367 -0.002535021 

Random effects: 
Formula: ~age | Lot 
Structure: General positive-definite 
      StdDev  Corr 
(Intercept) 0.082512196 (Intr) 
age   0.008092173 -1  
Residual 0.820627711  

Number of Observations: 237 
Number of Groups: 10 

Como IGF es groupedData, por lo tanto los códigos Son identicos. Estoy confundido por qué el primer código produce un error. Gracias por tu tiempo y ayuda.

+0

Eché un vistazo rápido a esto y nada salta a la vista. Puede que tengas mejor suerte en la lista de correo 'r-sig-mixed-models', que tiene una concentración mucho mayor de personas familiarizadas con este paquete ... –

+0

¿Has intentado aumentar los límites de iteración en el primer ejemplo? Ver '? LmeControl'. –

+0

Ver respuesta y comentarios a continuación. Su primer modelo no tiene la edad como un efecto fijo, ni las restricciones de efectos aleatorios que tiene el segundo modelo. –

Respuesta

5

Si traza los datos, puede ver que no hay ningún efecto de age, por lo que parece extraño intentar ajustar un efecto aleatorio de age a pesar de esto. No es de extrañar que no sea convergente.

library(nlme) 
library(ggplot2) 

dev.new(width=6, height=3) 
qplot(age, conc, data=IGF) + facet_wrap(~Lot, nrow=2) + geom_smooth(method='lm') 

enter image description here

Creo que lo que quiere hacer es modelo de un efecto aleatorio de Lot en la intersección. Podemos tratar incluyendo age como efecto fijo, pero vamos a ver que no es importante y puede ser expulsado:

> summary(lme(conc ~ 1 + age, data=IGF, random=~1|Lot)) 
Linear mixed-effects model fit by REML 
Data: IGF 
     AIC  BIC logLik 
    604.8711 618.7094 -298.4355 

Random effects: 
Formula: ~1 | Lot 
     (Intercept) Residual 
StdDev: 0.07153912 0.829998 

Fixed effects: conc ~ 1 + age 
       Value Std.Error DF t-value p-value 
(Intercept) 5.354435 0.10619982 226 50.41849 0.0000 
age   -0.000817 0.00396984 226 -0.20587 0.8371 
Correlation: 
    (Intr) 
age -0.828 

Standardized Within-Group Residuals: 
     Min   Q1   Med   Q3   Max 
-5.46774548 -0.43073893 -0.01519143 0.30336310 5.28952876 

Number of Observations: 237 
Number of Groups: 10 
+1

Su análisis sin duda responde a la pregunta de qué está pasando en los datos, pero todavía hay una pregunta interesante sobre cuáles son las diferencias en los modelos que realmente están equipados. Al observar los resultados del modelo exitoso anterior, puede ver que * no * se ajusta a un efecto aleatorio de la edad (aunque existe una correlación perfecta con la variación de intersección entre lotes, lo que indica que el modelo está sobreajustado ...) –

+0

modelo en la publicación de OP que funciona ajustando una pendiente de 'edad', con un efecto aleatorio de 'Lote' en esa pendiente. Eso es bueno si los datos lo respaldan.Para un buen ejemplo donde ese es el caso, haga 'lme (height ~ age, data = Oxboys, random = ~ 1 + age | Subject)'. Este también es el ejemplo en §4.9.3 en el libro 'ggplot2'. El primer modelo en la publicación de OP, que no funciona, tiene un efecto aleatorio para algo que no está especificado en la estructura de efectos fijos. Ni siquiera creo que tenga sentido. –

+1

Sí, pero esto también falla: 'lme (conc-age, data = IGF, random = ~ age | Lot)', lo que parecería ser un modelo idéntico. (No estoy * demasiado * inclinado a gastar mucho esfuerzo para seguir con esto, aunque estoy un poco curioso acerca de la respuesta, porque parece pertenecer a la categoría de: "Doctor, me duele cuando hago esto. . "" Bueno, entonces no hagas eso ... "') –

3

encuentro el otro, la respuesta más viejo aquí insatisfactoria. Yo distingo entre los casos donde, estadísticamente, la edad no tiene impacto y, a la inversa, encontramos un error de cálculo. Personalmente, he cometido errores de carrera al combinar estos dos casos. R ha señalado esto último y me gustaría profundizar en por qué es así.

El modelo que OP ha especificado es un modelo de crecimiento, con pendientes e intersecciones aleatorias. Se incluye una gran intersección, pero no una pendiente de Grand Age. Una restricción desagradable que se impone ajustando una pendiente aleatoria sin la adición de su término "grandioso" es que está forzando que la pendiente aleatoria tenga 0 significa, lo cual es muy difícil de optimizar. Los modelos marginales indican que la edad no tiene un valor diferente estadísticamente significativo de 0 en el modelo. Además, agregar edad como un efecto fijo no soluciona el problema.

> lme(conc~ age, random=~age|Lot, data=IGF) 
Error in lme.formula(conc ~ age, random = ~age | Lot, data = IGF) : 
    nlminb problem, convergence error code = 1 
    message = iteration limit reached without convergence (10) 

Aquí el error es obvio. Puede ser tentador establecer el número de iteraciones. lmeControl tiene muchas estimaciones iterativas. Pero incluso eso no funciona:

> fit <- lme(conc~ 1, random=~age|Lot, data=IGF, 
control = lmeControl(maxIter = 1e8, msMaxIter = 1e8)) 

Error in lme.formula(conc ~ 1, random = ~age | Lot, 
data = IGF, control = lmeControl(maxIter = 1e+08, : 
    nlminb problem, convergence error code = 1 
    message = singular convergence (7) 

Así que no es una cuestión de precisión, el optimizador se está quedando fuera de los límites.

Debe haber diferencias clave entre los dos modelos que ha propuesto y una forma de diagnosticar el error que ha encontrado. Un enfoque simple es especificar un ajuste "prolija" para el modelo problemática:

> lme(conc~ 1, random=~age|Lot, data=IGF, control = lmeControl(msVerbose = TRUE)) 
    0:  602.96050: 2.63471 4.78706 141.598 
    1:  602.85855: 3.09182 4.81754 141.597 
    2:  602.85312: 3.12199 4.97587 141.598 
    3:  602.83803: 3.23502 4.93514 141.598 
    (truncated) 
48:  602.76219: 6.22172 4.81029 4211.89 
49:  602.76217: 6.26814 4.81000 4425.23 
50:  602.76216: 6.31630 4.80997 4638.57 
50:  602.76216: 6.31630 4.80997 4638.57 

El primer término es la REML (creo). Los términos segundo a cuarto son los parámetros de un objeto llamado lmeSt de la clase lmeStructInt, lmeStruct y modelStruct. Si usa el depurador de Rstudio para inspeccionar los atributos de este objeto (el eje del problema), verá que es el componente de efectos aleatorios que explota aquí.coef(lmeSt) después de 50 iteraciones produce reStruct.Lot1 reStruct.Lot2 reStruct.Lot3 6.316295 4.809975 4638.570586

como se ve arriba y produce

> coef(lmeSt, unconstrained = FALSE) 

    reStruct.Lot.var((Intercept)) reStruct.Lot.cov(age,(Intercept)) 
         306382.7       2567534.6 
      reStruct.Lot.var(age) 
         21531399.4 

que es el mismo que el

Browse[1]> lmeSt$reStruct$Lot 
Positive definite matrix structure of class pdLogChol representing 
      (Intercept)  age 
(Intercept) 306382.7 2567535 
age   2567534.6 21531399 

Por lo tanto, está claro que la covarianza de los efectos aleatorios es algo que está explotando aquí para este optimizador particular. Las rutinas PORT en nlminb han sido criticadas por sus errores no informativos. El texto de David Gay (Bell Labs) está aquí http://ms.mcmaster.ca/~bolker/misc/port.pdf La documentación de PORT sugiere que nuestro error 7 al usar un billon iter max "x puede tener demasiados componentes libres. Ver §5". En lugar de corregir el algoritmo, nos corresponde preguntar si hay resultados aproximados que deberían generar resultados similares. Es, por ejemplo, son fáciles de instalar un objeto lmList para llegar a la intersección aleatoria y variación de pendiente aleatoria:

> fit <- lmList(conc ~ age | Lot, data=IGF) 
> cov(coef(fit)) 
      (Intercept)   age 
(Intercept) 0.13763699 -0.018609973 
age   -0.01860997 0.003435819 

aunque idealmente deberían ser ponderados por sus respectivos pesos de precisión:

Para utilizar el nlme paquete observo que la optimización no restringida utilizando BFGS no produce tal error y da resultados similares:

> lme(conc ~ 1, data=IGF, random=~age|Lot, control = lmeControl(opt = 'optim')) 
Linear mixed-effects model fit by REML 
    Data: IGF 
    Log-restricted-likelihood: -292.9675 
    Fixed: conc ~ 1 
(Intercept) 
    5.333577 

Random effects: 
Formula: ~age | Lot 
Structure: General positive-definite, Log-Cholesky parametrization 
      StdDev  Corr 
(Intercept) 0.9976 (Intr) 
age   0.005647296 -0.698 
Residual 0.820819785  

Number of Observations: 237 
Number of Groups: 10 

una declaración sintáctica alternativa de un modelo de este tipo se puede hacer con t él mucho más fácil lme4 paquete:

library(lme4) 
lmer(conc ~ 1 + (age | Lot), data=IGF) 

que produce:

> lmer(conc ~ 1 + (age | Lot), data=IGF) 
Linear mixed model fit by REML ['lmerMod'] 
Formula: conc ~ 1 + (age | Lot) 
    Data: IGF 
REML criterion at convergence: 585.7987 
Random effects: 
Groups Name  Std.Dev. Corr 
Lot  (Intercept) 0.056254  
      age   0.006687 -1.00 
Residual    0.820609  
Number of obs: 237, groups: Lot, 10 
Fixed Effects: 
(Intercept) 
     5.331 

Un atributo de lmer y su optimizador es que las correlaciones de efectos aleatorios que están muy cerca de 1, 0 o -1 se fijan simplemente a esos valores, ya que simplifica sustancialmente la optimización (y la eficiencia estadística de la estimación).

Tomados en conjunto, esto hace no sugieren que la edad no tiene un efecto, como se dijo anteriormente, y este argumento puede ser respaldado por los resultados numéricos.