InfluentialPoints.com
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)

Nested ANOVA

Worked example 1

Our first worked example uses data from Alexander Kerr (Guam Marine Lab) on a (possibly hypothetical) study on whether site exposure to waves affects coral size. Two transects were randomly located at each of four sites, which varied in their degree of wave exposure. Hence sites is a fixed factor, chosen to represent a certain treatment magnitude, and transects are a random factor nested within sites. Diameters of mound-forming faviid corals found along the transect were measured.

 Coral diameters in different sites Site: Urunao Tarague Hila'an Ritidian Transect: 1 2 3 4 5 6 7 8 3128332432 3633382937 585216 8118519 2927303324 3331343728 3846264233 4250304637

Note that if exposure could be assessed on a measurement scale, then the group (treatment) factor could be assessed directly using a regression-type design. As it is with the available information, and using ANOVA, we can only compare sites (rather than wave exposure).

1. Draw boxplots and assess normality

Plot out data to assess how appropriate (parametric) ANOVA is for each level of the data set. The transect means are the replicates for comparing sites so we need an aggregated dataset on which to base the boxplots. Unfortunately in this example we only have (b = 2) transects in each site so the boxplots are uninformative - neverthless do it below for the sake of completeness.

We then assess the distribution of observations within each transect.

The boxplots suggest that variances are fairly similar for each transect. We assess normality of observations using normal QQ plots for each subgroup.

These data seem remarkably 'well behaved' although we have single unusually high values in each of transects 3 and 4. At least we can be reasonably certain that the site means follow a normal distribution given central limit theorem.

 Some authors (notably Logan (2010)) suggest that the distribution of observations within subgroups is unimportant given that the comparison of fixed factor levels is based on subgroup means. Whilst there is an element of truth in this, it is still extremely important to assess within subgroup distributions for two reasons: In many cases (as with this one) the number of observations of subgroup means will be too small to adequately assess the distribution of means. If you also have evidence of normality within subgroups, then one has more grounds for believing that the subgroup means are normally distributed. Just relying on central limit theorem with small sample sizes is very unwise! If you are a 'pooler' (we are not, but have to accept that some people have still not seen the light), then you will use the P-value for subgroups within groups to decide whether or not to pool. But the validity of this test depends critically on the distribution of observations within subgroups.

2. Check homogeneity of variances

 Using R Bartlett test of homogeneity of variances data: siz by trn Bartlett's K-squared = 6.8008, df = 7, p-value = 0.4499

Since P = 0.450, there is no evidence for heterogeneity of variances.

3. Carry out analysis of variance

Sums of squares can be calculated manually as follows:

 SStotal = 36714 − (1104)2 = 6243.6 40

 SSsites = 3212 + 872 + 3062 + 3902 − (1104)2 10 10 10 10 40 = 35634.6 - 30470.4 = 5164.2

 SStransects within sites = 1482 + ..... + 2052 − 35634.6 5 5 = 35799.6 − 35634.6 = 165

 SSwithin transects = 6243.6 − 5164.2 − 165  =   914.4

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

 ANOVA table Source ofvariation df SS MS F- ratio P Sites 3 5164.2 1721.40 41.73 0.002 Transects w'in sites 4 165.0 41.25 1.44 0.2435 Within transects 32 914.4 28.57 Total 39 6243.6

The difference between sites was highly significant (P = 0.002). The difference between transects within sites was not significant (P = 0.244).

There are two ways to carry out a nested ANOVA in R. The most straightforward approach is to use the aov function with multiple error strata.

 Using RError: tran Df Sum Sq Mean Sq F value Pr(>F) site 3 5164.2 1721.40 41.731 0.001779 ** Residuals 4 165.0 41.25 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 32 914.4 28.575

This gives the P-value for sites as 0.0018 correctly using the MStransects w'in sites as the denominator for the F-ratio. The P-value for transects within sites is not given but can be readily calculated as MSresiduals:tran / MSresiduals:within.

Use of aov() is only appropriate for balanced designs and a lme () (a linear mixed-effects model) should be used for unbalanced designs. For a balanced design it will give the same mean square and F-ratio as the aov function.

 Using R numDF denDF F-value p-value (Intercept) 1 32 738.6758 <.0001 site 3 4 41.7309 0.0018

Both methods give an F-ratio between sites of 41.73 between sites with a P-value of 0.0018 as we obtained in the manual calculations.

4. Should we pool denominators?

Reference to the full ANOVA table above shows that transects within sites were not a significant source of variation. The provider of these data used this fact to justify pooling MStransects within sites and MSwithin transects to provide a new estimate of the Mean square errorwithin. We would regard this as pseudoreplication, and would not pool. In fact, since the P- value is < 0.25, even well known advocates of pooling (for example Quinn & Keough (2002) and Doncaster & Davey (2007)) would not pool in this case. The key point is that using a mean square error with 36 degrees of freedom implies that you have taken 10 random independent samples from each type of site - which you simply haven't! You have only sampled from two randomly located transects in each type of site - your analysis should reflect this fact!

5. Examine model diagnostics

The next step is to examine the fitted model diagnostics focusing on the first error stratum. The figure below shows residuals versus fitted values, and a qq plot to assess normality:

There is no evidence of any increase in the spread of residuals as fitted values increase so we accept that the model is adequate.

6. Comparison of means

Lastly we compare means. Kerr only made one comparison (site A versus Site B) but we will assume a planned full set of orthogonal comparisons, namely site A versus Site B, site C versus site D, and mean of site A & B versus mean of site C+D (whether this is a logical set of comparisons depends entirely on the particular study - we do it here only as an example). We do not give the manual calculations for this but only the R-code and output:

 Using Rcrossprod(contrasts(ecoker\$site)) [,1] [,2] [,3] [1,] 2 0 0 [2,] 0 2 0 [3,] 0 0 4 Error: tran Df Sum Sq Mean Sq F value Pr(>F) site 3 5164.2 1721.40 41.7309 0.001779 site: siteA vs siteB 1 2737.8 2737.80 66.3709 0.001235 site: siteC vs siteD 1 352.8 352.80 8.5527 0.043049 site: Mean site A & B vs Mean site C & D 1 2073.6 2073.60 50.2691 0.002089 Residuals 4 165.0 41.25 site ** site: siteA vs siteB ** site: siteC vs siteD * site: Mean site A & B vs Mean site C & D ** Residuals --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Error: Within Df Sum Sq Mean Sq F value Pr(>F) Residuals 32 914.4 28.575

In the first table the values above or below the cross-product matrix diagonal are all zero. This confirms that the set of planned contrasts is orthogonal. The next table decomposes the SSsites into the three planned contrasts - all prove to be significant.