# Mixed effects model or mixed design ANOVA in R

I'm wondering if you could quickly advise on the best statistical analysis for my data.

I have a designed experiment with 12 plots within one larger area. This area was blocked into 6 to overcome issues of natural gradients present there, so there are 2 plots in each block.

Treatment 1 is a two level treatment split between these 12 plots i.e. 1 plot in every block receives one level of the treatment ("High") and the other receives another level ("Low"). Nested within this, each plot is split into two (24 subplots), and given a second two level Treatment 2 (levels "A" and "B"). Every single plot receiving Treatment 1 has a Treatment Two Subplot receiving level "A", and a subplot receiving level "B".

As such, I have a two-way factorial design, with 4 types of treatment combination: High/A, High/B, Low/A and Low/B.

This experiment was run for three years, and species richness measured as a response (a simple count of the number of species present in each Treatment 2 subplot) each year.

I am interested to see what the effects of the two levels of Treatment Two are on Richness change with time, in the context of the level of Treatment 1 received.

I am not sure if I need:

• A two-way mixed ANOVA, with within subjects effects (Year?) and between subjects effects (Treatment 1 and Treatment 2?)
• A mixed effect model, with Treatment 1 and Treatment 2 as fixed effects, and Year and Plot as random effects. (Is Year a random effect? I am unsure...). Should subplot also be a random effect here?
• A mixed effect model with a time series.
• I am also unsure as to what to do with block. Block doesn't have levels of treatment, it was just an attempt to reduce the incidence of similarity between plots as a result of the natural gradients in the study area...
• In addition, because the richness data is a count of species, i'm not sure if this affects my choice of analysis...

If you have any thoughts on any of these points, I would be glad to hear them.

• I have trouble understanding what you mean by 'plots' and 'blocks'. Do you want to treat these/Do you think these can be treated as categorical variables in your analysis? – robert Jul 17 '13 at 9:25
• Hi Robert, thanks for your comment. The plots are the units in my experiment, akin to say subjects in a quetsionnaire. There are 12 plots in one area. This area was split in 3 "blocks", with 4 plots allocated in each block. Treatment 1 was randomly assigned to plots in each block in turn, to ensure High/Low treatments were seen in plots in every block i.e. in every section of the area. This is an attempt to reduce the non-independence of plots due to similarities that occur between them as a result of the natural gradients in the area (e.g. water gradients, light gradients etc). – Sarah Jul 17 '13 at 9:46
• Ecah plot was split into two (creating two subplots per larger plot) and Treatment Two randomly assigned, such that one subplot per plot received level "A" of Treatment 2 and one subplot received level "B". – Sarah Jul 17 '13 at 9:47

So I've done a lot of reading and chatting to people and I have a solution.

My experimental design is a split plot design, which is quite different from a nested or hierarchical design. I was originally confusing the terms. As Robert correctly states in his answer, what is needed is a mixed effects model. Thus:

Fixed effects: Year, Treatment1, Treatment2

Random effects: Year, Block, Treatment1

The model is specified thus:

mod<- lmer(Richness~Treatment1*Treatment2*Year+(1|Block/Treatment1)+(1|Year),data=dat,poisson)


The fixed effects are the terms specified in the brackets. Since none of these are continuous (the effect of Year doesn't necessarily increase each year in a linear fashion so I have classed it as a categorical fixed effect), they are specified 1|fixed effect, where 1 represents the intercept.

If Block were actually a continuous fixed effect (obviously hypothetical!) then the fixed effects might be specified +(Block|Treatment1)+(1|Year).

The model can then be simplified as appropriate.

Several things to note:

1) When specified as a random effect, Year is listed separately from Block and Treatment1, since it doesn't have an intuitive "level" at which to be nested between them (Year isn't any different at any plot size of the experiment: for every block, plot and subplot Year is the same.

2) Treatment 2 does not need to be specified as a random effect since it represents the highest level of replication in the experiment and therefore will not be psuedoreplicated

3) In mixed effects models it is possible to specify an error distribution if errors are not normal. I have specified poisson here, since my response data are counts - this improved the distribution of the model residuals.

My suggestion is to use a mixed effects model with

• Treatment 1 and 2 as factors (fixed effects)
• Year as an interval-scaled factor
• Block as a random factor

Q: Should year be treated as a random effect?

A: Depends on what you expect based on your theory. I suppose the number of species might increase over time, and this increase might be linear. Then you should use year as a factor. However, if there is reason to assume that, say, in year one there is little increase, in year two a great increase, and again little increase in year three, you might consider using year as a fixed factor.

Q: How should 'block' be treated?

A: My understanding is that environmental conditions differ between blocks, but not between plots in the same block. Then I would suggest adding blockas a random factor.

Does that help in answering your questions?

• Hi Robert, when you say interval-scaled factor, do you mean an ordered variable? Would I create this by simply saying as.factor(Year)? And then ordering this? – Sarah Jul 17 '13 at 11:39
• And what is the difference between Year as a factor and Year as a fixed factor in this context? – Sarah Jul 17 '13 at 11:40
• By interval-scaled I mean year should be a number. as.factor() in R would make it a categorical variable/factor. If your dataframe in R "knows" that year is numerical variable all is fine, if not you could use as.numeric(as.character(year)) to make sure it is. – robert Jul 17 '13 at 11:46
• And regarding the difference between factor and fixed factor, I would use these as synonyms here. But there's a chance I'm wrong about the terminology here. – robert Jul 17 '13 at 11:49