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.