Solved – Random effect nested under fixed effect model in R

rrandom-effects-model

Looking for help to create a model for my data gathered from a reciprocal transplant experiment: I have 2 populations of fish (A and B) and 2 temperature regimes (warm and cold) that were crossed with each other (= 4 treatment groups- A:warm, A:cold, B:warm, B:cold). Each treatment group consists of two replicate tanks (8 tanks total). There are 100 fish in each tank. The response variable of interest is growth rate and I have recorded growth rate for each of 800 fish.I want to test the effects of population and temperature and their interaction. I also want to include tank effect. I think I should nest tank effect within the interaction effect (Pop:Temp which is equivalent to the treatment group) and include nest effect as a random effect. So I should have Population (fixed effect) and Temp (fixed effect) and tank effect (random effect) nested within the interaction of Pop:Temp. Is it possible to create a random effect nested within a fixed effect interaction? If so, how do I do this in R?

Best Answer

It doesn't make sense to both include tank as a random effect and nest tank within the pop/temp fixed effect. You only need one of these, depending on how tank is coded.

If tank is coded 1-8, you only need the tank random effect. Nesting it within the pop/temp fixed effect results in the same 8 units, so is not necessary.

If tank is coded 1-2 (that is, which rep it was), you only need to nest tank within the pop/temp fixed effect, because that gives you your 8 unique tanks. Including the tank random effect is only desired if the tanks were first divided into two groups and then randomized to treatment; if the eight tanks were completely randomized to treatment, this is not necessary.

You could do this with likelihood based solutions such those in nlme and lme4 but if everything is balanced, it might be simpler to use the traditional ANOVA approach using aov.

Creating some sample data:

set.seed(5)
d <- within(expand.grid(pop=factor(c("A","B")),
                        temp=factor(c("warm", "cold")),
                        rep=1:2,
                        fish=1:100), {
                          tank <- factor(paste(pop, temp, rep, sep="."))
                          tanke <- round(rnorm(nlevels(tank))[unclass(tank)],1)
                          e <- round(rnorm(length(pop)),1)
                          m <- 10 + 2*as.numeric(pop)*as.numeric(temp)
                          growth <- m + tanke + e
                        })

Using aov like this:

a0 <- aov(growth ~ pop*temp + Error(tank), data=d)
summary(a0)

or lme like this:

library(nlme)
m1 <- lme(growth ~ pop*temp, random=~1|tank, data=d)
anova(m1)
Related Question