Problemas para encontrar um bom modelo adequado para dados de contagem com efeitos mistos - ZINB ou algo mais?

Eu tenho um conjunto de dados muito pequeno sobre a abundância de abelhas solitárias que estou tendo problemas para analisar. São dados de contagem, e quase todas as contagens estão em um tratamento com a maioria dos zeros no outro tratamento. Também existem alguns valores muito altos (um em cada um dos seis locais), portanto, a distribuição das contagens tem uma cauda extremamente longa. Estou trabalhando na R. Eu usei dois pacotes diferentes: lme4 e glmmADMB.

Os modelos mistos de Poisson não se encaixavam: os modelos eram superdispersos quando efeitos aleatórios não eram ajustados (modelo glm) e subdispersos quando efeitos aleatórios eram ajustados (modelo glmer). Eu não entendo por que isso é. O design experimental exige efeitos aleatórios aninhados, por isso preciso incluí-los. Uma distribuição de erro lognormal de Poisson não melhorou o ajuste. Tentei a distribuição negativa de erros binomiais usando glmer.nb e não consegui ajustá-lo - o limite de iteração foi atingido, mesmo quando a tolerância foi alterada usando glmerControl (tolPwrss = 1e-3).

Como muitos dos zeros devem-se ao fato de eu simplesmente não ver as abelhas (geralmente são pequenas coisas negras), tentei um modelo inflado a zero. O ZIP não se encaixou bem. O ZINB foi o melhor modelo até agora, mas ainda não estou muito feliz com o modelo. Estou sem saber o que tentar em seguida. Tentei um modelo de obstáculos, mas não consegui ajustar uma distribuição truncada para resultados diferentes de zero - acho que porque muitos dos zeros estão no tratamento de controle (a mensagem de erro era "Erro no model.frame.default (formula = s.bee ~ tmt + lu +: comprimentos variáveis ​​diferem (encontrado para 'tratamento') ”).

Além disso, acho que a interação que incluí está fazendo algo estranho aos meus dados, pois os coeficientes são irrealisticamente pequenos - embora o modelo que contém a interação tenha sido melhor quando comparamos modelos usando o AICctab no pacote bbmle.

Estou incluindo algum script R que praticamente reproduzirá meu conjunto de dados. As variáveis ​​são as seguintes:

d = data juliana, df = data juliana (como fator), d.sq = df ao quadrado (o número de abelhas aumenta e depois cai durante o verão), st = local, s.bee = contagem de abelhas, tmt = tratamento, lu = tipo de uso da terra, hab = porcentagem de habitat semi-natural na paisagem circundante, ba = área limite dos campos redondos.

Todas as sugestões sobre como eu posso obter um bom ajuste do modelo (distribuições de erro alternativas, diferentes tipos de modelo etc.) serão muito bem recebidas!

Obrigado.

d <- c(80,  80,  121, 121, 180, 180, 86,  86,  116, 116, 144, 144, 74,  74, 143, 143, 163, 163, 71, 71,106, 106, 135, 135, 162, 162, 185, 185, 83,  83,  111, 111, 133, 133, 175, 175, 85,  85,  112, 112,137, 137, 168, 168, 186, 186, 64,  64,  95,  95,  127, 127, 156, 156, 175, 175, 91,  91, 119, 119,120, 120, 148, 148, 56, 56)
df <- as.factor(d)
d.sq <- d^2
st <- factor(rep(c("A", "B", "C", "D", "E", "F"), c(6,12,18,10,14,6)))
s.bee <- c(1,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,4,0,0,0,0,1,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,3,0,0,0,0,5,0,0,2,0,50,0,10,0,4,0,47,3)
tmt <- factor(c("AF","C","C","AF","AF","C","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","C","AF","AF","C","AF","C","AF","C","AF","C","AF","C","AF","C",
"C","AF","AF","C","AF","C","AF","C","AF","C","C","AF","C","AF","C","AF","AF","C","AF","C",
"AF","C","AF","C","AF","C"))
lu <- factor(rep(c("p","a","p","a","p"), c(6,12,28,14,6)))
hab <- rep(c(13,14,13,14,3,4,3,4,3,4,3,4,3,4,15,35,37,35,37,35,37,35,37,0,2,1,2,1,2,1), 
        c(1,2,2,1,1,1,1,2,2,1,1,1,1,1,18,1,1,1,2,2,1,1,1,14,1,1,1,1,1,1))
ba <-  c(480,6520,6520,480,480,6520,855,1603,855,1603,1603,855,855,12526,855,5100,855,5100,2670,7679,7679,2670,
2670,7679,2670,7679,7679,2670,2670,7679,2670,7679,2670,7679,2670,7679,1595,3000,1595,3000,3000,1595,1595,3000,1595
,3000,4860,5460,4860,5460,5460,4860,5460,4860,5460,4860,4840,5460,4840,5460,3000,1410,3000,1410,3000,1410)
data <- data.frame(st,df,d.sq,tmt,lu,hab,ba,s.bee)
with(data, table(s.bee, tmt) )

# below is a much abbreviated summary of attempted models:

library(MASS)
library(lme4)
library(glmmADMB)
library(coefplot2)

###
### POISSON MIXED MODEL

    m1 <- glmer(s.bee ~ tmt + lu + hab + (1|st/df), family=poisson)
    summary(m1)

    resdev<-sum(resid(m1)^2)
    mdf<-length(fixef(m1)) 
    rdf<-nrow(data)-mdf
    resdev/rdf
# 0.2439303
# underdispersed. ???



###
### NEGATIVE BINOMIAL MIXED MODEL

    m2 <- glmer.nb(s.bee ~ tmt + lu + hab + d.sq + (1|st/df))
# iteration limit reached. Can't make a model work.



###
### ZERO-INFLATED POISSON MIXED MODEL

    fit_zipoiss <- glmmadmb(s.bee~tmt + lu + hab + ba + d.sq +
                    tmt:lu +
                    (1|st/df), data=data,
                    zeroInflation=TRUE,
                    family="poisson")
# has to have lots of variables to fit
# anyway Poisson is not a good fit



###
### ZERO-INFLATED NEGATIVE BINOMIAL MIXED MODELS

## BEST FITTING MODEL SO FAR:

    fit_zinb <- glmmadmb(s.bee~tmt + lu + hab +
                    tmt:lu +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")
    summary(fit_zinb)
# coefficients are tiny, something odd going on with the interaction term
# but this was best model in AICctab comparison

# model check plots
    qqnorm(resid(fit_zinb))
    qqline(resid(fit_zinb))

    coefplot2(fit_zinb)

    resid_zinb <- resid(fit_zinb , type = "pearson")
    hist(resid_zinb)

    fitted_zinb <- fitted (fit_zinb)
    plot(resid_zinb ~ fitted_zinb)


## MODEL WITHOUT INTERACTION TERM - the coefficients are more realistic:

    fit_zinb2 <- glmmadmb(s.bee~tmt + lu + hab +
                    (1|st/df),data=data,
                    zeroInflation=TRUE,
                    family="nbinom")

# model check plots
    qqnorm(resid(fit_zinb2))
    qqline(resid(fit_zinb2))

    coefplot2(fit_zinb2)

    resid_zinb2 <- resid(fit_zinb2 , type = "pearson")
    hist(resid_zinb2)

    fitted_zinb2 <- fitted (fit_zinb2)
    plot(resid_zinb2 ~ fitted_zinb2)



# ZINB models are best so far
# but I'm not happy with the model check plots

Este post tem quatro anos, mas eu queria seguir o que a fsociety disse em um comentário. O diagnóstico de resíduos nos GLMMs não é simples, pois os gráficos residuais padrão podem mostrar não normalidade, heterocedasticidade etc., mesmo que o modelo esteja especificado corretamente. Existe um pacote R DHARMa, especificamente adequado para diagnosticar esse tipo de modelo.

O pacote é baseado em uma abordagem de simulação para gerar resíduos em escala a partir de modelos mistos lineares generalizados ajustados e gera diferentes gráficos de diagnóstico facilmente interpretáveis. Aqui está um pequeno exemplo com os dados da postagem original e do primeiro modelo ajustado (m1):

library(DHARMa)
sim_residuals <- simulateResiduals(m1, 1000)
plotSimulatedResiduals(sim_residuals)

O gráfico à esquerda mostra um gráfico QQ dos resíduos dimensionados para detectar desvios da distribuição esperada, e o gráfico à direita representa resíduos versus valores previstos, enquanto executa regressão quantílica para detectar desvios da uniformidade (as linhas vermelhas devem ser horizontais e a 0,25 0,50 e 0,75).

Além disso, o pacote também possui funções específicas para testar dispersão acima / abaixo e inflação zero, entre outras:

testOverdispersionParametric(m1)

Chisq test for overdispersion in GLMMs

data:  poisson
dispersion = 0.18926, pearSS = 11.35600, rdf = 60.00000, p-value = 1
alternative hypothesis: true dispersion greater 1

testZeroInflation(sim_residuals)

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


data:  sim_residuals
ratioObsExp = 0.98894, p-value = 0.502
alternative hypothesis: more