I have some questions regarding specification and interpretation of GLMMs. 3 questions are definitely statistical and 2 are more specifically about R. I am posting here because ultimately I think the issue is interpretation of GLMM results.
I am currently trying to fit a GLMM. I'm using US census data from the Longitudinal Tract Database. My observations are census tracts. My dependent variable is the number of vacant housing units and I'm interested in the relationship between vacancy and socio-economic variables. The example here is simple, just using two fixed effects: percent non-white population (race) and median household income (class), plus their interaction. I would like to include two nested random effects: tracts within decades and decades, i.e. (decade/tract). I am considering these random in an effort to control for spatial (i.e. between tracts) and temporal (i.e. between decades) autocorrelation. However, I am interested in decade as a fixed effect also, so I am including it as a fixed factor also.
Since my independent variable is a non-negative integer count variable, I've been trying to fit poisson and negative binomial GLMMs. I am using the log of total housing units as an offset. This means coefficients are interpreted as the effect on vacancy rate, not total number of vacant houses.
I currently have results for a Poisson and a negative binomial GLMM estimated using glmer and glmer.nb from lme4. The interpretation of coefficients makes sense to me based on my knowledge of the data and study area.
Here are my results:
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)
Negative binomial model
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
Poisson DHARMa tests
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
Negative binomial DHARMa tests
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
Since I am still figuring out GLMMs I am feeling insecure about specification and interpretation. I have some questions:
It appears my data do not support using a Poisson model and therefore I am better off with negative binomial. However, I consistently get warnings that my negative binomial models reach their iteration limit, even when I increase the maximum limit. "In theta.ml(Y, mu, weights = object@resp$weights, limit = limit, : iteration limit reached." This happens using quite a few different specifications (i.e. minimial and maximal models for both fixed and random effects). I have also tried removing outliers in my dependent (gross, I know!), since the top 1% of values are very much outliers (bottom 99% range from 0-1012, top 1% from 1013-5213). That didn't have any effect on iterations and very little effect on coefficients either. I don't include those details here. Coefficients between Poisson and negative binomial are also pretty similar. Is this lack of convergence a problem? Is the negative binomial model a good fit? I have also run the negative binomial model using AllFit and not all optimizers throw this warning (bobyqa, Nelder Mead, and nlminbw did not).
The variance for my decade fixed effect is consistently very low or 0. I understand this could mean the model is overfit. Taking decade out of the fixed effects does increase the decade random effect variance to 0.2620 and doesn't have much of an effect on fixed effect coefficients. Is there anything wrong with leaving it in? I am fine interpreting it as simply not being needed to explain between observation variance.
Do these results indicate I should try zero-inflated models? DHARMa seems to suggest zero-inflation may not be the issue. If you think I should try anyway, see below.
I would be willing to try zero-inflated models, but I am not sure which package impliments nested random effects for zero-inflated Poisson and negative binomial GLMMs. I would use glmmADMB to compare AIC with zero-inflated models, but it is restricted to a single random effect so doesn't work for this model. I could try MCMCglmm, but I do not know Bayesian statistics so that is also not attractive. Any other options?
Can I display exponentiated coefficients within summary(model), or do I have to do it outside summary as I've done here?