Resolvendo a heterocedasticidade no Poisson GLMM

8

Eu tenho dados de coleta de longo prazo e gostaria de testar se o número de animais coletados é influenciado pelos efeitos climáticos. Meu modelo se parece abaixo:

glmer(SumOfCatch ~ I(pc.act.1^2) +I(pc.act.2^2) + I(pc.may.1^2) + I(pc.may.2^2) + 
                   SampSize + as.factor(samp.prog) + (1|year/month), 
      control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=1e9,npt=5)), 
      family="poisson", data=a2)

Explicação das variáveis ​​utilizadas:

  • SumOfCatch: número de animais coletados
  • pc.act.1, pc.act.2: eixos de um componente principal que representam as condições climáticas durante a amostragem
  • pc.may.1, pc.may.2: eixos de um PC representando as condições meteorológicas em maio
  • SampSize: número de armadilhas de interceptação ou coleta de transectos de comprimentos padrão
  • samp.prog: método de amostragem
  • ano: ano da amostragem (de 1993 a 2002)
  • mês: mês da amostragem (de agosto a novembro)

Os resíduos do modelo ajustado mostram considerável heterogeneidade (heterocedasticidade?) Quando plotados em relação aos valores ajustados (ver Fig.1):

Resíduos vs. valores ajustados - modelo principal

Minha principal pergunta é: esse é um problema que torna questionável a confiabilidade do meu modelo? Nesse caso, o que posso fazer para resolvê-lo?

Até agora, tentei o seguinte:

  • controle de sobredispersão, definindo efeitos aleatórios no nível de observação, ou seja, usando um ID exclusivo para cada observação e aplicando essa variável de ID como efeito aleatório; embora meus dados mostrem uma excessiva dispersão considerável, isso não ajudou, pois os resíduos se tornaram ainda mais feios (veja a Fig. 2)

Resíduos vs. valores ajustados - modelo com controle de DO

  • Eu ajustei modelos sem efeitos aleatórios, com quase-Poisson glm e glm.nb; também produziu parcelas residuais vs. ajustadas similares ao modelo original

Até onde eu sei, pode haver maneiras de estimar erros padrão consistentes em heterocedasticidade, mas não consegui encontrar nenhum método para os GLMMs de Poisson (ou qualquer outro tipo de) em R.


Em resposta a @FlorianHartig: o número de observações no meu conjunto de dados é N = 554, acho que isso é uma obs. número para esse modelo, mas, é claro, quanto mais, melhor. Postei duas figuras, a primeira delas é a plotagem residual em escala DHARMa (sugerida por Florian) do modelo principal.

insira a descrição da imagem aqui

A segunda figura é de um segundo modelo, no qual a única diferença é que ele contém o efeito aleatório no nível de observação (o primeiro não).

insira a descrição da imagem aqui

ATUALIZAR

Figura da relação entre uma variável climática (como preditor, ou seja, eixo x) e sucesso da amostragem (resposta):

Weather-PC e sucesso de amostragem

ATUALIZAÇÃO II.

Figuras mostrando valores preditivos vs. resíduos:

Preditores vs. Residuais

Z. Radai
fonte
Você já pensou em executar um estimador não paramétrico? Ou comparando um ols a uma regressão mediana? Eu percebo que Poisson é o modelo dominante em bio, mas GLMs são inconsistentes sob heterocedasticidade e OLS não é.
Superpronker
1
Às vezes, a super-dispersão é causada pela inflação zero. Nesse caso, você pode tentar um modelo de poisson com parâmetro de inflação zero ou um modelo de barreira. O pacote glmmADMB possui ótimos recursos para lidar com isso: glmmadmb.r-forge.r-project.org/glmmADMB.html
Niek
Caro @Superpronker, obrigado pela sugestão, não verifiquei o OLS, não sabia que essa abordagem seria flexível o suficiente para lidar com meus dados. Vou dar uma olhada nele
Z. Radai
Caro @Niek nos meus dados, não há observações de zero - caso contrário, pensei nos modelos zeroinfl e hurdle (no pacote 'pscl') devido ao seu bom manuseio de sobredispersão, mas eles só são utilizáveis ​​em dados com zeros na resposta . Há alguns meses, tentei o glmmADMB, mas não produziu melhores resultados. Cheers, ZR
Z. Radai
1
@mdewey A lógica por trás disso é que a relação entre efeitos climáticos e sucesso na amostragem segue ótima: a probabilidade e a extensão do sucesso na amostragem são mais altas em um determinado intervalo (neste caso, zero e sua vizinhança) do preditor. Quando o valor do preditor estiver mais distante desse ótimo, o sucesso da amostragem será menor, correspondendo a um subóptimo. Eu uso o termo quadrático, porque (1) os preditores são redimensionados e atualizados em zero e (2) isso fornece uma melhor aproximação a uma conexão linear. Cheers, ZR
Z. Radai

Respostas:

9

É difícil avaliar o ajuste do Poisson (ou qualquer outro GLM com valor inteiro) com resíduos de Pearson ou desvio, porque também um Poisson GLMM perfeitamente adequado exibirá resíduos de desvio não homogêneos.

Isso é especialmente verdade se você fizer GLMMs com REs no nível de observação, porque a dispersão criada pelos OL-REs não é considerada pelos resíduos de Pearson.

Para demonstrar o problema, o código a seguir cria dados Poisson superdispersos, que são então equipados com um modelo perfeito. Os resíduos de Pearson se parecem muito com o seu enredo - portanto, pode ser que não haja nenhum problema.

Esse problema é resolvido pelo pacote DHARMa R , que simula a partir do modelo ajustado para transformar os resíduos de qualquer GL (M) M em um espaço padronizado. Feito isso, é possível avaliar / testar visualmente problemas residuais, como desvios da distribuição, dependência residual de um preditor, heterocedasticidade ou autocorrelação da maneira normal. Veja a vinheta do pacote para obter exemplos detalhados. Você pode ver no gráfico inferior que o mesmo modelo agora parece bom, como deveria.

Se você ainda perceber heterocedasticidade após plotar com o DHARMa, terá que modelar a dispersão em função de algo, o que não é um grande problema, mas provavelmente exigiria a mudança para JAGs ou outro software bayesiano.

library(DHARMa)
library(lme4)

testData = createData(sampleSize = 200, overdispersion = 1, randomEffectVariance = 1, family = poisson())

fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) + (1|ID), family = "poisson", data = testData, control=glmerControl(optCtrl=list(maxfun=20000) ))

# standard Pearson residuals
plot(fittedModel, resid(., type = "pearson") ~ fitted(.) , abline = 0)

# DHARMa residuals
plot(simulateResiduals(fittedModel))

insira a descrição da imagem aqui

insira a descrição da imagem aqui

Florian Hartig
fonte
Dear @FlorianHartig! Obrigado por suas idéias, tentei tramar com a DHARMa. Com base no gráfico, ainda há algo acontecendo, fazendo com que o quantil inferior seja modelado como uma curva recíproca, em vez de uma linha reta. Você mencionou que, nesse caso, uma solução pode ser modelar a dispersão em função de algo - você poderia ajudar exatamente como posso avaliar essa função? Cheers, ZR
Z. Radai
Você pode me enviar ou publicar a trama? Alguma variabilidade menor é esperado quando o N é pequeno
Florian Hartig
Caro @FlorianHartig, a pergunta foi editada, agora também mostrando os gráficos da DHARMa!
Z. Radai
@ Z.Radai - o que vejo nas plotagens é que seus resíduos são sistematicamente altos demais para previsões de modelo baixas. Isso parece mais um problema de estrutura do modelo (falta de preditores?) Do que um problema de distribuição - eu tentaria plotar os resíduos em relação a possíveis e possivelmente ausentes preditores.
Florian Hartig
1
Eu não me preocuparia com a heterocedasticidade, no seu caso é moderada e o efeito na inferência deve ser leve - a única questão que vejo é a subpredição sistemática de valores pequenos, que não serão resolvidos modelando a variância. Mas se você precisa saber, veja aqui stats.stackexchange.com/questions/247183/…
Florian Hartig