Pregunta:
¿Por qué este simple modelo mixto no logra converger?
Philippe Massicotte
2019-02-06 19:16:16 UTC
view on stackexchange narkive permalink

Tengo datos sobre los que me gustaría comparar la media de la variable valor entre el día 28 y el día 83. Debido a que la experiencia implica una pseudo-replicación ( cultura ), estaba pensando en utilizar un modelo mixto con un efecto aleatorio en cultura Tengo dos preguntas. El primero trata sobre este conjunto de datos:

 biblioteca  (lme4)
#> Cargando el paquete requerido: Matrix
  biblioteca (lmerTest)
#>
#> Adjuntando paquete: 'lmerTest'
#> El siguiente objeto está enmascarado del 'paquete: lme4':
#>
#> lmer
#> El siguiente objeto está enmascarado de 'paquete: stats':
#>
#> paso

df <- estructura (lista (
  día_verdadero = c (28, 83, 28, 83, 28, 83), valor = c (
    758453.333333333,
    575133.333333333, 684160, 656933.333333333, 816840, 734700
  ),
  cultivo = c (1L, 1L, 2L, 2L, 3L, 3L)
), fila.nombres = c (NA, -6L), clase = c ("data.frame"))
 

Si lo coloco de la siguiente manera, no tengo advertencias.

  lme4 :: lmer (value ~ factor (day_true) + (1 | culture), data = df)
#> Ajuste de modelo lineal mixto por REML ['lmerMod']
#> Fórmula: valor ~ factor (día_verdadero) + (1 | cultura)
#> Datos: df
#> Criterio REML en la convergencia: 102.7974
#> Efectos aleatorios:
#> Nombre de grupos Desv. Estándar
#> culture (Intercepción) 47535
#> residual 55990
#> Número de obs: 6, grupos: cultura, 3
Efectos fijos #>:
#> (Intercepción) factor (day_true) 83
#> 753151 -97562
 

Sin embargo, si lo encajo con lmerTest tengo una advertencia. ¿Debería preocuparme por eso?

  lmerTest :: lmer (value ~ factor (day_true) + (1 | culture), data = df)
#> Advertencia en as_lmerModLT (modelo, devfun): el modelo puede no haber convergido con 1
#> autovalor cercano a cero: 2.6e-09
#> Ajuste de modelo lineal mixto por REML ['lmerModLmerTest']
#> Fórmula: valor ~ factor (día_verdadero) + (1 | cultura)
#> Datos: df
#> Criterio REML en la convergencia: 102.7974
#> Efectos aleatorios:
#> Nombre de grupos Desv. Estándar
#> culture (Intercepción) 47535
#> residual 55990
#> Número de obs: 6, grupos: cultura, 3
Efectos fijos #>:
#> (Intercepción) factor (day_true) 83
#> 753151 -97562
 

Según el modelo resultante, no existe una diferencia significativa. Pero visualmente parece tener uno. Solo quiero asegurarme de entender el problema de convergencia.

  boxplot (valor ~ factor (day_true), data = df)
 

Mi segunda pregunta es sobre estos datos, para los cuales tengo un ajuste singular mensaje. No encuentro el motivo. ¿Es porque tengo muy pocos puntos (n = 3 por grupo)? Alternativamente, ¿existe un mejor análisis para comparar el media de medidas repetidas entre estos dos grupos?

  df <- structure (list (day_true = c (0, 28, 0, 28, 0, 28), value = c (
  34.6732447526395,
  31.5635584852635, 34.2763264775584, 32.1719125021771, 35.0747566289866,
  31.7318622838194
), cultura = c (1L, 1L, 2L, 2L, 3L, 3L)), fila.nombres = c (
  N / A,
  -6L
), class = c ("data.frame"))

  lme4 :: lmer (valor ~ factor (día_verdadero) + (1 | cultura), datos = df)
#> singular fit
#> Ajuste de modelo lineal mixto por REML ['lmerMod']
#> Fórmula: valor ~ factor (día_verdadero) + (1 | cultura)
#> Datos: df
#> Criterio REML en la convergencia: 5.3578
#> Efectos aleatorios:
#> Nombre de grupos Desv. Estándar
#> culture (Intercepción) 0.0000
#> Residual 0.3592
#> Número de obs: 6, grupos: cultura, 3
Efectos fijos #>:
#> (Intercepción) factor (day_true) 28
#> 34.675 -2.852
#> código de convergencia 0; 1 advertencias del optimizador; 0 advertencias lme4
 

Creado el 2019-02-06 por el paquete reprex (v0.2.1)

Información de la sesión

  devtools :: session_info ()
#> ─ Información de la sesión ───────────────────────────────────────────── ─────────────
#> valor de configuración
#> versión R versión 3.5.2 (2018-12-20)
#> para Linux Mint 19.1
#> system x86_64, linux-gnu
#> ui X11
#> language en_CA: en
#> collate en_CA.UTF-8
#> ctype en_CA.UTF-8
#> tz America / Toronto
#> fecha 2019-02-06
#>
#> ─ Paquetes ────────────────────────────────────────────── ─────────────────
#> paquete * versión fecha lib fuente
#> asevera que 0.2.0 2017-04-11 [1] CRAN (R 3.5.1)
#> backports 1.1.3 2018-12-14 [1] CRAN (R 3.5.1)
#> callr 3.1.1 2018-12-21 [1] CRAN (R 3.5.2)
#> cli 1.0.1 2018-09-25 [1] CRAN (R 3.5.1)
#> colorspace 1.4-0 2019-01-13 [1] CRAN (R 3.5.2)
#> crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
#> curl 3.3 2019-01-10 [1] CRAN (R 3.5.2)
#> desc 1.2.0 2018-05-01 [1] CRAN (R 3.5.1)
#> devtools 2.0.1 2018-10-26 [1] CRAN (R 3.5.1)
#> digest 0.6.18 2018-10-10 [1] CRAN (R 3.5.1)
#> dplyr 0.8.0 2019-02-01 [1] Github (tidyverse / dplyr @ e9902f1)
#> evalúa 0.12 2018-10-09 [1] CRAN (R 3.5.1)
#> fs 1.2.6 2018-08-23 [1] CRAN (R 3.5.1)
#> ggplot2 3.1.0 2018-10-25 [1] CRAN (R 3.5.1)
#> pegamento 1.3.0 2018-07-17 [1] CRAN (R 3.5.1)
#> gtable 0.2.0 2016-02-26 [1] CRAN (R 3.5.1)
#> highr 0.7 2018-06-09 [1] CRAN (R 3.5.1)
#> htmltools 0.3.6 2017-04-28 [1] CRAN (R 3.5.1)
#> httr 1.4.0 2018-12-11 [1] CRAN (R 3.5.1)
#> knitr 1.21 2018-12-10 [1] CRAN (R 3.5.1)
#> celosía 0.20-38 2018-11-04 [1] CRAN (R 3.5.1)
#> lazyeval 0.2.1 2017-10-29 [1] CRAN (R 3.5.1)
#> lme4 * 1.1-20 2019-02-04 [1] CRAN (R 3.5.2)
#> lmerTest * 3.0-1 2018-04-23 [1] CRAN (R 3.5.2)
#> magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
#> MASS 7.3-51.1 2018-11-01 [1] CRAN (R 3.5.1)
#> Matrix * 1.2-15 2018-11-01 [1] CRAN (R 3.5.1)
#> memoise 1.1.0 2017-04-21 [1] CRAN (R 3.5.1)
#> mime 0.6 2018-10-05 [1] CRAN (R 3.5.1)
#> minqa 1.2.4 2014-10-09 [1] CRAN (R 3.5.1)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
#> nlme 3.1-137 2018-04-07 [1] CRAN (R 3.5.2)
#> nloptr 1.2.1 2018-10-03 [1] CRAN (R 3.5.1)
#> numDeriv 2016.8-1 2016-08-27 [1] CRAN (R 3.5.1)
#> pilar 1.3.1.9000 2019-02-01 [1] Github (r-lib / pillar @ 3a54b8d)
#> pkgbuild 1.0.2 2018-10-16 [1] CRAN (R 3.5.1)
#> pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
#> pkgload 1.0.2 2018-10-29 [1] CRAN (R 3.5.1)
#> plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
#> prettyunits 1.0.2 2015-07-13 [1] CRAN (R 3.5.1)
#> processx 3.2.1 2018-12-05 [1] CRAN (R 3.5.1)
#> ps 1.3.0 2018-12-21 [1] CRAN (R 3.5.2)
#> purrr 0.3.0 2019-01-27 [1] CRAN (R 3.5.2)
#> R6 2.3.0 2018-10-04 [1] CRAN (R 3.5.1)
#> Rcpp 1.0.0 2018-11-07 [1] CRAN (R 3.5.1)
#> remotes 2.0.2 2018-10-30 [1] CRAN (R 3.5.1)
#> rlang 0.3.1.9000 2019-02-01 [1] Github (tidyverse / rlang @ 7243c6d)
#> rmarkdown 1.11 2018-12-08 [1] CRAN (R 3.5.1)
#> rprojroot 1.3-2 2018-01-03 [1] CRAN (R 3.5.1)
#> scale 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
#> sessioninfo 1.1.1 2018-11-05 [1] CRAN (R 3.5.1)
#> stringi 1.2.4 2018-07-20 [1] CRAN (R 3.5.1)
#> stringr 1.3.1 2018-05-10 [1] CRAN (R 3.5.1)
#> testthat 2.0.1 2018-10-13 [1] CRAN (R 3.5.1)
#> tibble 2.0.1 2019-01-12 [1] CRAN (R 3.5.2)
#> tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
#> usethis 1.4.0 2018-08-14 [1] CRAN (R 3.5.1)
#> withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
#> xfun 0.4 2018-10-23 [1] CRAN (R 3.5.1)
#> xml2 1.2.0 2018-01-24 [1] CRAN (R 3.5.1)
#> yaml 2.2.0 2018-07-25 [1] CRAN (R 3.5.1)
#>
#> [1] /home/pmassicotte/R/x86_64-pc-linux-gnu-library/3.5
#> [2] / usr / local / lib / R / site-library
#> [3] / usr / lib / R / site-library
#> [4] / usr / lib / R / library
 

No puedo replicar este problema, funciona sin errores para mí.¿Actualizaste tus paquetes?
Parece que tengo todo actualizado.Pegaré la información de mi sesión.
Creo que es de lmer in lmerTest.lmer de lme4 me da los mismos resultados, pero sin la advertencia.
Si utilizo la función de resumen en la función lmerTest lmer obtengo información adicional `Error en el cálculo de la aproximación de Satterthwaite.Se devuelve la salida del paquete lme4 se devuelve el resumen de lme4 se ha producido algún error de cálculo en lmerTest`
¿Tiene seis puntos de datos en total, o este es un ejemplo reducido para publicar?
Estos son datos reales.Entonces, tal vez haya un mejor análisis para usar.
Dos respuestas:
#1
+13
Robert Long
2019-02-06 20:30:46 UTC
view on stackexchange narkive permalink

Esto, con toda probabilidad, no es una advertencia de la que deba preocuparse. Como puede ver, las estimaciones de los parámetros son las mismas en ambos casos. La versión de lmer en lmertest aparentemente tiene una verificación de convergencia más conservadora que la versión actual de lme4 .

El problema en lmertest :: lmer se debe a que las variables están en escalas muy diferentes, lo que puede hacer que algunas de las rutinas de optimización sean inestables y, por lo tanto, generar la advertencia. Si vuelve a escalar, el problema desaparece:

  > df  $ valor <- df $  valor / 10000

> lmerTest :: lmer (valor ~ factor (día_verdadero) + (1 | cultura), datos = gl)

Ajuste de modelo lineal mixto por REML ['lmerModLmerTest']
Fórmula: valor ~ factor (día_verdadero) + (1 | cultura)
   Datos: df
Criterio REML en convergencia: 29.1147
Efectos aleatorios:
 Grupos Nombre Desv. Estándar
 cultura (intersección) 4.753
 Residual 5.599
Número de obs: 6, grupos: cultura, 3
Efectos fijos:
       (Intercepción) factor (día_verdadero) 83
            75,315 -9,756
 

Tenga en cuenta que, dado que solo tiene 3 grupos, puede ser mejor modelar la cultura como fija:

  > lm (valor ~ factor (día_verdadero) + cultura, datos = df)
Coeficientes:
   (Intercepción) factor (day_true) 83 cultura
        64,417 -9,756 5,449
 

Claramente, la estimación del efecto fijo de day_true es la misma en ambos análisis.

La razón por la que no se encuentra una estimación estadísticamente significativa es porque el tamaño de la muestra es muy pequeño. Es muy preferible ejecutar un "análisis de potencia" antes de recopilar datos y ajustar el modelo.

Gracias Robert por esta buena explicación.Usar "cultura" como fijo podría ser una solución.Sin embargo, según los resultados, parece que no hay diferencia entre el día 28 y el día 83. ¿Es correcto decir eso?Si ves la trama que hice arriba, parece que tengo una diferencia.
@PhilippeMassicotte No, no está diciendo que no haya diferencia.La diferencia se estima en -9,756 en ambos modelos, por lo que ambos dicen lo mismo.Si se refiere al valor p de la estimación o al intervalo de confianza, no preste mucha atención a esto; con una muestra tan pequeña, podría ser difícil obtener un resultado "significativo" y no tendrá suficiente poder estadístico.para detectar el tamaño del efecto (esto debería haber sido verificado antes de recolectar los datos).
¿Crees que está bien ajustar un modelo mixto para n = 6?
@amoeba no, no, por eso sugerí el modelo lineal con efectos fijos para la variable de agrupación, en mi respuesta.
#2
+8
Wayne
2019-02-06 21:25:42 UTC
view on stackexchange narkive permalink

No puedo hablar sobre el problema del cálculo, pero seis puntos de datos no son mucho, mucho menos si desea ajustar efectos aleatorios, de los cuales solo tiene 2 y solo 3 ejemplos de cada uno.

El resultado no estadísticamente significativo tiene mucho sentido aquí: tiene una pequeña cantidad de datos que apuntan a algo, pero no lo suficiente para calcular nada lo suficientemente bien como para sacar conclusiones generales.

Tus diagramas de caja te están engañando: se trata de calcular medias, cuartiles, etc. con solo 3 números para cada caja.En este caso, hay más funciones visualizadas que datos.

Gracias por tu valioso comentario.Según tengo entendido, dada la escasa cantidad de observaciones, es difícil sacar conclusiones claras.¿Ves otros lugares para abordar esto?
@PhilippeMassicotte Deshágase del diagrama de caja y simplemente trace tres puntos para los valores reales y conéctelos por líneas.En términos de prueba, es una prueba t pareada, pero tiene poco sentido ejecutar pruebas para n = 3.Solo indique que el valor disminuyó para las 3 culturas.
@PhilippeMassicotte: Lo que dice la ameba.Básicamente, no tiene suficientes datos para generalizar y hacer una declaración general sobre el mundo, pero puede hablar sobre los datos que tiene y lo que sucedió en esos casos.Si las tres líneas bajan, eso comienza a sugerir que es posible que haya un efecto allí.Lo que podría justificar la realización de un nuevo experimento con suficientes datos para determinar si ha descubierto algo sobre cómo funcionan las cosas.
Sería bueno que estos comentarios se conservaran en una respuesta (ya sea ediciones en la respuesta de @Wayne's o como una respuesta separada).Nunca puede haber suficientes explicaciones sensatas sobre qué hacer con pequeños conjuntos de datos.
Puedo agregarlos, o @amoeba también podría agregar una respuesta, o podemos dejarlos como comentarios.De cualquier manera.
¡Continúe y actualice su respuesta!


Esta pregunta y respuesta fue traducida automáticamente del idioma inglés.El contenido original está disponible en stackexchange, a quien agradecemos la licencia cc by-sa 4.0 bajo la que se distribuye.
Loading...