19
$\begingroup$

Sorry if I'm missing something very obvious here but I am new to mixed effect modelling.

I am trying to model a binomial presence/absence response as a function of percentages of habitat within the surrounding area. My fixed effect is the percentage of the habitat and my random effect is the site (I mapped 3 different farm sites). 

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

When 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

My output is this:

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

I basically do not understand why my random effect is 0? Is it because the random effect only has 3 levels? I don't see why this would be the case. I have tried it with lots of different models and it always comes out as 0.

It cant be because the random effect doesn't explain any of the variation because I know the habitats are different in the different farms.

Here is an example set of data using 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)

It goes on with around 10 more columns for different types of habitat (like built, setaside is one of them) with percentages in it.

Improve this question
$\endgroup$
3
  • $\begingroup$ dput(territory) gives a lot of text - shall I put it all up here? I will post the results of verbose=TRUE $\endgroup$
    – Cec.g
    Commented Aug 23, 2012 at 20:07
  • 1
    $\begingroup$ This output says that the estimated random effect variance is 0. One way of interpreting this result is that observations of treat with the same value for farm are uncorrelated. Does this result make sense in the context of the application? Some explanation of what the treat variable is would help clarify this. If it's a treatment variable, and treatment was randomly assigned within farm, then the population intraclass correlation would be zero, so we'd expect a small estimate. $\endgroup$
    – Macro
    Commented Aug 23, 2012 at 21:35
  • $\begingroup$ Very related: stats.stackexchange.com/questions/115090. $\endgroup$
    – amoeba
    Commented Jun 21, 2017 at 7:55

2 Answers 2

Reset to default
21
$\begingroup$

With just three farms, there is no point in trying to pretend that you can fit a Gaussian distribution to three points. Analyze this simply as lm(response~as.factor(farm) + treat+other stuff), and don't bother with lmer; you won't be able to do much better than ANOVA, anyway.

Generally, hitting exactly zero is not that unusual. The variance estimate is a nonlinear function of the data, the difference between the overall variance and the within-site variance. If the true variance is zero, this nonlinear statistic has a distribution that puts non-zero mass to the left of zero (this will also be true if the true value is a small positive quantity, but the sampling variability is large enough to overshoot below zero). Due to the way the estimator is programmed, however (Cholesky factorization), it can only take non-negative values. So whenever the unattainably best estimate would have been at zero (as in your balanced-by-design situation) or below it, the log-likelihood will be maximized at zero, with a negative gradient to the right of it. Self & Liang (1987) is the standard biostat reference for the problem; I better like Andrews (1999) which is even more general.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks very much for your help. I think I will do either a binomial response glm or as you say, switch it around and do an ANOVA with an error term for farm. $\endgroup$
    – Cec.g
    Commented Aug 23, 2012 at 23:14
7
$\begingroup$

Looks like there was probably no effect due to Farm built in from the experimental design; each farm has exactly half treated and half not.

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

It can also be instructive to fit farm as a fixed effect and see what happens; we see that the Farm effect is very, very small compared with the built effect, so I wouldn't be too surprised that the fitted variance in the mixed model is zero.

> 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
Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Oh ok, thanks for your help! I think this may be down to my confusion over how to assign the random effects. Treat is rather confusingly named I'm afraid. Treat is whether a territory is occupied or not. So I have territory locations on each farm treat==1 and have selected the same number of random points to act as unoccupied territories treat==0. So when I include farm as a random effect in the model, I'm not grouping the data into the three levels but asking the model if the number of occupied or unnoccupied territories differs per farm? (which by design is exactly the same). $\endgroup$
    – Cec.g
    Commented Aug 23, 2012 at 22:09
  • 2
    $\begingroup$ Yes, that's exactly what you are asking. You should have a real response variable on the LHS (e.g., the number of individuals found on each site, biomass production, whatever... something biological) $\endgroup$
    – StasK
    Commented Aug 23, 2012 at 22:26
  • $\begingroup$ Thanks. I was doing an incidence function response where the incidence is a function of the habitat in the surrounding area. But got very confused with how to combine this with random effects. I dont suppose I can do this in this type of model - and make the variance of built randomly affected by farm without changing the model around so incidence is a function of the surrounding habitat? $\endgroup$
    – Cec.g
    Commented Aug 23, 2012 at 22:51

Not the answer you're looking for? Browse other questions tagged or ask your own question.