Random effect nested under fixed effect model in R

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?

• How can there be an effect of POP (or interaction) of any kind if there are always 100 fish in each tank?
– John
Commented Jul 13, 2011 at 20:49
• I think you could think about this particular fish population as a "species", or a "group" (even though they may be the same species, but from a different region, or a different breed). So in effect, I think, katy d. might be looking for an effect of this "group" and temperature on growth rate. Commented Aug 15, 2011 at 7:30

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)

• Nice answer! Could you explain the reasoning behind this: "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."? Commented Feb 6, 2017 at 14:25

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.

• +1 Nice job thinking about (and answering) the OP's underlying question. Commented Aug 16, 2011 at 1:24

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

• I'm not very good at the whole nested interpretation that the OP wants... but to test my limited understanding of this lmer model it is: The intercept is the sum of the fixed common intercept plus a random, per-tank intercept. The slope is the sum of a fixed treatment slope and a random, treatment-per-tank slope. The per-tank intercept and slope are not correlated. Is that correct? Commented Aug 15, 2011 at 17:17
• yep, that was it Commented Aug 25, 2011 at 21:59