Solved – Mixed effects model or mixed design ANOVA in R

anovamixed modelrrepeated measures

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.

Best Answer

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.

Related Solutions

Solved – Does the blocking factor confound the fixed effect in the model

Here's what I would do. I expanded your dataset a bit so we have something to work with:

dat <- structure(list(year = c(2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 2001L, 
2001L, 2001L, 2001L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 2002L, 
2002L, 2002L, 2002L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 2003L, 
2003L, 2003L, 2003L), plot = c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 
3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 
9L, 9L, 9L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 
5L, 5L, 6L, 6L, 6L, 7L, 7L, 7L, 8L, 8L, 8L, 9L, 9L, 9L, 1L, 1L, 
1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 
7L, 7L, 7L, 8L, 8L, 8L, 9L, 9L, 9L), block = structure(c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L
), .Label = c("b1", "b2", "b3"), class = "factor"), treatment = structure(c(1L, 
1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 
2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 
2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 
3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 
1L, 3L, 3L, 3L, 2L, 2L, 2L, 1L, 1L, 1L, 3L, 3L, 3L, 2L, 2L, 2L
), .Label = c("C", "H", "L"), class = "factor"), tree = structure(c(1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 
3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L
), .Label = c("a", "b", "c"), class = "factor"), elevation = c(10.3, 
10.3, 10.3, 12.2, 12.2, 12.2, 14.4, 14.4, 14.4, 1.5, 1.5, 1.5, 
3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 5.5, 5.5, 5.5, 6.4, 6.4, 6.4, 6.9, 
6.9, 6.9, 10.3, 10.3, 10.3, 12.2, 12.2, 12.2, 14.4, 14.4, 14.4, 
1.5, 1.5, 1.5, 3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 5.5, 5.5, 5.5, 6.4, 
6.4, 6.4, 6.9, 6.9, 6.9, 10.3, 10.3, 10.3, 12.2, 12.2, 12.2, 
14.4, 14.4, 14.4, 1.5, 1.5, 1.5, 3.1, 3.1, 3.1, 2.4, 2.4, 2.4, 
5.5, 5.5, 5.5, 6.4, 6.4, 6.4, 6.9, 6.9, 6.9), growth = c(11L, 
7L, 11L, 4L, 5L, 12L, 7L, 8L, 11L, 6L, 10L, 10L, 5L, 10L, 8L, 
5L, 5L, 8L, 11L, 12L, 10L, 9L, 7L, 6L, 10L, 9L, 10L, 6L, 11L, 
6L, 8L, 5L, 6L, 11L, 9L, 12L, 9L, 8L, 8L, 12L, 5L, 8L, 4L, 12L, 
10L, 6L, 8L, 9L, 6L, 5L, 6L, 4L, 4L, 9L, 10L, 9L, 4L, 11L, 9L, 
11L, 5L, 10L, 11L, 4L, 10L, 4L, 10L, 6L, 10L, 4L, 9L, 10L, 9L, 
6L, 5L, 8L, 7L, 11L, 12L, 4L, 12L)), .Names = c("year", "plot", 
"block", "treatment", "tree", "elevation", "growth"), class = "data.frame", row.names = c(NA, 
-81L))

Now I am not super familiar with lme() so I am using lmer() from the lme4 package. Given your design, I would remove year as a fixed effect. It has 10 levels, and since it is a growth experiment, it will probably be significant because trees grow over time. Making sense of a model output with 3 main effects and their interactions when year has 10 levels might be a bit too much. That reduces your model to treatment and block (I wouldn't use elevation unless these small differences in elevation within the plots were purposely chosen).

As for the random part, I would include year (to account for multiple measurements on the same trees) and block:treatment, which represents the plots, because the 3 trees in the plots may not be spatially independent(?!), which should be accounted for.

library(lme4)
m1 <- lmer(growth ~ treatment * block + (1|year/block:treatment), data = dat)
summary(m1)

Linear mixed model fit by REML ['lmerMod']
Formula: growth ~ treatment * block + (1 | year/block:treatment)
   Data: dat

REML criterion at convergence: 364

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.6991 -0.8580  0.2137  0.6845  1.6875 

Random effects:
 Groups               Name        Variance Std.Dev.
 block:treatment:year (Intercept) 0.4115   0.6415  
 year                 (Intercept) 0.0000   0.0000  
 Residual                         6.6914   2.5868  
Number of obs: 81, groups:  block:treatment:year, 27; year, 3

Fixed effects:
                   Estimate Std. Error t value
(Intercept)          8.3333     0.9384   8.880
treatmentH           1.0000     1.3271   0.753
treatmentL          -0.4444     1.3271  -0.335
blockb2             -0.6667     1.3271  -0.502
blockb3              0.1111     1.3271   0.084
treatmentH:blockb2  -1.2222     1.8769  -0.651
treatmentL:blockb2   1.0000     1.8769   0.533
treatmentH:blockb3  -1.2222     1.8769  -0.651
treatmentL:blockb3  -0.7778     1.8769  -0.414

To get p-values for the fixed effects you could use the mixed() function from the afex package:

library(afex)
mixed(growth~treatment*block+(1|year/block:treatment), data=dat)

Contrasts set to contr.sum for the following variables: treatment, block
Fitting 4 (g)lmer() models:
[....]
Obtaining 3 p-values:
[...]
           Effect    df F.scaling    F p.value
1       treatment 2, 16      1.00 0.27     .76
2           block 2, 16      1.00 0.51     .61
3 treatment:block 4, 16      1.00 0.51     .73

However, note that this may not be the only way to analyze your dataset. Depending on the actual questions you have, you can build multiple models and compare them. To generally get a better idea about fixed and random effects, you could start reading here.

Best Answer

Related Solutions

Solved – Does the blocking factor confound the fixed effect in the model

Related Question