#Example of a one-way MANOVA.  Packages needed: candisc (heplots comes with it), and Rcmdr or car for the Manova function

#Basic Manova
defattract <- read.table("http://www.unt.edu/rss/class/mike/data/juryattract.dat",header=TRUE)
defattract$pa = factor(defattract$pa,levels=c(1, 2, 3), labels=c("Attractive", "Average", "Unattractive"))
attach(defattract)
Y = cbind(warm, kind, sophist)
manova.mod <- manova(Y ~ pa, data=defattract)
Manova(manova.mod)  #Add test.statistic="Wilks" for example.  Also available are "Hotelling-Lawley" and "Roy"

#DV and Post hocs
library(candisc)
modb=candisc(manova.mod)  #Test in terms of the canonical correlation
modb
(defattract2=modb$scores)  #If you want to do post hocs on the actual function of the DVs, here are the scores.
modb$structure #DV loadings on the composites
posthocmod1=aov(Can1~pa, data=defattract2)
TukeyHSD(posthocmod1)


#Graphics for better post hoc
plot(modb)  #circles represent 95% CIs for the canonical group means.  In the history menu check 'recording' to save graph.
heplot(manova.mod, verbose=T)  	#shows directly how the group means differ on two DVs jointly, in relation to the within-group variation; a "significant" hypothesis ellipse extends outside of the error ellipse
#heplot(manova.mod, size="effect.size")	#puts the hypothesis ellipse on the same scale as the error ellipse
heplot(modb)  #I find the plots related to the canonical variates more interpretable
#heplot(modb, size="effect.size")  
heplot3d(manova.mod) #maximize the window and use the mouse wheel for zoom
pairs(manova.mod)  #for some reason this may spit out warnings, but the graphs are fine
#pairs(manova.mod, size="effect.size")