Should a mixed effects model be used?

Question

This post provides an excellent example of the inner workings of a mixed effects model: http://emhart.github.com/blog/2012/11/16/making-sense-of-random-effects/

In a hypothetical study, I have measured the wing length of a bird species at 10 different locations:

df <- data.frame(wing.length=rnorm(30),
                 location=paste('location', rep(1:10, 3)), 
                 county=paste('county', sort(rep(1:3, 10))))

First 15 rows of dataframe:

   wing.length    location   county
1     29.29024  location 1 county 1
2     30.06387  location 2 county 1
3     27.72127  location 3 county 1
4     29.74502  location 4 county 1
5     29.85506  location 5 county 1
6     30.26669  location 6 county 1
7     30.58748  location 7 county 1
8     30.47608  location 8 county 1
9     31.86882  location 9 county 1
10    28.87578 location 10 county 1
11    30.00726  location 1 county 2
12    30.52488  location 2 county 2
13    30.64339  location 3 county 2
14    29.93695  location 4 county 2
15    27.86217  location 5 county 2

However, unlike the mixed effects model in the linked blog post where each individual is measured 5 times, each location is only measured once.

In my hypothetical example, I want to test the hypothesis that wing length varies among counties. County is a fixed effect and location is a random effect. Should a type of mixed effects model be applied, or would a two-way ANOVA be more appropriate? If a mixed effects model is needed, then what type?

Answer 1

To test a random effects model of a continuous DV using nominal, nested predictors in R, Lowe (2013; A.K.A. @conjugateprior) recommends the lmer function in the lme4 package (Bates, Maechler, Bolker, & Walker, 2013). For your particular data frame and hypothesis, here's the code:

lmer(wing.length~location+(1|location:county),data=df)

With wing.length=rnorm(30), location:county has only one observation per group. Therefore this produces the following warning message (and no other output):

Error in checkNlevels(reTrms$flist, n = n, control) : 
  number of levels of each grouping factor must be < number of observations

However, it works with wing.length=rnorm(150) (location and county just repeat to fill in the rows). If you wrap the code in summary(), you'll get a $t$ value for your intercept and fixed effects. The $t$ value for the intercept achieves significance when I add rep(1:3,50) to wing.length, so all seems well. :)

References

<sup>- Bates, D., Maechler, M., Bolker, B., & Walker, S. (2013). lme4: Linear mixed-effects models using Eigen and S4. R package version 1.0-4. http://CRAN.R-project.org/package=lme4.
- Lowe, (2013, January 10). Formulae in R: ANOVA and other models, mixed and fixed. Conjugate Prior. Retrieved from http://conjugateprior.org/2013/01/formulae-in-r-anova/.
- R Core Team (2012). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/.</sup>

Answer 2

I guess that you did not mean that each location belongs to the three counties, so location in your database should be (1:30) and not rep(1:10, 3). If true, then you only have one measurement per location. Then, your design as given only has one factor of interest (fixed factor county) and it can be solved with a one-way three-level ANOVA. Though you actually measured lenghts from wings from birds, your samples here are the locations ("random", as samples are), of which you happened to measure ten in each county.

However, I will follow as if you really meant to have several measured birds in each location. Then, if I understand correctly, you consider that wing lenghts from birds captured in the same location are related by some (unknown) factor affecting their wing size (e.g., same parents, same resources). And you are not interested in that but on evaluating differences in average wing length among (only those) three counties (your fixed factor). So birds from the same location are considered not independent samples (of wing length) or pseudo-replicates (for that hypothesis of yours). If above is OK, then location is a random factor nested in county. I call that a nested design, with some samples per location and ten locations per county. A proper analysis is that which considers that you only have ten independent locations to evaluate county differences and not N independent samples of wing length.

For the analysis I would use:

With lme4 package: lmer(wing.length ~ county + (1|location), data=df)

With nlme package: lme(fixed = wing.length ~ county, random = ~1 | location, data = df)

or, if under assumptions (balanced, normal errors, you know...), you can simply go to a nested ANOVA: aov(wing.length ~ county + Error(location), data=df)