Multivariate
ANOVA (MANOVA)
We'll use a simple example adapted from 'Discovering
Statistics with R' by Andy Fields (by the way – it's a pretty good book
explaining statistics).
First, download the data file into your workspace
from OCD and load
it:
> ocdData=read.delim("OCD.dat",
header = TRUE)
If you want, have a look at the data (note one
categorical IV and two numerical DVs)
> ocdData
You should also run some descriptive stats.
Let's code Group levels with shorter names:
> ocdData$Group=factor(ocdData$Group, levels = c("CBT", "BT",
"No Treatment Control"), labels = c("CBT", "BT",
"NT"))
The formula is the same as for ANOVA (i.e. DV ~ IV)
but this time DV is not a vector but a matrix where each column is a separate
DV.
First, we need to create this matrix of DVs:
> outcome=cbind(ocdData$Actions, ocdData$Thoughts)
Then we use manova() function to generate the model:
> ocdModel=manova(outcome ~ Group, data = ocdData)
Let's look at the output of the model and the test
statistic:
> ocdModel
> summary(ocdModel,
intercept = TRUE)
Next, you can inspect the differences between
various tests which can be applied to the model:
> summary(ocdModel,
intercept = TRUE, test = "Wilks")
> summary(ocdModel,
intercept = TRUE, test = "Hotelling")
> summary(ocdModel,
intercept = TRUE, test = "Roy")
Watch the difference in p-values!
We can also test for effects on two DVs separately:
> summary.aov(ocdModel)
Alternatively, the same results can be achieved by
explicitly testing two one-way ANOVA models:
> actionModel=aov(Actions ~ Group, data = ocdData)
> thoughtsModel=aov(Thoughts ~ Group, data = ocdData)
> summary(actionModel)
> summary(thoughtsModel)
In our case those two effects are not significant,
so you cannot infer anything from planned comparisons. Otherwise, the
procedure for the contrasts is the same as for one-way ANOVA (since actually
you've just run two separate one-way ANOVAs)
> summary.lm(actionModel)
> summary.lm(thoughtsModel)
Finally, there are non-parametric alternatives to
MANOVA mulrank() and cmanova().
They are both part of WRS package, if you have spare time try installing
it (instructions are in Tuesday session) and making it work with this dataset.