Ayuda para interpretar datos de recuento GLMM usando lme4 glmer y glmer.nb - Binomio negativo versus Poisson

Question

Tengo algunas preguntas sobre la especificación e interpretación de los GLMM. 3 preguntas son definitivamente estadísticas y 2 son más específicamente sobre R. Estoy publicando aquí porque en última instancia, creo que la cuestión es la interpretación de los resultados GLMM.

Actualmente estoy intentando montar una GLMM. Estoy utilizando los datos del censo de EE.UU. del Base de datos de trazos longitudinales . Mis observaciones son secciones censales. Mi variable dependiente es el número de viviendas vacías y estoy interesado en la relación entre las vacantes y las variables socioeconómicas. El ejemplo es sencillo, sólo se utilizan dos efectos fijos: el porcentaje de población no blanca (raza) y la renta media de los hogares (clase), más su interacción. Me gustaría incluir dos efectos aleatorios anidados: tramos dentro de las décadas y décadas, es decir, (década/tramo). Estoy considerando estos efectos aleatorios en un esfuerzo por controlar la autocorrelación espacial (es decir, entre tramos) y temporal (es decir, entre décadas). Sin embargo, también me interesa la década como efecto fijo, por lo que también la incluyo como factor fijo.

Dado que mi variable independiente es una variable de recuento de enteros no negativos, he intentado ajustar GLMMs de poisson y binomial negativa. Estoy utilizando el logaritmo de las unidades de vivienda totales como compensación. Esto significa que los coeficientes se interpretan como el efecto sobre la tasa de vacantes, no sobre el número total de viviendas vacías.

Actualmente tengo los resultados de un GLMM de Poisson y uno binomial negativo estimados con glmer y glmer.nb de lme4 . La interpretación de los coeficientes tiene sentido para mí en base a mi conocimiento de los datos y del área de estudio.

Si quiere el datos y guión están en mi Github . El guión incluye más de las investigaciones descriptivas que hice antes de construir los modelos.

Aquí están mis resultados:

Modelo de Poisson

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) +      (1 | decade/TRTID10)
   Data: scaled.mydata

     AIC      BIC   logLik deviance df.resid 
 34520.1  34580.6 -17250.1  34500.1     3132 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.24211 -0.10799 -0.00722  0.06898  0.68129 

Random effects:
 Groups         Name        Variance Std.Dev.
 TRTID10:decade (Intercept) 0.4635   0.6808  
 decade         (Intercept) 0.0000   0.0000  
Number of obs: 3142, groups:  TRTID10:decade, 3142; decade, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -3.612242   0.028904 -124.98  < 2e-16 ***
decade1980       0.302868   0.040351    7.51  6.1e-14 ***
decade1990       1.088176   0.039931   27.25  < 2e-16 ***
decade2000       1.036382   0.039846   26.01  < 2e-16 ***
decade2010       1.345184   0.039485   34.07  < 2e-16 ***
P_NONWHT         0.175207   0.012982   13.50  < 2e-16 ***
a_hinc          -0.235266   0.013291  -17.70  < 2e-16 ***
P_NONWHT:a_hinc  0.093417   0.009876    9.46  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980  -0.693                                            
decade1990  -0.727  0.501                                     
decade2000  -0.728  0.502  0.530                              
decade2010  -0.714  0.511  0.517  0.518                       
P_NONWHT     0.016  0.007 -0.016 -0.015  0.006                
a_hinc      -0.023 -0.011  0.023  0.022 -0.009  0.221         
P_NONWHT:_h  0.155  0.035 -0.134 -0.129  0.003  0.155   -0.233
convergence code: 0
Model failed to converge with max|grad| = 0.00181132 (tol = 0.001, component 1)

Modelo binomial negativo

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: Negative Binomial(25181.5)  ( log )
Formula: R_VAC ~ decade + P_NONWHT + a_hinc + P_NONWHT * a_hinc + offset(HU_ln) +      (1 | decade/TRTID10)
   Data: scaled.mydata

     AIC      BIC   logLik deviance df.resid 
 34522.1  34588.7 -17250.1  34500.1     3131 

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-2.24213 -0.10816 -0.00724  0.06928  0.68145 

Random effects:
 Groups         Name        Variance  Std.Dev. 
 TRTID10:decade (Intercept) 4.635e-01 6.808e-01
 decade         (Intercept) 1.532e-11 3.914e-06
Number of obs: 3142, groups:  TRTID10:decade, 3142; decade, 5

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -3.612279   0.028946 -124.79  < 2e-16 ***
decade1980       0.302897   0.040392    7.50 6.43e-14 ***
decade1990       1.088211   0.039963   27.23  < 2e-16 ***
decade2000       1.036437   0.039884   25.99  < 2e-16 ***
decade2010       1.345227   0.039518   34.04  < 2e-16 ***
P_NONWHT         0.175216   0.012985   13.49  < 2e-16 ***
a_hinc          -0.235274   0.013298  -17.69  < 2e-16 ***
P_NONWHT:a_hinc  0.093417   0.009879    9.46  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) dc1980 dc1990 dc2000 dc2010 P_NONWHT a_hinc
decade1980  -0.693                                            
decade1990  -0.728  0.501                                     
decade2000  -0.728  0.502  0.530                              
decade2010  -0.715  0.512  0.517  0.518                       
P_NONWHT     0.016  0.007 -0.016 -0.015  0.006                
a_hinc      -0.023 -0.011  0.023  0.022 -0.009  0.221         
P_NONWHT:_h  0.154  0.035 -0.134 -0.129  0.003  0.155   -0.233

Pruebas de Poisson DHARMa

    One-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.044451, p-value = 8.104e-06
alternative hypothesis: two-sided

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model

data:  simulationOutput
ratioObsExp = 1.3666, p-value = 0.159
alternative hypothesis: more

Pruebas binomiales negativas DHARMa

    One-sample Kolmogorov-Smirnov test

data:  simulationOutput$scaledResiduals
D = 0.04263, p-value = 2.195e-05
alternative hypothesis: two-sided

    DHARMa zero-inflation test via comparison to expected zeros with simulation under H0 = fitted model

data:  simulationOutput2
ratioObsExp = 1.376, p-value = 0.174
alternative hypothesis: more

Parcelas DHARMa

Poisson

Binomio negativo

Preguntas de estadística

Dado que todavía estoy descubriendo los GLMM, me siento inseguro en cuanto a la especificación y la interpretación. Tengo algunas preguntas:

1. Parece que mis datos no admiten el uso de un modelo de Poisson y, por lo tanto, es mejor utilizar el modelo binomial negativo. Sin embargo, recibo constantemente advertencias de que mis modelos binomiales negativos alcanzan su límite de iteración, incluso cuando aumento el límite máximo. "En theta.ml(Y, mu, pesos = objeto@resp$pesos, límite = límite, : límite de iteración alcanzado". Esto ocurre utilizando bastantes especificaciones diferentes (es decir, modelos mínimos y máximos tanto para efectos fijos como aleatorios). También he probado a eliminar los valores atípicos en mi dependiente (¡bruto, lo sé!), ya que el 1% superior de los valores son muy atípicos (el 99% inferior va de 0 a 1012, el 1% superior de 1013 a 5213). Eso no tuvo ningún efecto en las iteraciones y tampoco tuvo mucho efecto en los coeficientes. No incluyo esos detalles aquí. Los coeficientes entre Poisson y la binomial negativa también son bastante similares. ¿Es esta falta de convergencia un problema? ¿Es el modelo binomial negativo un buen ajuste? También he ejecutado el modelo binomial negativo utilizando AllFit y no todos los optimizadores lanzan esta advertencia (bobyqa, Nelder Mead y nlminbw no lo hicieron).

2. La varianza de mi efecto fijo de década es sistemáticamente muy baja o 0. Entiendo que esto podría significar que el modelo está sobreajustado. Sacar la década de los efectos fijos aumenta la varianza del efecto aleatorio de la década a 0,2620 y no tiene mucho efecto en los coeficientes del efecto fijo. ¿Hay algo malo en dejarlo? Me parece bien interpretar que simplemente no es necesario para explicar la varianza entre observaciones.

3. ¿Indican estos resultados que debería probar con modelos de inflación cero? DHARMa parece sugerir que la inflación cero puede no ser el problema. Si cree que debería intentarlo de todos modos, vea más abajo.

Preguntas de la R

4. Estaría dispuesto a probar modelos con inflación cero, pero no estoy seguro de qué paquete implanta efectos aleatorios anidados para GLMMs de Poisson y binomial negativa con inflación cero. Utilizaría glmmADMB para comparar el AIC con modelos de inflación cero, pero está restringido a un solo efecto aleatorio, por lo que no funciona para este modelo. Podría probar con MCMCglmm, pero no sé de estadística bayesiana, así que tampoco es atractivo. ¿Alguna otra opción?

5. ¿Puedo mostrar los coeficientes exponenciados dentro de summary(model), o tengo que hacerlo fuera de summary como he hecho aquí?

Answer 1

Creo que hay algunos problemas importantes que hay que abordar con su estimación.

Por lo que he deducido al examinar sus datos, sus unidades no están agrupadas geográficamente, es decir, las secciones censales dentro de los condados. Por lo tanto, el uso de tramos como factor de agrupación no es apropiado para capturar la heterogeneidad espacial, ya que esto significa que usted tiene el mismo número de individuos como grupos (o dicho de otro modo, todos sus grupos tienen una sola observación cada uno). El uso de una estrategia de modelización multinivel nos permite estimar la varianza a nivel individual, controlando al mismo tiempo la varianza entre grupos. Dado que sus grupos sólo tienen un individuo cada uno, su varianza entre grupos es la misma que su varianza a nivel individual, lo que anula el propósito del enfoque multinivel.

Por otro lado, el factor de agrupación puede representar mediciones repetidas a lo largo del tiempo. Por ejemplo, en el caso de un estudio longitudinal, las puntuaciones de "matemáticas" de un individuo pueden recortarse anualmente, por lo que tendríamos un valor anual para cada estudiante durante n años (en este caso, el factor de agrupación es el estudiante, ya que tenemos n número de observaciones "anidadas" dentro de los estudiantes). En su caso, tiene medidas repetidas de cada tramo censal por decade . Así, podría utilizar su TRTID10 como factor de agrupación para captar la "varianza entre décadas". Esto lleva a 3142 observaciones anidadas en 635 tramos, con aproximadamente 4 y 5 observaciones por tramo censal.

Como se ha mencionado en un comentario, el uso de decade como factor de agrupación no es muy adecuado, ya que sólo se dispone de unas 5 décadas para cada sección censal, y su efecto puede captarse mejor introduciendo decade como covariable.

En segundo lugar, para determinar si sus datos deben modelarse mediante un modelo poisson o binomial negativo (o un enfoque de inflación cero). Considere la cantidad de sobredispersión en sus datos. La característica fundamental de una distribución de Poisson es la equidispersión, lo que significa que la media es igual a la varianza de la distribución. Observando tus datos, está bastante claro que hay mucha sobredispersión. Las varianzas son mucho mayores que las medias.

library(dplyr)    
 dispersionstats <- scaled.mydata %>%
 + group_by(decade) %>%
 + summarise(
 + means = mean(R_VAC),
 + variances = var(R_VAC),
 + ratio = variances/means)

##   dispersionstats
##   # A tibble: 5 x 5
##   decade     means variances     ratio 
##    <int>     <dbl>     <dbl>     <dbl> 
## 1   1970  45.43513   4110.89  90.47822 
## 2   1980 103.52365  17323.34 167.33707 
## 3   1990 177.68038  62129.65 349.67087 
## 4   2000 190.23150  91059.60 478.67784 
## 5   2010 247.68246 126265.60 509.78821 

No obstante, para determinar si la binomial negativa es más apropiada desde el punto de vista estadístico, un método estándar es hacer una prueba de razón de verosimilitud entre un modelo de Poisson y una binomial negativa, que sugiere que la negbin se ajusta mejor.

library(MASS)
library(lmtest)

modelformula <- formula(R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln))

poismodel <- glm(modelformula, data = scaled.mydata, family = "poisson")   
nbmodel <- glm.nb(modelformula, data = scaled.mydata)

lrtest(poismodel, nbmodel)

## Likelihood ratio test

##  Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)  
## Model 2: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -154269
## 2   9  -17452  1 273634  < 2.2e-16 ***
##  ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Después de establecer esto, una siguiente prueba podría considerar si el enfoque multinivel (modelo mixto) se justifica utilizando un enfoque similar, lo que sugiere que la versión multinivel proporciona un mejor ajuste. (Una prueba similar podría utilizarse para comparar un ajuste de glmer asumiendo una distribución poisson con el objeto glmer.nb, siempre que los modelos sean por lo demás iguales).

library(lme4)

glmmformula <- update(modelformula, . ~ . + (1|TRTID10))

nbglmm <- glmer.nb(glmmformula, data = scaled.mydata)

lrtest(nbmodel, nbglmm)

## Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
## Model 2: R_VAC ~ factor(decade) + P_NONWHT + a_hinc + (1 | TRTID10) +
##     P_NONWHT:a_hinc + offset(HU_ln)
##   #Df LogLik Df Chisq Pr(>Chisq)
## 1   9 -17452
## 2  10 -17332  1 239.3  < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

En cuanto a las estimaciones de los modelos poisson y nb, en realidad se espera que sean muy similares entre sí, siendo la principal distinción los errores estándar, es decir, si hay sobredispersión, el modelo poisson tiende a proporcionar errores estándar sesgados. Tomando sus datos como ejemplo:

poissonglmm <- glmer(glmmformula, data = scaled.mydata)
summary(poissonglmm)

## Random effects:
##  Groups  Name        Variance Std.Dev.
## TRTID10 (Intercept) 0.2001   0.4473
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.876013   0.020602 -139.60   <2e-16 ***
## factor(decade)1980  0.092597   0.007602   12.18   <2e-16 ***
## factor(decade)1990  0.903543   0.007045  128.26   <2e-16 ***
## factor(decade)2000  0.854821   0.006913  123.65   <2e-16 ***
## factor(decade)2010  0.986126   0.006723  146.67   <2e-16 ***
## P_NONWHT           -0.125500   0.014007   -8.96   <2e-16 ***
## a_hinc             -0.107335   0.001480  -72.52   <2e-16 ***
## P_NONWHT:a_hinc     0.160937   0.003117   51.64   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(nbglmm)
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  TRTID10 (Intercept) 0.09073  0.3012
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.797861   0.056214  -49.77  < 2e-16 ***
## factor(decade)1980  0.118588   0.039589    3.00  0.00274 **
## factor(decade)1990  0.903440   0.038255   23.62  < 2e-16 ***
## factor(decade)2000  0.843949   0.038172   22.11  < 2e-16 ***
## factor(decade)2010  1.068025   0.037376   28.58  < 2e-16 ***
## P_NONWHT            0.020012   0.089224    0.22  0.82253
## a_hinc             -0.129094   0.008109  -15.92  < 2e-16 ***
## P_NONWHT:a_hinc     0.149223   0.018967    7.87 3.61e-15 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Obsérvese cómo las estimaciones de los coeficientes son todas muy similares, siendo la principal diferencia únicamente la importancia de una de sus covariables, así como la diferencia en la varianza de los efectos aleatorios, lo que sugiere que la varianza a nivel de unidad capturada por el parámetro de sobredispersión en el modelo nb (el theta en el objeto glmer.nb) captura parte de la varianza entre tramos capturada por los efectos aleatorios.

En cuanto a los coeficientes exponenciales (y los intervalos de confianza asociados), puede utilizar lo siguiente:

fixed <- fixef(nbglmm)
confnitfixed <- confint(nbglmm, parm = "beta_", method = "Wald") # Beware: The Wald method is less accurate but much, much faster.

# The exponentiated coefficients are also known as Incidence Rate Ratios (IRR)
IRR <- exp(cbind(fixed, confintfixed)
IRR
##                         fixed      2.5 %     97.5 %
## (Intercept)        0.06094028 0.05458271 0.06803835
## factor(decade)1980 1.12590641 1.04184825 1.21674652
## factor(decade)1990 2.46807856 2.28979339 2.66024515
## factor(decade)2000 2.32553168 2.15789585 2.50619029
## factor(decade)2010 2.90962703 2.70410073 3.13077444
## P_NONWHT           1.02021383 0.85653208 1.21517487
## a_hinc             0.87889172 0.86503341 0.89297205
## P_NONWHT:a_hinc    1.16093170 1.11856742 1.20490048

Reflexiones finales, sobre la inflación cero. No hay ninguna implementación multinivel (al menos que yo conozca) de un modelo poisson o negbin de inflación cero que permita especificar una ecuación para el componente de inflación cero de la mezcla. el glmmADMB permite estimar un parámetro de inflación constante y nula. Un enfoque alternativo sería utilizar el zeroinfl en la función pscl aunque éste no admite modelos multinivel. Por lo tanto, podría comparar el ajuste de una binomial negativa de un solo nivel con la binomial negativa de un solo nivel con inflación cero. Lo más probable es que si la inflación cero no es significativa para los modelos de un solo nivel, es probable que no sea significativa para la especificación multinivel.

Apéndice

Si le preocupa la autocorrelación espacial, podría controlarla utilizando alguna forma de regresión geográfica ponderada (aunque creo que ésta utiliza datos puntuales, no áreas). Alternativamente, podría agrupar sus secciones censales por un factor de agrupación adicional (estados, condados) e incluir esto como un efecto aleatorio. Por último, y no estoy seguro de que esto sea totalmente factible, puede ser posible incorporar la dependencia espacial utilizando, por ejemplo, el recuento medio de R_VAC en los vecinos de primer orden como covariable. En cualquier caso, antes de estos enfoques, sería sensato determinar si la autocorrelación espacial está realmente presente (utilizando la I de Moran global, las pruebas LISA y enfoques similares).