Tengo algunas preguntas con respecto a la especificación e interpretación de 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 el problema es la interpretación de los resultados de GLMM.

Actualmente estoy tratando de ajustar un GLMM. Estoy usando datos del censo de EE. UU . De la Base de datos del tracto longitudinal . Mis observaciones son secciones censales. Mi variable dependiente es el número de unidades de vivienda vacantes y estoy interesado en la relación entre la vacante y las variables socioeconómicas. El ejemplo aquí es simple, solo usando dos efectos fijos: porcentaje de población no blanca (raza) e ingreso familiar promedio (clase), más su interacción. Me gustaría incluir dos efectos aleatorios anidados: tratados dentro de décadas y décadas, es decir (década / tratado). Los estoy considerando al azar en un esfuerzo por controlar la autocorrelación espacial (es decir, entre tractos) y temporal (es decir, entre décadas). Sin embargo, estoy interesado en la década como un efecto fijo también, así que también lo estoy incluyendo como un factor fijo.

Como mi variable independiente es una variable de recuento de enteros no negativa, he estado tratando de ajustar poisson y GLMM binomiales negativos. Estoy usando el registro del total de unidades de vivienda como compensación. Esto significa que los coeficientes se interpretan como el efecto sobre la tasa de vacantes, no el número total de casas vacantes.

Actualmente tengo resultados para un Poisson y un GLMM binomial negativo estimado usando glmer y glmer.nb de lme4 . La interpretación de los coeficientes tiene sentido para mí según mi conocimiento de los datos y el área de estudio.

Si quieres los datos y el script , 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 de 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 estadísticas

Como todavía estoy descubriendo GLMMs, me siento inseguro sobre las especificaciones 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, estoy mejor con binomio negativo. Sin embargo, constantemente recibo 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,: se alcanzó el límite de iteración". Esto ocurre usando bastantes especificaciones diferentes (es decir, modelos mínimos y máximos para efectos fijos y aleatorios). También he intentado eliminar valores atípicos en mi dependiente (¡bruto, lo sé!), Ya que el 1% superior de los valores son muy atípicos (rango inferior del 99% de 0-1012, superior 1% de 1013-5213). No tiene ningún efecto en las iteraciones y muy poco efecto en los coeficientes tampoco. No incluyo esos detalles aquí. Los coeficientes entre Poisson y el binomio negativo también son bastante similares. ¿Es esta falta de convergencia un problema? ¿Es el modelo binomial negativo un buen ajuste? También ejecuté el modelo binomial negativo usandoAllFit y no todos los optimizadores lanzan esta advertencia (bobyqa, Nelder Mead y nlminbw no lo hicieron).

2. La variación para mi efecto fijo de la década es consistentemente muy baja o 0. Entiendo que esto podría significar que el modelo está sobreajustado. Eliminar la década de los efectos fijos aumenta la varianza de efectos aleatorios de la década a 0.2620 y no tiene mucho efecto sobre los coeficientes de efectos fijos. ¿Hay algo malo en dejarlo? Estoy bien interpretándolo como simplemente no es necesario explicar entre la varianza de observación.

3. ¿Estos resultados indican que debería probar modelos con inflación cero? DHARMa parece sugerir que la inflación cero puede no ser el problema. Si crees que debería intentarlo de todos modos, mira a continuación.

Preguntas R

4. Estaría dispuesto a probar modelos inflados a cero, pero no estoy seguro de qué paquete implica efectos aleatorios anidados para Poisson inflado a cero y GLMM binomiales negativos. Usaría glmmADMB para comparar AIC con modelos inflados a cero, pero está restringido a un solo efecto aleatorio, por lo que no funciona para este modelo. Podría probar MCMCglmm, pero no conozco las estadísticas bayesianas, por lo que tampoco es atractivo. ¿Alguna otra opción?

5. ¿Puedo mostrar coeficientes exponenciados dentro del resumen (modelo), o tengo que hacerlo fuera del resumen como lo hice aquí?

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

Por lo que obtuve al examinar sus datos, sus unidades no están agrupadas geográficamente, es decir, secciones censales dentro de los condados. Por lo tanto, el uso de tratados como un factor de agrupación no es apropiado para capturar la heterogeneidad espacial, ya que esto significa que tiene la misma cantidad de individuos que grupos (o dicho de otra manera, todos sus grupos tienen solo una observación cada uno). El uso de una estrategia de modelado multinivel nos permite estimar la varianza a nivel individual, mientras controlamos la varianza entre grupos. Como sus grupos tienen solo un individuo cada uno, su variación entre grupos es la misma que su variación 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, los puntajes de "matemáticas" de un individuo se pueden registrar 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 sección censal por decade. Por lo tanto, podría usar su TRTID10variable como un factor de agrupación para capturar "la variación entre décadas". Esto lleva a 3142 observaciones anidadas en 635 secciones, con aproximadamente 4 y 5 observaciones por sección censal.

Como se menciona en un comentario, usarlo decadecomo factor de agrupación no es muy apropiado, ya que solo tiene alrededor de 5 décadas para cada sección censal, y su efecto puede capturarse mejor introduciéndose decadecomo una covariable.

En segundo lugar, para determinar si sus datos deben modelarse utilizando un modelo de Poisson o binomial negativo (o un enfoque inflado a 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. Mirando sus datos, es bastante claro que hay mucha sobredispersión. Las variaciones son mucho más grandes que los medios.

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 el binomio negativo es estadísticamente más apropiado, un método estándar es hacer una prueba de razón de probabilidad entre un modelo de Poisson y un modelo binomial negativo, lo que sugiere que el negbin es un mejor ajuste.

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 próxima prueba podría considerar si el enfoque multinivel (modelo mixto) está garantizado utilizando un enfoque similar, lo que sugiere que la versión multinivel proporciona un mejor ajuste. (Se podría usar una prueba similar para comparar un ajuste glmer suponiendo una distribución de Poisson al objeto glmer.nb, siempre que los modelos sean los mismos).

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

Con respecto a las estimaciones de los modelos de poisson y nb, en realidad se espera que sean muy similares entre sí, siendo la principal diferencia los errores estándar, es decir, si hay una sobredispersión, el modelo de 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

Observe cómo las estimaciones de coeficientes son muy similares, la diferencia principal es solo la importancia de una de sus covariables, así como la diferencia en la variación de efectos aleatorios, lo que sugiere que la variación de nivel de unidad capturada por el parámetro de sobredispersión en el nb modelo (el thetavalor en el objeto glmer.nb) captura parte de la varianza entre tractos capturada por los efectos aleatorios.

Con respecto a los coeficientes exponenciados (y los intervalos de confianza asociados), puede usar 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, con respecto a la inflación cero. No existe una implementación multinivel (al menos que yo sepa) de un modelo de poisson o negbin inflado a cero que le permita especificar una ecuación para el componente inflado a cero de la mezcla. El glmmADMBmodelo le permite estimar un parámetro de inflación cero constante. Un enfoque alternativo sería utilizar la zeroinflfunción en el psclpaquete, aunque esto no es compatible con modelos multinivel. Por lo tanto, puede comparar el ajuste de un binomio negativo de un solo nivel, con el binomio negativo inflado de un solo nivel. 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, puede controlar esto utilizando alguna forma de regresión ponderada geográfica (aunque creo que esto utiliza datos de puntos, no áreas). Alternativamente, puede 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 si esto es completamente factible, es posible incorporar la dependencia espacial utilizando, por ejemplo, el recuento promedio de R_VACvecinos de primer orden como covariable. En cualquier caso, antes de tales enfoques, sería sensato determinar si la autocorrelación espacial está realmente presente (usando las pruebas I, LISA y enfoques similares de Global Moran).