How to perform a planned comparison test (in R, using gmodels)
Task
Suppose that ANOVA reveals a significant difference between treatment levels, and you wish to explore further through post hoc analysis by comparing two specific treatment levels. How can we perform perform planned comparisons, also called a contrast test?
Solution
Usually, you have data you wish to compare, but we will use example data here.
We load the “oats” dataset from R’s MASS package, about the yield of oats
from a split-plot field trial using three varieties (V) and four levels of
manurial treatment (N). The experiment was laid out in 6 blocks (B) of 3 main
plots, each split into 4 sub-plots. The varieties were applied to the main plots
and the manurial treatments to the sub-plots.
1
2
3
# install.package('MASS') # if you have not already done so, and want this data
library(MASS)
df <- oats
Before we perform the contrast test, let’s verify that the yield of oats Y
depends on the nitrogen manurial treatment given to it N.
1
2
aov1 <- aov(Y ~ N, data = df)
summary(aov1)
1
2
3
4
5
Df Sum Sq Mean Sq F value Pr(>F)
N 3 20020 6673 14.2 2.78e-07 ***
Residuals 68 31965 470
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
The $p$-value in the Pr(>F) column is below $\alpha=0.05$.
So at the 5% significance level, the yield differs according to the
nitrogen manurial treatment. We assume that the model assumptions are met but
do not verify them here.
We now want to perform a planned comparison test (or contrast test) on the data
to see whether there is a difference between the $N<0.5$ levels and the
$N>0.5$ levels. We will use the fit.contrast function in the gmodels package.
Since the order of the levels is 0, 0.2, 0.4 and 0.6, the
contrast coefficients will be $-0.5$, $-0.5$, $0.5$, $0.5$, respectively.
1
2
3
# install.package('gmodels') # if you have not already done so
library(gmodels)
fit.contrast(aov1, "N", coeff=c(-1/2,-1/2,1/2,1/2))
1
2
3
4
Estimate Std. Error t value Pr(>|t|)
N c=( -0.5 -0.5 0.5 0.5 ) 29.66667 5.110338 5.805225 1.855598e-07
attr(,"class")
[1] "fit_contrast"
The $p$-value in the Pr(>|t|) column is below $\alpha=0.05$.
This tells us that there is a significant difference between the average yields
of the $N<0.5$ and $N>0.5$ levels.
Content last modified on 24 July 2023.
See a problem? Tell us or edit the source.
Contributed by:
- Krtin Juneja (KJUNEJA@falcon.bentley.edu)
- Nathan Carter (ncarter@bentley.edu)