Worked example 1
Our first worked example uses interpolated data
from Richardson et al. (2008) on the effect of different oyster aquaculture techniques on
microbial aerobic activity at Willapa Bay. The author described the design as a
randomized complete block with three locations (= blocks) and six habitat types (fixed
factor) within each block. However, since it is an observational study, treatments
(habitat types) were clearly not randomly allocated, and it would be better described as an
unreplicated two factor design. It can be analysed using a randomized block ANOVA but
some aspects are problematical.
| Aerobic sole-source carbon use (SCU)
Draw boxplots and assess normality
Plot out data to get a visual assessment of the treatment and block effects, and assess how
appropriate (parametric) ANOVA is for the set of data.
First impressions are that the block effect seems fairly small but there may well be a
difference between habitats. There is nothing too horrendous in terms of distributions - at
least nothing that a transformation would improve - so we proceed to an interaction plot.
Check for location × habitat type interaction
We start with a simple location × habitat type interaction plot. If there were no
interaction between locations and habitat types, microbial activity (SCU) should follow
parallel trends across habitat types for each of the different locations.
Microbial activity seems to follow similar trends across habitats in Nemah and Long
Island, but not in Nahcotta. Hummocks have the highest SCU readings in Nahcotta, but
one of the lowest in Nemah. We cannot make a general test on whether this interaction is
significant because of the lack of replication. However, we can test for whether there is
significant interaction of the multiplicative form using the Tukey non-additivity test.
There appear to be two versions of this test available in R - the version in the package
alr3 which is tukey.nonadd.test() and a version provided by
Unfortunately they give quite different results! We have used the latter version since it
gives the same results as those given by Underwood
In this case P = 0.182. Using the conventional P < 0.05 level, we
would say there is no evidence for a multiplicative interaction. However, using more
liberal significance level of P = 0.25 (as widely recommended), there are clearly
grounds for suspecting interaction. For now we will continue with the analysis, but will
return to the point later.
Get table of means
Carry out analysis of variance
Sums of squares can be calculated manually as follows:
||16569.14 - 16568.07|| = 1.07
|| 16763.64 - 16568.07
|| = 195.57
||288.62 - 195.57 - 1.07 = 91.98
These values are inserted into the ANOVA table and the mean squares
The habitat factor came out significant at P = 0.025.
Assuming blocks are fixed, SEhabitat mean = √(9.198/3) = 1.751
Assuming blocks are random, SEhabitat mean = √[(9.198+0.534)/3)] = 1.801 (Note taking this into account
makes little difference in this case because the block effect is very slight)
SEdifference between habitat means = √
[(2 × 9.198)/3)] = 2.476
Testing between blocks is inappropriate if blocks is a random factor.
In R there are several different ways to do this analysis.
- Use lm() with the model defined as model=lm(resp~B+A).
- Use aov() with the model defined as model=aov(resp~B+A)
- Use aov() with the model defined as model=aov(resp~A, error=B).
Check model diagnostics
Since there is only one error term in this analysis, we can use the standard residual
plots produced by R to examine diagnostics.
The residuals versus fitted plot shows no clear trend, and the residuals are reasonably
normally distributed. If we had random allocation of treatments within blocks, we could
proceed to comparing mean. But in this case we are dealing with an observational study
with no random allocation of treatments. Hence it would make sense to check the
homogeneity of covariances assumption by estimating sphericity
and (just for good measure) the variance-covariance matrix.
Unfortunately epsilon comes out at only 0.19 indicating a major deviation from
sphericity. This casts serious doubt on the validity of the analysis. We attempted to fit
the linear mixed effect model (using lme()) but the model with unstructured covariance
structure would not converge, and none of the other covariance structures were
significantly better than the compound symmetry model which gave the same P-
value for the treatment factor of 0.0247.
Hence we proceed to test differences between habitats. In this case we will use
Tukey's HSD test. One way to do it is to obtain a list of all differences between means
and compare these differences with the honestly significant difference:
Only differences greater than 8.6 are significant at P = 0.05. Hence the only
significant difference in microbial activity is between the eelgrass habitat and the picked
We get the same result if we use the glht() function in the multcomp package:
So what can we conclude??
- The evidence for any consistent difference in microbial activity between habitat
types is weak.
- Eelgrass sites appear to have a higher microbial aerobic activity than picked sites
(P = 0.0254). There is also a suggestion that longlines and bare sites have a higher
microbial aerobic activity than picked sites (P < 0.10).
- But there are serious threats to validity of these conclusions, namely the possibility of
a habitat × location interaction (meaning that habitat differences are not consistent
across locations) and the unreliability of the P-values (given that homogeneity of
covariance assumptions are not met).