Solved – mixed effects and lme4: Do I need nesting

ecologyexperiment-designlme4-nlmenested datarandom-effects-model

I am analyzing data from a field experiment, and I am interested in the effects of fauna and altitude (fixed). Altitude has two levels, and at each site I have 5 blocks for the three levels of fauna (the blocks were set up to account for the heterogeneity in the soil).
I guess that the model treating block as random, in lme4, would be:

lmer(response ~ altitude  + fauna + (1|block))

However, I replicated this study in 5 mountains. So can I treat mountain as another random factor and then:

lmer(response ~ altitude  + fauna + (1|block) + (1|mountain))

Or should I nest my blocks within mountains?

lmer(response ~ altitude  + fauna + (1|mountain:block))

I see people using nested designs when they have pseudoreplicates, but that is not my really my case as all my replicates are independent.

Here's how my data looks like: 

mountain   altitude   treatment block response  
m1         high       t1        b1    124.77  
m1         high       t2        b1    55.77  
m1         high       t3        b1    88.99  
m1         high       t1        b2    88.99  
m1         high       t2        b2    88.99  
m1         high       t3        b2    88.99  
m1         low        t1        b6    124.77  
m1         low        t2        b6    55.77  
m1         low        t3        b6    78.99  
m1         low        t1        b7    89.99  
m1         low        t2        b7    33.99  
m1         low        t3        b7    22.87  
m2
.
.
.

Best Answer

Since the blocks are explicitly nested, the random parts are (1|block) and (1|mountain).

altitude is a fixed term. If you want a single estimate for each level of altitude (except, of course, for the reference level), then just include altitude a single term as you did in your second suggestion above. However, you can let the effect of altitude vary between mountains, ie. a random slope model.

lmer(response ~ fauna + (1|block) + (altitude|mountain))

(You do not mention treatment, but it is part of your dataset so perhaps it should be in the model too?)

I think that a loglikelihood ratio test can be used to determine if the random slope model performs significantly better than the simpler model.

Eg.

fm1 <- lmer(response ~ altitude + fauna + (1|block) + (1|mountain))
fm2 <- lmer(response ~ fauna + (1|block) + (altitude|mountain))
anova(fm1, fm2)
Related Question