Random effect nested under fixed effect model in R

Question

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?

Answer 1

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)

Answer 2

Because the treatment (pop*temp) does not vary within levels of the random effect (tank), this is a simple nested design (I think). Unless you are particularly interested in the amount of between-tank variation, it will be much easier just to aggregate the data to the level of tank and then run a simple (non-mixed) linear model -- the statistical inferences should be identical to the tests of the fixed effects in the mixed ANOVAs above.

Using the sample data above (I'm sure there's an easier way to do the aggregation):

d2 <- rename(cast(subset(d,select=c(pop,tank,temp,growth)),
                  pop+temp+tank~.,fun.agg=mean,value="growth"),
             c("(all)"="meangrowth"))

m1 <- lm(meangrowth~pop*temp,data=d2)
summary(m1)

See Murtaugh (2007) "Simplicity and complexity in ecological data analysis", Ecology, 88(1), 2007, pp. 56–62.

Answer 3

Did you try: lmer(growth~1+treatment+(1|tank)+(0+treatment|tank))