Advertencia "El modelo no converge" en lmer()

Question

Con el siguiente conjunto de datos, quería ver si la respuesta (efecto) cambia con respecto a los sitios, la estación, la duración y sus interacciones. Algunos foros en línea sobre estadística me sugirieron que siguiera con los modelos lineales de efectos mixtos, pero el problema es que, dado que las réplicas son aleatorias dentro de cada estación, tengo pocas posibilidades de recoger la muestra exactamente en el mismo lugar en estaciones sucesivas (por ejemplo, la réplica 1 de s1 del postmonzón puede no ser la misma que la del monzón). No es como los ensayos clínicos (con diseño dentro de un sujeto), en los que se mide al mismo sujeto repetidamente a lo largo de las estaciones. Sin embargo, considerando los sitios y la estación como un factor aleatorio, ejecuté los siguientes comandos y recibí un mensaje de advertencia:

Warning messages:
1: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), optpar,ctrl=controlpar,ctrl=controlcheckConv, 
: Model failed to converge: degenerate Hessian with 1 negative eigenvalues

¿Puede alguien ayudarme a resolver el problema? Los códigos son los siguientes:

library(lme4)
read.table(textConnection("duration season  sites   effect
                          4d    mon s1  7305.91
                          4d    mon s2  856.297
                          4d    mon s3  649.93
                          4d    mon s1  10121.62
                          4d    mon s2  5137.85
                          4d    mon s3  3059.89
                          4d    mon s1  5384.3
                          4d    mon s2  5014.66
                          4d    mon s3  3378.15
                          4d    post    s1  6475.53
                          4d    post    s2  2923.15
                          4d    post    s3  554.05
                          4d    post    s1  7590.8
                          4d    post    s2  3888.01
                          4d    post    s3  600.07
                          4d    post    s1  6717.63
                          4d    post    s2  1542.93
                          4d    post    s3  1001.4
                          4d    pre s1  9290.84
                          4d    pre s2  2199.05
                          4d    pre s3  1149.99
                          4d    pre s1  5864.29
                          4d    pre s2  4847.92
                          4d    pre s3  4172.71
                          4d    pre s1  8419.88
                          4d    pre s2  685.18
                          4d    pre s3  4133.15
                          7d    mon s1  11129.86
                          7d    mon s2  1492.36
                          7d    mon s3  1375
                          7d    mon s1  10927.16
                          7d    mon s2  8131.14
                          7d    mon s3  9610.08
                          7d    mon s1  13732.55
                          7d    mon s2  13314.01
                          7d    mon s3  4075.65
                          7d    post    s1  11770.79
                          7d    post    s2  4254.88
                          7d    post    s3  753.2
                          7d    post    s1  11324.95
                          7d    post    s2  5133.76
                          7d    post    s3  2156.2
                          7d    post    s1  12103.76
                          7d    post    s2  3143.72
                          7d    post    s3  2603.23
                          7d    pre s1  13928.88
                          7d    pre s2  3208.28
                          7d    pre s3  8015.04
                          7d    pre s1  11851.47
                          7d    pre s2  6815.31
                          7d    pre s3  8478.77
                          7d    pre s1  13600.48
                          7d    pre s2  1219.46
                          7d    pre s3  6987.5
                          "),header=T)->dat1

m1 = lmer(effect ~ duration + (1+duration|sites) +(1+duration|season),
          data=dat1, REML=FALSE)

Answer 1

"Solucionar" el problema que experimentas en el sentido de no recibir avisos sobre la convergencia fallida es bastante sencillo: no utilizas el BOBYQA pero en su lugar se opta por utilizar el Nelder-Mead rutina de optimización utilizada por defecto en anteriores 1.0.x versiones anteriores. O bien instala el paquete optimx por lo que puede directamente un L-BFGS-B de rutina o nlminb (igual que lme4 versiones anteriores a la ver. 1 ). Por ejemplo:

m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(optimizer ="Nelder_Mead")
library(optimx)
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='L-BFGS-B')))
m1 = lmer(effect~duration+(1+duration|sites)+(1+duration|season), 
          data = dat1, REML = FALSE, 
          control = lmerControl(
                           optimizer ='optimx', optCtrl=list(method='nlminb')))

todo funciona bien (sin advertencias). Las preguntas interesantes son:

1. por qué tienes estas advertencias para empezar y
2. por qué cuando usaste REML = TRUE no tienes ninguna advertencia.

Sucintamente, 1. usted recibió esas advertencias porque definió duration tanto como efecto fijo como pendiente aleatoria para el factor sites así como season . El modelo efectivamente se quedó sin los grados de libertad para estimar las correlaciones entre las pendientes y los interceptos que usted definió. Si usted utilizó un modelo marginalmente más simple como:

m1 = lmer(effect~duration+ (1+duration|sites) + (0+duration|season) + (1|season),
          data=dat1, REML = FALSE)

no experimentará problemas de convergencia. Este modelo estimaría efectivamente interceptos aleatorios no correlacionados y pendientes aleatorias para cada season .

Además, 2. al definir REML = FALSE has utilizado la Máxima Verosimilitud estimada en lugar de la Máxima Verosimilitud Restringida. Las estimaciones REML intentan "eliminar" la influencia de los efectos fijos $X$ antes de pasar a encontrar la estructura óptima de la varianza de los efectos aleatorios (véase el hilo " ¿Qué es la "máxima probabilidad restringida" y cuándo debe utilizarse? " para obtener información más detallada al respecto). Desde el punto de vista computacional, este procedimiento se realiza esencialmente multiplicando las dos partes de la ecuación original del modelo LME $y = X\beta + Z\gamma + \epsilon$ por una matriz $K$ tal que $KX = 0$ es decir, cambias tanto el original $y$ a $Ky$ así como el $Z$ a $KZ$ . Sospecho fuertemente que esto afectó el número de condición de la matriz de diseño $Z$ y, por tanto, le ayudarán a salir de la difícil situación numérica en la que se encontraba en primer lugar.

Una última nota es que no estoy seguro de si tiene sentido utilizar season como efecto aleatorio para empezar. Al fin y al cabo, sólo hay un número determinado de temporadas, así que es mejor tratarlas como efectos fijos.

Answer 2

La cuestión es más estadística que técnica. En realidad, utilicé un modelo de efectos aleatorios en lugar de uno de efectos fijos. Creo que ninguno de los factores debería tratarse como factor aleatorio, ya que necesitamos al menos 5 o 6 niveles o réplicas para tratar un factor como efecto aleatorio (véase aquí ¿Cuál es el número mínimo recomendado de grupos para un factor de efectos aleatorios? ).

El conjunto de datos anterior sólo contiene muestras por triplicado/sitio/temporada, lo que es insuficiente para un modelo de efectos aleatorios. En el conjunto de datos, la duración, de 4 días y de 7 días, pertenece a dos experimentos paralelos separados realizados en el mismo tiempo. Por lo tanto, dividir el conjunto de datos por duración (4 y 7 días) y realizar un anova de 2 vías para cada duración con la estación y los sitios como factores sería suficiente para modelar el efecto (variable de respuesta) aquí. El modelo debería ser el siguiente

lm(day_4_effect~sites*season, data=dat1)

lm(day_7_effect~sites*season, data=dat1)

Gracias a Bodo Winter ( http://www.bodowinter.com/tutorial/bw_LME_tutorial2.pdf ) y @usεr11852 que me ayudó a resolver el problema.

Answer 3

He tenido éxito con el cambio de optimizador sugerido aquí en muestras más grandes, pero esta solución no siempre funciona en muestras relativamente pequeñas en mi experiencia. Matuschek et al., 2017 (Journal of Memory and Language) señalan que la convergencia fallida del modelo también puede ser resultado de una especificación excesiva de la estructura de efectos aleatorios cuando los datos no pueden soportar modelos máximos o casi máximos. En este caso, tendrá que simplificar su estructura de efectos aleatorios, para lo cual Matuschek et al., 2017 ofrecen algunas sugerencias. Al menos asumiendo que he interpretado su trabajo correctamente, los modelos de efectos mixtos son confusos.