Efecto aleatorio igual a 0 en el modelo lineal mixto generalizado

Question

Lo siento si me estoy perdiendo algo muy obvio aquí, pero soy nuevo en el modelado de efectos mixtos.

Estoy tratando de modelar una respuesta binomial de presencia/ausencia en función de los porcentajes de hábitat dentro del área circundante. Mi efecto fijo es el porcentaje del hábitat y mi efecto aleatorio es el sitio (he mapeado 3 sitios de granjas diferentes).

glmmsetaside <- glmer(treat~setas+(1|farm),
       family=binomial,data=territory)

En verbose=TRUE :

0:     101.32427: 0.333333 -0.0485387 0.138083 
1:     99.797113: 0.000000 -0.0531503 0.148455  
2:     99.797093: 0.000000 -0.0520462 0.148285  
3:     99.797079: 0.000000 -0.0522062 0.147179  
4:     99.797051: 7.27111e-007 -0.0508770 0.145384  
5:     99.797012: 1.45988e-006 -0.0495767 0.141109  
6:     99.797006: 0.000000 -0.0481233 0.136883  
7:     99.797005: 0.000000 -0.0485380 0.138081  
8:     99.797005: 0.000000 -0.0485387 0.138083  

Mi salida es esta:

Generalized linear mixed model fit by the Laplace approximation 
Formula: treat ~ setasidetrans + (1 | farm) 

AIC   BIC logLik deviance
105.8 112.6  -49.9     99.8
Random effects:
 Groups Name        Variance Std.Dev.
farm   (Intercept)  0        0  
Number of obs: 72, groups: farm, 3

Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.04854    0.44848  -0.108    0.914
setasidetrans  0.13800    1.08539   0.127    0.899

Correlation of Fixed Effects:
            (Intr)
setasidtrns -0.851

Básicamente, no entiendo por qué mi efecto aleatorio es 0. ¿Es porque el efecto aleatorio sólo tiene 3 niveles? No veo por qué este sería el caso. Lo he probado con muchos modelos diferentes y siempre sale 0.

No puede ser porque el efecto aleatorio no explica nada de la variación porque sé que los hábitats son diferentes en las distintas granjas.

A continuación se muestra un conjunto de datos de ejemplo utilizando dput :

list(territory = c(1, 7, 8, 9, 10, 11, 12, 13, 14, 2, 3, 4, 5, 
6, 15, 21, 22, 23, 24, 25, 26, 27, 28, 16, 17, 18, 19, 20, 29, 
33, 34, 35, 36, 37, 38, 39, 40, 30, 31, 32, 41, 45, 46, 47, 48, 
49, 50, 51, 52, 42, 43, 44, 53, 55, 56, 57, 58, 59, 60, 61, 62, 
54, 63, 65, 66, 67, 68, 69, 70, 71, 72, 64), treat = c(1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 
0, 0, 0, 0, 0, 0, 0), farm = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3), 
built = c(5.202332763, 1.445026852, 2.613422283, 2.261705833, 
2.168842186, 1.267473928, 0, 0, 0, 9.362387965, 17.55433115, 
4.58020626, 4.739300829, 8.638442377, 0, 1.220760647, 7.979990338, 
13.30789514, 0, 8.685544976, 3.71617163, 0, 0, 6.802926951, 
8.925512803, 8.834006678, 4.687723044, 9.878232478, 8.097800267, 
0, 0, 0, 0, 5.639651095, 9.381654651, 8.801754791, 5.692392532, 
3.865304919, 4.493438554, 4.826277798, 3.650995554, 8.20818417, 
0, 8.169597157, 8.62030666, 8.159474015, 8.608979238, 0, 
8.588288678, 7.185700856, 0, 0, 3.089524893, 3.840381223, 
31.98103158, 5.735501995, 5.297691011, 5.17141191, 6.007539933, 
2.703345394, 4.298077606, 1.469986793, 0, 4.258511595, 0, 
21.07029581, 6.737664009, 14.36176373, 3.056631919, 0, 32.49289428, 
0)

Continúa con unas 10 columnas más para diferentes tipos de hábitat (como built , setaside es uno de ellos) con porcentajes.

Answer 1

Con sólo tres explotaciones, no tiene sentido pretender que se puede ajustar una distribución gaussiana a tres puntos. Analícelo simplemente como lm(response~as.factor(farm) + treat+other stuff) y no te molestes en lmer de todos modos, no podrá hacerlo mucho mejor que ANOVA.

Por lo general, llegar exactamente a cero no es tan inusual. La estimación de la varianza es una función no lineal de los datos, la diferencia entre la varianza global y la varianza dentro del sitio. Si la varianza real es cero, este estadístico no lineal tiene una distribución que sitúa la masa distinta de cero a la izquierda de cero (esto también será cierto si el valor real es una cantidad positiva pequeña, pero la variabilidad del muestreo es lo suficientemente grande como para sobrepasar por debajo de cero). Sin embargo, debido a la forma en que está programado el estimador (factorización de Cholesky), sólo puede tomar valores no negativos. Por tanto, siempre que la mejor estimación inalcanzable hubiera estado en cero (como en su situación de equilibrio por diseño) o por debajo de él, la log-verosimilitud se maximizará en cero, con un gradiente negativo a la derecha de él. Self y Liang (1987) es la referencia biostática estándar para el problema; me gusta más Andrews (1999) que es aún más general.

Answer 2

Parece que probablemente no hubo efecto debido a la Granja incorporada desde el diseño experimental; cada granja tiene exactamente la mitad tratada y la otra mitad no.

> xtabs(~treat+farm, territory)
     farm
treat  1  2  3
    0 14 12 10
    1 14 12 10

También puede ser instructivo ajustar Granja como efecto fijo y ver qué ocurre; vemos que el efecto Granja es muy, muy pequeño comparado con el efecto construido, así que no me sorprendería demasiado que la varianza ajustada en el modelo mixto fuera cero.

> m2<-glm(treat~built+factor(farm),family=binomial,data=territory)
> library(car)
> Anova(m2)

Analysis of Deviance Table (Type II tests)

Response: treat
             LR Chisq Df Pr(>Chisq)
built         0.50685  1     0.4765
factor(farm)  0.02008  2     0.9900