Eu tenho algumas perguntas sobre especificação e interpretação de GLMMs. Definitivamente, três perguntas são estatísticas e duas são mais específicas sobre R. Estou publicando aqui porque, em última análise, acho que a questão é a interpretação dos resultados do GLMM.

Atualmente, estou tentando ajustar um GLMM. Estou usando os dados do censo dos EUA do Longitudinal Tract Database . Minhas observações são setores censitários. Minha variável dependente é o número de unidades habitacionais vagas e estou interessado na relação entre vaga e variáveis ​​socioeconômicas. O exemplo aqui é simples, basta usar dois efeitos fixos: porcentagem da população não branca (raça) e renda média da família (classe), além de sua interação. Eu gostaria de incluir dois efeitos aleatórios aninhados: tratados dentro de décadas e décadas, ou seja (década / setor). Estou considerando essas aleatórias em um esforço para controlar a autocorrelação espacial (isto é, entre setores) e temporal (isto é, entre décadas). No entanto, também estou interessado na década como efeito fixo, por isso também o incluo como fator fixo.

Como minha variável independente é uma variável de contagem inteira não negativa, eu tenho tentado ajustar poisson e GLMMs binomiais negativos. Estou usando o log do total de unidades habitacionais como um deslocamento. Isso significa que os coeficientes são interpretados como o efeito na taxa de desocupação, não no número total de casas desocupadas.

Atualmente, tenho resultados para um Poisson e um GLMM binomial negativo estimado usando glmer e glmer.nb do lme4 . A interpretação dos coeficientes faz sentido para mim com base no meu conhecimento da área de dados e estudo.

Se você deseja que os dados e o script estejam no meu Github . O script inclui mais investigações descritivas que eu fiz antes de construir os modelos.

Aqui estão os meus 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

Testes 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

Testes DHARMa binomiais negativos

    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

DHARMa parcelas

Poisson

Binomial negativo

Questões estatísticas

Como ainda estou descobrindo os GLMMs, estou me sentindo insegura quanto à especificação e interpretação. Eu tenho algumas questões:

1. Parece que meus dados não suportam o uso de um modelo de Poisson e, portanto, estou melhor com o binômio negativo. No entanto, recebo constantemente avisos de que meus modelos binomiais negativos atingem seu limite de iteração, mesmo quando eu aumento o limite máximo. "Em theta.ml (Y, mu, pesos = objeto @ resp $ pesos, limite = limite,: limite de iteração atingido." Isso ocorre usando algumas especificações diferentes (por exemplo, modelos mínimo e máximo para efeitos fixos e aleatórios). Eu também tentei remover discrepantes no meu dependente (bruto, eu sei!), Já que o 1% superior dos valores é muito discrepante (99% inferior variam de 0 a 1012, 1% superior a 1013-5213). não tem nenhum efeito nas iterações e muito pouco efeito nos coeficientes.Não incluo esses detalhes aqui. Os coeficientes entre Poisson e o binômio negativo também são bastante semelhantes. Essa falta de convergência é um problema? O modelo binomial negativo é um bom ajuste? Também executei o modelo binomial negativo usandoO AllFit e nem todos os otimizadores emitem esse aviso (bobyqa, Nelder Mead e nlminbw não).

2. A variação para o efeito fixo da minha década é consistentemente muito baixa ou 0. Entendo que isso pode significar que o modelo está super ajustado. Tirar uma década dos efeitos fixos aumenta a variação do efeito aleatório da década para 0,2620 e não afeta muito os coeficientes de efeitos fixos. Há algo de errado em deixá-lo dentro? Estou bem interpretando-o como simplesmente não sendo necessário para explicar a variação da observação.

3. Esses resultados indicam que eu deveria tentar modelos inflados com zero? A DHARMa parece sugerir que a inflação zero pode não ser o problema. Se você acha que eu deveria tentar de qualquer maneira, veja abaixo.

Perguntas R

4. Eu gostaria de experimentar modelos inflados com zero, mas não tenho certeza de qual pacote implementa efeitos aleatórios aninhados para Poisson inflado com zero e GLMMs binomiais negativos. Eu usaria o glmmADMB para comparar o AIC com modelos inflados por zero, mas ele é restrito a um único efeito aleatório, portanto não funciona para este modelo. Eu poderia tentar o MCMCglmm, mas não conheço as estatísticas bayesianas, o que também não é atraente. Alguma outra opção?

5. Posso exibir coeficientes exponenciados no resumo (modelo) ou preciso fazê-lo fora do resumo, como fiz aqui?

Acredito que existem alguns problemas importantes a serem abordados com sua estimativa.

Pelo que reuni examinando seus dados, suas unidades não estão geograficamente agrupadas, ou seja, setores censitários nos condados. Portanto, o uso de setores como um fator de agrupamento não é apropriado para capturar a heterogeneidade espacial, pois isso significa que você tem o mesmo número de indivíduos que os grupos (ou, em outras palavras, todos os seus grupos têm apenas uma observação cada). O uso de uma estratégia de modelagem multinível nos permite estimar a variação no nível individual, enquanto controlamos a variação no grupo. Como seus grupos têm apenas um indivíduo cada, sua variação entre grupos é a mesma que sua variação no nível individual, derrotando assim o objetivo da abordagem multinível.

Por outro lado, o fator de agrupamento pode representar medições repetidas ao longo do tempo. Por exemplo, no caso de um estudo longitudinal, as pontuações de "matemática" de um indivíduo podem ser registradas anualmente, portanto, teríamos um valor anual para cada aluno por n anos (nesse caso, o fator de agrupamento é o aluno como temos n número de observações "aninhadas" nos alunos). No seu caso, você repetiu as medidas de cada setor censitário até decade. Assim, você pode usar sua TRTID10variável como um fator de agrupamento para capturar "variação entre décadas". Isso leva a 3142 observações aninhadas em 635 setores, com aproximadamente 4 e 5 observações por setor censitário.

Como mencionado em um comentário, o uso decadecomo fator de agrupamento não é muito apropriado, pois você tem apenas cerca de 5 décadas para cada setor censitário, e seu efeito pode ser melhor capturado decadequando apresentado como covariável.

Segundo, para determinar se seus dados devem ser modelados usando um modelo binomial de poisson ou negativo (ou uma abordagem inflada a zero). Considere a quantidade de sobredispersão em seus dados. A característica fundamental de uma distribuição de Poisson é a equidispersão, o que significa que a média é igual à variação da distribuição. Olhando para seus dados, é bastante claro que há muita sobredispersão. As variações são muito maiores do que as médias.

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 entanto, para determinar se o binômio negativo é mais apropriado estatisticamente, um método padrão é fazer um teste de razão de verossimilhança entre um modelo de Poisson e um modelo binomial negativo, o que sugere que o negbin é um ajuste melhor.

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

Após estabelecer isso, um próximo teste poderá considerar se a abordagem multinível (modelo misto) é justificada usando uma abordagem semelhante, o que sugere que a versão multinível fornece um melhor ajuste. (Um teste semelhante pode ser usado para comparar um ajuste glmer assumindo uma distribuição de poisson para o objeto glmer.nb, desde que os modelos sejam os mesmos.)

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

Com relação às estimativas dos modelos de poisson e nb, espera-se que eles sejam muito semelhantes entre si, com a principal distinção os erros padrão, ou seja, se houver excesso de dispersão, o modelo de poisson tende a fornecer erros padrão tendenciosos. Tomando seus dados como exemplo:

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 como as estimativas do coeficiente são todas muito semelhantes, a principal diferença sendo apenas o significado de uma de suas covariáveis, bem como a diferença na variação de efeitos aleatórios, o que sugere que a variação no nível da unidade capturada pelo parâmetro de super-dispersão no nb O modelo (o thetavalor no objeto glmer.nb) captura algumas das variações entre áreas capturadas pelos efeitos aleatórios.

Em relação aos coeficientes exponenciados (e intervalos de confiança associados), você pode usar o seguinte:

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

Considerações finais sobre inflação zero. Não há implementação multinível (pelo menos que eu saiba) de um modelo de poisson ou negbin inflado com zero que permita especificar uma equação para o componente inflado com zero da mistura. o glmmADMBmodelo permite estimar um parâmetro de inflação zero constante. Uma abordagem alternativa seria usar a zeroinflfunção no psclpacote, embora isso não suporte modelos multiníveis. Assim, você pode comparar o ajuste de um binômio negativo de nível único com o binomial negativo inflado com nível único de zero. As chances são de que, se a inflação zero não for significativa para modelos de nível único, é provável que não seja significativa para a especificação multinível.

Termo aditivo

Se você estiver preocupado com a autocorrelação espacial, poderá controlar isso usando alguma forma de regressão geográfica ponderada (embora eu acredite que isso use dados pontuais, não áreas). Como alternativa, você pode agrupar seus setores censitários por um fator de agrupamento adicional (estados, municípios) e incluí-lo como um efeito aleatório. Por fim, e não tenho certeza se isso é totalmente viável, pode ser possível incorporar dependência espacial usando, por exemplo, a contagem média de R_VACvizinhos de primeira ordem como covariável. Em qualquer caso, antes de tais abordagens, seria sensato determinar se a autocorrelação espacial está realmente presente (usando os testes I, LISA da Global Moran e abordagens semelhantes).