Biology, images, analysis, design...
Use/Abuse Principles How To Related
"It has long been an axiom of mine that the little things are infinitely the most important" (Sherlock Holmes)

Search this site



ANOVA for blocked designs

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)
HabitatSite (=block)
Long Island
Eelgrass (eel)
Bare (bar)
Longlines (lon)
Picked (pik)
Dredged (drd)
Hummocks (hum)

  1. 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.

  2. 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 Guido Wyseure. Unfortunately they give quite different results! We have used the latter version since it gives the same results as those given by Underwood (1997)

    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.

  3. Get table of means

  4. Carry out analysis of variance

    Sums of squares can be calculated manually as follows:

    SStotal  =  16856.69  −  (546.1)2  =  288.62

    SSlocations  =  1812  +   1812  +   184.12  −  546.12
       =  16569.14 - 16568.07 = 1.07

    SShabitats  =  97.32  +  82.92  +  102.32  +  94.32  +  96.12  +  73.22  −  (546.1)2  
       =  16763.64 - 16568.07   = 195.57

    SSresidual  =  288.62 - 195.57 - 1.07  =   91.98

    These values are inserted into the ANOVA table and the mean squares calculated.

    ANOVA table
    Source of
    df SS MS F-
    Locations 2 1.07 0.534 (0.058) (0.944)
    Habitats  5  195.57 39.12 4.25 0.025
    Residual 10 91.98 9.198    
    Total 17 288.62      

    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).
  5. 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 habitat.

    We get the same result if we use the glht() function in the multcomp package:


    So what can we conclude??

    1. The evidence for any consistent difference in microbial activity between habitat types is weak.
    2. 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).
    3. 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).