¿Por qué los valores predichos no coinciden con el resumen del modelo binomial de GLMM?

Question

Estoy usando un GLMM para determinar si los cierres de negocios por COVID-19 afectaron la actividad de las ratas en la ciudad. La variable de respuesta es binomial (sin actividad/actividad), medida en estaciones de cebo que el consejo utiliza para desplegar cebo envenenado. El modelo tiene un predictor continuo DAY, y un predictor categórico Period, y un intercepción aleatoria para el factor aleatorio Station. Todos los predictores están dentro de la Station, pero los datos no están completos, ya que las estaciones no fueron revisadas todos los días, sino entre 5 y 10 días. El factor categórico tiene 3 niveles (Pre-cierre, Cierre, y Post-cierre).

Esto es lo que hice:

Ajusté un GLMM binomial usando glmer del paquete lme4 e hice selección de modelos usando drop1 aunque ninguna variable fue eliminada. Aquí está el modelo final:

M <- glmer(ActivityBi ~ Period + Period:DAY + (1|MBP_ID2),family=binomial(logit), data = (DATA), control = glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e5)))
drop1(M)
summary(M)

Este es el resultado:

Modelo lineal mixto generalizado ajustado por máximo verosimilitud (Aproximación de Laplace) ['glmerMod']
 Familia: binomial  ( logit )
Fórmula: ActivityBi ~ Period + Period:DAY + (1 | MBP_ID2)
   Datos: (DATA)
Control: glmerControl(optimizer = "nloptwrap", optCtrl = list(maxfun = 2e+05))

     AIC      BIC   logLik deviance df.resid 
   15028    15080    -7507    15014    12417 

Residuos escalados: 
    Min      1Q  Median      3Q     Máx 
-3.2887 -0.7699 -0.2530  0.7782  4.7020 

Efectos aleatorios:
 Grupos  Nombre        Varianza Desv.Est.
 MBP_ID2 (Intercepción) 1.007    1.004   
Número de observaciones: 12424, grupos:  MBP_ID2, 659

Efectos fijos:
                           Estimación Error estándar valor z valor p    
(Intercepto)               0.1047790  0.0666821   1.571   0.1161    
Period2.Cierre            8.5385820  1.2521507   6.819 9.16e-12 ***
Period3.Postcierre       -0.0029200  0.7599282  -0.004   0.9969    
Period1.Precierre:DAY     0.0007355  0.0004266   1.724   0.0847 .  
Period2.Cierre:DAY       -0.0408955  0.0053535  -7.639 2.19e-14 ***
Period3.Postcierre:DAY   -0.0066198  0.0025231  -2.624   0.0087 ** 
---
Códigos de significancia:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlación de Efectos Fijos:
            (Interc) Prd2.L Prd3.P P1.P:D P2.L:D
Prd2.Crr -0.033                            
Prd3.Pstc -0.059  0.005                     
Prd1.Pr:Dy -0.685  0.043  0.069              
Prd2.Crr:Dy -0.002 -0.997 -0.001 -0.002       
Prd3.Pstc:Dy -0.002 -0.003 -0.993 -0.002  0.003
código de convergencia: 0
El modelo no convergió con max|grad| = 0.0633013 (tol = 0.002, componente 1)
El modelo es casi no identificable: eigenvalor muy grande
 - ¿Reescalar variables?
El modelo es casi no identificable: grande relación de eigenvalores
 - ¿Reescalar variables?

Luego uso Anova de car para calcular los valores p para el factor fijo y las interacciones, y emmeans para calcular las estimaciones para cada nivel de periodo

> Anova(M, test="Chisq", type=3)
Tabla de análisis de desviación (pruebas de Chi-cuadrado tipo III Wald)

Respuesta: ActivityBi
              Chisq Df Pr(>Chisq)    
(Intercepto)  2.4691  1     0.1161    
Periodo      46.5020  2  7.984e-11 ***
Periodo:DAY  68.0141  3  1.136e-14 ***
---
Códigos de significancia:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> marginal = emmeans(M, ~ Periodo|DAY, type="response") 
NOTA: Los resultados pueden ser engañosos debido a la participación en interacciones
> pairs(marginal2,adjust="tukey")
 contraste                       razón de momios    EE  df z.ratio valor p
 1.Precierre / 2.Cierre         0.0844  0.04 Inf -5.210  <.0001 
 1.Precierre / 3.Postcierre     2.9291  1.15 Inf  2.728  0.0175 
 2.Cierre / 3.Postcierre       34.7132 21.36 Inf  5.766  <.0001 

Ajuste de valor p: método tukey para comparar una familia de 3 estimaciones 
Las pruebas se realizan en la escala de razón de momios logarítmica 
> cld(marginal2,alpha=0.05, Letras=letters, ajuste="tukey")
 Periodo          prob     EE  df intervalo.l.inf asint.UCL grupo
 3.Postlockdown 0.297 0.0825 Inf     0.141     0.520  a    
 1.Prelockdown  0.553 0.0127 Inf     0.522     0.583   b   
 2.Cierre    0.936 0.0285 Inf     0.825     0.979    c 

Como se puede ver en el resumen del modelo y los resultados de emmeans, la actividad de las ratas durante el cierre fue mayor que antes y después del cierre.

Finalmente, debido al efecto de DAY quiero calcular las estimaciones para graficar. Para eso, utilicé predict() pero sigo recibiendo este error:

> EstM<-predict(M, newdata=data, nsim=100, type="response", interval="confidence", re.form=NA, se.fit=TRUE, na.rm = TRUE)
Mensaje de advertencia:
En lme4:::predict.merMod(x, ...) : argumentos no utilizados ignorados

Luego intenté usar ```predictInterval``` en su lugar

EstM2 <- predictInterval(M, newdata = data, which="fixed", n.sims=1000, type = "linear.prediction")

Pero los valores predichos no coinciden con el resumen del modelo, con la actividad de las ratas antes del cierre siendo mayor que durante el cierre.

Incluso si uso "probabilidad"

EstM2 <- predictInterval(M, newdata = data, which="fixed", n.sims=1000, type = "probability")

En ambos casos, los valores antes del cierre parecen ser más altos que durante el cierre. ¿Alguna idea de por qué está sucediendo esto?

Answer 1

Creo que el problema aquí es que hay una dependencia entre Período y DÍA. Dado que DÍA es una covariable, emmeans() lo reduce a su valor promedio cuando determina la cuadrícula de referencia en la que se basan los EMMs.

Aquí hay un ejemplo más simple (y reproducible) que ilustra esto:

set.seed(1.2021)
fake = data.frame(Día = 1:50, 
                  Período = factor(c(rep("A",20), rep("B",10), rep("C",20))),
                  y = rnorm(50) + 100 - 1:50)
plot(y ~ Día, col = Período, data = fake)

Aquí podemos ver que hay una tendencia lineal con Día, y los períodos comprenden los primeros 20 días, los siguientes 10 días, y los últimos 20 días, respectivamente.

Ahora vamos a ajustar un modelo y obtener los EMMs:

fake.lm = lm(y ~ Período + Día:Período, data = fake)

library(emmeans)

emmeans(fake.lm, "Período")
#> NOTA: Los resultados pueden ser engañosos debido a la participación en interacciones
#>  Período emmean    SE df lower.CL upper.CL
#>  A        75.0 0.531 44     73.9     76.1
#>  B        74.4 0.270 44     73.8     74.9
#>  C        74.5 0.531 44     73.4     75.6
#> 
#> Nivel de confianza utilizado: 0.95

Todos los EMMs son aproximadamente 75. Esto es porque estamos prediciendo todo en el mismo valor promedio de Día:

ref_grid(fake.lm)@grid
#>   Período  Día .wgt.
#> 1      A 25.5    20
#> 2      B 25.5    10
#> 3      C 25.5    20

Pero ahora vamos a tener en cuenta la dependencia entre Período y Día requiriendo que calculemos el promedio de Día para cada Período:

ref_grid(fake.lm, cov.reduce = Día ~ Período)@grid
#>   Período  Día .wgt.
#> 1      A 10.5    20
#> 2      B 25.5    10
#> 3      C 40.5    20

Con este cambio, los EMMs reflejan los promedios que vemos en el gráfico:

emmeans(fake.lm, "Período", cov.reduce = Día ~ Período)
#> NOTA: Los resultados pueden ser engañosos debido a la participación en interacciones
#>  Período emmean     SE df lower.CL upper.CL
#>  A       89.69 0.1906 44    89.31    90.07
#>  B       74.37 0.2696 44    73.82    74.91
#>  C       59.63 0.1906 44    59.24    60.01
#> 
#> Nivel de confianza utilizado: 0.95

Creado el 2021-01-25 por el paquete reprex (v0.3.0)

En este ejemplo, los efectos de Período e interacción son muy cercanos a cero porque los valores generados de y son simplemente una función lineal de Día. Eso hace que todos los EMMs en el primer conjunto sean aproximadamente iguales. En tu ejemplo, realmente hay efectos de Período e interacción, lo que hace que los resultados sean más irregulares.