GLMM para datos de recuento utilizando el enlace de la raíz cuadrada en lme4

Question

Tengo datos de un estudio de campo. El objetivo del estudio es relacionar el número de plántulas (variable de respuesta, datos de recuento), la forma del terreno (variable exploratoria, variable categórica con 3 niveles) y el porcentaje de cobertura del dosel (variable exploratoria, cuantitativa). En cada hábitat, tengo datos de cinco parcelas de 25x25 metros. Dentro de cada parcela utilicé tres subparcelas de 2x2 metros anidadas dentro de la parcela mayor, y el número de plántulas se contó a partir de estas subparcelas. El número total de observaciones es de 60; 20 parcelas x 3 subparcelas. Sólo un tipo de terreno contenía plántulas. Las otras dos formas de terreno no contenían plántulas, por lo que son parcelas vacías:

data.frame':    60 obs. of  5 variables:
 $ plot         : Factor w/ 20 levels "HD2","LC1","LC2",..: 16 16 16 17 17 17 12 12 12 9 ...
 $ subplot      : Factor w/ 60 levels "HD2.1","HD2.2",..: 46 47 48 49 50 51 34 35 36 25 ...
 $ av.canopy    : num  92.2 92.2 92.2 92.3 92.3 ...
 $ landform     : Factor w/ 3 levels "abandoned","ridge",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ seedling.2016: int  6 7 5 2 5 4 4 6 4 0 ...

El problema es cuando comparé el número de plántulas entre formas de terreno utilizando este código:

seedling <- glmer(seedling.2016 ~ landform +(1|plot), family = poisson)

El resultado no tiene sentido para mí: no hay ninguna diferencia significativa entre las formas del terreno, ya que sólo hay una forma del terreno (cresta) que tiene plántulas, mientras que otras no tienen plántulas. El resultado es el siguiente:

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: poisson  ( log )
Formula: seedling.2016 ~ landform + (1 | plot)
   Data: streblus.subplots

     AIC      BIC   logLik deviance df.resid 
   294.9    303.3   -143.5    286.9       56 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-6.3619 -0.0001 -0.0001  0.0000  8.7305 

Random effects:
 Groups Name        Variance Std.Dev.
 plot   (Intercept) 2.637    1.624   
Number of obs: 60, groups:  plot, 20

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)    -20.412   1461.267  -0.014    0.989
landformridge   22.250   1461.265   0.015    0.988
landformtemp     1.066    390.540   0.003    0.998

Cuando cambié la función de enlace a la raíz cuadrada como este código:

Seedling2 <- glmer(seedling.2016 ~ landform +  (1|plot), family = poisson(link = sqrt))
Fixed effects:
#Estimate      Std. Error z value Pr(>|z|)    
#(Intercept)   -1.220e-05  5.296e-01   0.000        1    
#landformridge  3.250e+00  7.429e-01   4.376 1.21e-05 ***
# landformtemp   1.018e-05  7.795e-01   0.000        1  

Ahora el número de plántulas en la cresta es significativamente mayor que el otro, y tiene más sentido para mí.

Mi pregunta es: ¿Es válido en términos de estadística utilizar el enlace de la raíz cuadrada con la distribución de Poisson, hay algún método mejor disponible con mejor base de estadística?

Answer 1

Parece que tienes un caso de separación completa :

• sólo hay una forma de terreno (cresta) que tiene plántulas, mientras que otras no tienen plántulas al

• grandes estimaciones ( $|\hat \beta|>10$ ), y estimaciones de error estándar ridículamente grandes.

Básicamente, lo que ocurre es que el nivel de referencia ("abandonado") tiene un número esperado de recuentos igual a cero para todas las parcelas, por lo que el intercepto $\beta_0$ - que es el log(recuento) esperado para el nivel de referencia - debe estimarse como $-\infty$ ... que estropea la estimación de Wald de la incertidumbre (el método aproximado y rápido que summary() usos).

Puede leer más sobre la separación completa en otro lugar; se suele hablar de ella en el contexto de regresión logística (en parte porque la regresión logística se utiliza más que la regresión de recuento...)

Soluciones:

• su solución de enlace de raíz cuadrada es razonable (en este caso el intercepto se espera $\sqrt{\textrm{counts}}$ en el nivel de referencia, que es cero en lugar de $-\infty$ ); cambiará ligeramente la distribución asumida de los efectos aleatorios (es decir, Normal en la raíz cuadrada en lugar de en la escala logarítmica), pero eso no me preocuparía demasiado. Si tuviera covariables continuas o interacciones en el modelo, también cambiaría la interpretación de los efectos fijos.
• se podría utilizar algún tipo de penalización (muy convenientemente en un marco bayesiano), como se describe en mi respuesta a la pregunta enlazada (y aquí buscar "separación completa") para que los parámetros sean razonables.

Answer 2

De hecho, parece que se trata de un problema de separación. Para tener en cuenta estos casos, en mi GLMMadaptive puede incluir una penalización para los coeficientes de efectos fijos en forma de densidad Students-t (es decir, para un df suficientemente grande, equivalente a la regresión ridge). Para ver un ejemplo práctico, consulte la última sección de esta viñeta .