Resolviendo la heterocedasticidad en Poisson GLMM

8

Tengo datos de recolección a largo plazo y me gustaría probar si el número de animales recolectados está influenciado por los efectos del clima. Mi modelo se ve a continuación:

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)

Explicación de las variables utilizadas:

  • SumOfCatch: número de animales recolectados
  • pc.act.1, pc.act.2: ejes de un componente principal que representan las condiciones climáticas durante el muestreo
  • pc.may.1, pc.may.2: ejes de una PC que representan las condiciones climáticas en mayo
  • SampSize: número de trampas de caída o recolección de transectos de longitudes estándar
  • samp.prog: método de muestreo
  • año: año de muestreo (de 1993 a 2002)
  • mes: mes de muestreo (de agosto a noviembre)

Los residuos del modelo ajustado muestran una considerable falta de homogeneidad (¿heterocedasticidad?) Cuando se grafican contra los valores ajustados (ver Fig.1):

Residuales frente a valores ajustados - Modelo principal

Mi pregunta principal es: ¿es este un problema que hace que la confiabilidad de mi modelo sea cuestionable? Si es así, ¿qué puedo hacer para resolverlo?

Hasta ahora he intentado lo siguiente:

  • controlar la sobredispersión definiendo efectos aleatorios a nivel de observación, es decir, utilizando una identificación única para cada observación y aplicando esta variable de identificación como efecto aleatorio; Aunque mis datos muestran una sobredispersión considerable, esto no ayudó ya que los residuos se volvieron aún más feos (ver Fig. 2)

Residuos frente a valores ajustados: modelo con control OD

  • Ajusté modelos sin efectos aleatorios, con cuasi-Poisson glm y glm.nb; También arrojó parcelas residuales versus ajustadas similares al modelo original

Hasta donde sé, puede haber formas de estimar los errores estándar consistentes con la heterocedasticidad, pero no he podido encontrar ningún método para GLMM de Poisson (o cualquier otro tipo de) en R.


En respuesta a @FlorianHartig: el número de observaciones en mi conjunto de datos es N = 554, creo que esto es bastante obs. número para tal modelo, pero, por supuesto, cuantos más mejor. Publico dos figuras, la primera de las cuales es la gráfica residual escalada de DHARMa (sugerida por Florian) del modelo principal.

ingrese la descripción de la imagen aquí

La segunda figura es de un segundo modelo, en el que la única diferencia es que contiene el efecto aleatorio a nivel de observación (el primero no).

ingrese la descripción de la imagen aquí

ACTUALIZAR

Figura de la relación entre una variable meteorológica (como predictor, es decir, eje x) y el éxito del muestreo (respuesta):

Weather-PC y muestreo exitoso

ACTUALIZACIÓN II.

Figuras que muestran valores predictores vs. residuales:

Predictores vs. Residuos

Z. Radai
fuente
¿Has considerado ejecutar un estimador no paramétrico? ¿O comparar un ols con una regresión mediana? Me doy cuenta de que Poisson es el modelo dominante en bio, pero los GLM son inconsistentes bajo la heterocedasticidad y OLS no lo es.
Superpronker
1
A veces, la sobredispersión es causada por la inflación cero. En ese caso, puede probar un modelo de Poisson con parámetro de inflación cero o un modelo de obstáculo. El paquete glmmADMB tiene excelentes características para lidiar con esto: glmmadmb.r-forge.r-project.org/glmmADMB.html
Niek
Estimado @Superpronker, gracias por la sugerencia, no verifiqué OLS, no me di cuenta de que este enfoque sería lo suficientemente flexible como para manejar mis datos. Lo investigaré
Z. Radai
Estimado @Niek en mis datos, no hay observaciones de cero; de lo contrario, pensé en zeroinfl y los modelos de obstáculo (en el paquete 'pscl') debido a su buen manejo de la sobredispersión, pero solo se pueden usar en datos con ceros en la respuesta . Hace unos meses intenté glmmADMB, pero no arrojó mejores resultados. Saludos, ZR
Z. Radai
1
@mdewey la razón detrás de esto, es que la relación entre los efectos del clima y el éxito del muestreo sigue un óptimo: la probabilidad y el alcance del éxito del muestreo es mayor en un rango dado (en este caso, cero y su proximidad) del predictor. Cuando el valor del predictor está más alejado de este óptimo, el éxito del muestreo será menor, lo que corresponde a un subóptimo. Utilizo el término cuadrático, porque (1) los predictores se reescalan y se vuelven a centrar en cero, y (2) esto proporciona una mejor aproximación a una conexión lineal. Saludos, ZR
Z. Radai

Respuestas:

9

Es difícil evaluar el ajuste del Poisson (o cualquier otro GLM con valor entero para ese caso) con los residuos de Pearson o de desviación, porque también un GLMM de Poisson perfectamente ajustado exhibirá residuos de desviación no homogéneos.

Esto es especialmente cierto si realiza GLMM con REs de nivel de observación, porque la dispersión creada por OL-REs no es considerada por los residuos de Pearson.

Para demostrar el problema, el siguiente código crea datos de Poisson sobredispersos, que luego se ajustan con un modelo perfecto. Los residuos de Pearson se parecen mucho a su trama, por lo tanto, puede ser que no haya ningún problema.

Este problema se resuelve con el paquete DHARMa R , que simula a partir del modelo ajustado para transformar los residuos de cualquier GL (M) M en un espacio estandarizado. Una vez hecho esto, puede evaluar / evaluar visualmente problemas residuales, como desviaciones de la distribución, dependencia residual de un predictor, heterocedasticidad o autocorrelación de la manera normal. Consulte la viñeta del paquete para ver ejemplos prácticos. Puede ver en la gráfica inferior que el mismo modelo ahora se ve bien, como debería.

Si todavía ve heteroscedasticidad después de trazar con DHARMa, tendrá que modelar la dispersión en función de algo, lo cual no es un gran problema, pero es probable que requiera pasar a JAG u otro 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))

ingrese la descripción de la imagen aquí

ingrese la descripción de la imagen aquí

Florian Hartig
fuente
Estimado @FlorianHartig! Gracias por sus ideas, traté de planear con DHARMa. Según la gráfica, todavía hay algo que hace que el cuantil inferior tenga la forma de una curva recíproca, en lugar de una línea recta. Usted ha mencionado que en este caso, una solución podría ser modelar la dispersión en función de algo. ¿Podría ayudarme exactamente cómo puedo evaluar esa función? Saludos, ZR
Z. Radai
¿Me pueden enviar o publicar la trama? Se espera alguna variabilidad menor cuando su N es pequeña
Florian Hartig
Estimado @FlorianHartig, la pregunta fue editada, ¡ahora también muestra las parcelas de DHARMa!
Z. Radai
@ Z.Radai: lo que veo en las parcelas es que sus residuos son sistemáticamente demasiado altos para predicciones de modelo bajas. Esto se parece más a un problema de la estructura del modelo (¿predictores faltantes?) Que a un problema de distribución: trataría de graficar los residuos contra predictores posibles y potencialmente faltantes.
Florian Hartig
1
No me preocuparía la heterocedasticidad, en su caso es moderada y el efecto sobre la inferencia debería ser leve: el único problema que veo es la predicción sistemática de valores pequeños, que no se resolverá modelando la varianza. Pero si debe saberlo, vea aquí stats.stackexchange.com/questions/247183/…
Florian Hartig el