I want to apply an ANOVA analysis to my nested data. My data structure is as follows:

I have 2 sites (Bruuk en Ketelbroek). Each site has 2 types of habitat; a main habitat (Struweel for Ketelbroek and Nat_bos for Bruuk) and a secondary habitat (grassland for both). In the main habitat I randomly deployed 2 pitfall series, existing of 5 pitfalls each. In the secondary habitat I deployed 1 series of 5 pitfalls, also randomly. The 5 individual pitfalls per series were not randomly placed, but according a systematic design (like the number 5 on a dice). The data looks as follows:

pitfall series  habitat         site            Rich
1       1       Struweel        Ketelbroek      6.628
2       1       Struweel        Ketelbroek      4.959
3       1       Struweel        Ketelbroek      7.205
4       1       Struweel        Ketelbroek      5.629
5       1       Struweel        Ketelbroek      5.793
6       2       Grasland        Ketelbroek      8.195
7       2       Grasland        Ketelbroek      6.917
8       2       Grasland        Ketelbroek      7.803
9       2       Grasland        Ketelbroek      6.395
10      2       Grasland        Ketelbroek      7.147
11      3       Struweel        Ketelbroek      7.011
12      3       Struweel        Ketelbroek      5.959
13      3       Struweel        Ketelbroek      6
14      3       Struweel        Ketelbroek      7.966
15      3       Struweel        Ketelbroek      6.458
16      4       Nat_bos         Bruuk           6.372
17      4       Nat_bos         Bruuk           6.393
18      4       Nat_bos         Bruuk           5.031
19      4       Nat_bos         Bruuk           5.735
20      4       Nat_bos         Bruuk           4.679
21      5       Grasland        Bruuk           5.372
22      5       Grasland        Bruuk           4.917
23      5       Grasland        Bruuk           6
24      5       Grasland        Bruuk           5.713
25      5       Grasland        Bruuk           5.359
26      6       Nat_bos         Bruuk           5.24
27      6       Nat_bos         Bruuk           5.992
28      6       Nat_bos         Bruuk           5.966
29      6       Nat_bos         Bruuk           6.924
30      6       Nat_bos         Bruuk           4.574

I want to test for the effect of the site and the habitat on the rarefied richness ('Rich') of carabid beetles that I encountered in the pitfalls. I know the setup is unbalanced (because the secondary habitat has only 1 series and not 2 like the main habitat), but I would still like to go with a mixed effect model. I was thinking of the following formula:

lmer(richness ~ site + habitat +(1|habitat/series),data=dat)

I am rather uncertain though about the random part


I already looked at the many related posts on random effect structuring of mixed models, but I could not figure what was best, since in my case only the series were randomly selected within each habitat type, and the pitfalls within each series were not. Also the two sites share their secondary habitat type (grasland) but are different in their main habitat. I am not even sure if I should apply a random effect at all, given such few series per habitat (according to the GLMM FAQ one should at least have 5 to 6 repetitions for a random effect).

What is the best way to model the effect of the site (my main interest) and of habitat (my second interest)? Thank you!


1 Answer 1


I would agree with you that modelling a random effect using just 2 levels would be difficult. I would consider using a strict fixed effect model.

Because you have partial nesting of your factors, you could consider creating a third variable that is a combination of the other two, then use that third variable in your model. Your third variable combo would consist of the following four levels: {Struweel_Ketelbroek, Grasland_Ketelbroek, Nat_bos_Bruuk, Grasland_Bruuk}. Then your model would be lm(Rich ~ combo, data=dat), which is the classical ANOVA setup. Since you don't have a natural baseline, you could use contrasts to compare the changes in richness for each level of the combo variable from the mean of all levels (using a sum-to-zero contrasts, or #3 here). For interpretation, you'd interpret the changes of the mean of each level with the grand mean of all levels, or combination of the levels if you want a broader site-to-site comparison or a habitat-to-habitat comparison.

EDIT: Example

Here's an example of defining contrasts for combinations of potential interest. Note that the order of the levels is the same as defined above. I follow the procedure outlined here. First, the full code:

#Read in dat from dataset in question
dat$combo <- paste(dat$habitat, dat$site, sep='_')
dat$combo <- factor(dat$combo, levels=c(unique(dat$combo)))

#Set contrasts of interest
mat <- matrix(c(1/4, 1/4, 1/4, 1/4, 1/2, 1/2, -1/2, -1/2,
                -1/2, 1/2, -1/2, 1/2, 0, 0, -1, 1), ncol = 4)
cont.mat <- solve(t(mat))
contrasts(dat$combo) <- cont.mat[,-1]

#Fit model
summary(lm(Rich ~ combo, data=dat))

Looking at the initial matrix, I've set the intercept to be the grand mean, coefficient 1 to be the difference between means of the sites, coefficient 2 to be the difference between the mean of Grasland from the mean of the other two combined, and the third to be the difference between the means of Grasland and Nat_bos within the Bruuk site only.

> mat
     [,1] [,2] [,3] [,4]
[1,] 0.25  0.5 -0.5    0
[2,] 0.25  0.5  0.5    0
[3,] 0.25 -0.5 -0.5   -1
[4,] 0.25 -0.5  0.5    1

Therefore, when we look at the summary, we see only the intercept and the first coefficient are significantly different from zero. This indicates that the mean richness from all sites is statistically not zero (usually an uninteresting finding) and that the difference between means of the sites is statistically not zero. For the other two contrasts, there isn't enough evidence to reject the null hypothesis (that the difference in means is zero).

> summary(lm(Rich ~ combo, data=dat))

lm(formula = Rich ~ combo, data = dat)

    Min      1Q  Median      3Q     Max 
-1.4018 -0.5291 -0.0279  0.5238  1.6052 

            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6.2037     0.1479  41.935  < 2e-16 ***
combo1        1.2447     0.2959   4.207 0.000272 ***
combo2        0.3561     0.2959   1.204 0.239606    
combo3       -0.2184     0.4184  -0.522 0.606122    
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7639 on 26 degrees of freedom
Multiple R-squared:  0.4282,    Adjusted R-squared:  0.3623 
F-statistic: 6.491 on 3 and 26 DF,  p-value: 0.001999

Careful selection of your contrasts will allow some precise testing of different hypotheses. One key thing to keep in mind when defining your own contrasts is to keep track of the signs and the sum of the contrasts. They influence the sign and magnitude of the coefficient.

  • $\begingroup$ Thank you! This is a new approach, good to know. I got the results, and I am left with two minor questions. I want to combine the levels for a broader site-to-site comparison, do I just average their means and then compare to the grand mean? And how can I combine the means of the 2nd and 4th level for this purpose, since the 4th does not have an estimate? $\endgroup$
    – B_rood
    Feb 1, 2017 at 13:26
  • $\begingroup$ The estimate of the 4th is just the deviance of its mean to the grand mean, I got that now. But what I would like to know most is if the combination of two levels of the same site (for a site-by-site comparison) significantly deviates from the grand mean. I cannot just combine p-values of individual levels, right? $\endgroup$
    – B_rood
    Feb 1, 2017 at 13:52
  • $\begingroup$ Correct, you can't just combine $p$-values. In the link above, #9 is user-defined contrasts. These are extremely powerful in defining comparisons between levels of a factor in a model. When I get a chance today, I'll edit my answer to include an example. $\endgroup$
    – Ashe
    Feb 1, 2017 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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