#Factorial MANOVA.  Requires candisc and car or Rcmdr libraries.

manova3 <- read.table("http://www.unt.edu:8080/rss/class/mike/data/manova3.dat",header=TRUE)
manova3$Achieve=factor(manova3$Achieve,1:2, labels=c("High","Low")) #make the numeric values factor levels with labels
manova3$Method=factor(manova3$Method,1:4)
attach(manova3)

Y=cbind(Test1,Test2)  #bind the DVs
manovamod=lm(Y~Achieve + Method +Achieve:Method, data=manova3)
Manova(manovamod)

library(candisc)
#Achieve main effect
(term1=candisc(manovamod, term= "Achieve"))  #only 1 canonical variate since k-1 groups equals 1
term1$structure
plot(term1)  #still new to these functions, I'm not exactly sure what the structure part is showing as the loadings are very similar but the plot would suggest a huge difference (it makes more sense in with 2 variates)
#Method main effect
(term2=candisc(manovamod, term= "Method"))  #The method main effect is large but mostly along the first composite
term2$structure #both strongly load on the first function
plot(term2)   #The plus signs are the mean canonical score for each group (for each variate).  The circles are created to encompass 95% CIs for the groups canonical means.  
		  #The arrows are the projections of the original variables into the canonical 'space'.  They start at the origin 0,0 representing the mean on the two canonical variates.  
              #Note Test1 loadings then notice how the Test 1 projection moves into positive space for the first dimension but is associated with negative values for the 2nd.
		  #Test 2 is positive for both, but does not extend far (i.e. does not share much variance with) the 2nd dimension
#Interaction
(term3=candisc(manovamod, term= "Achieve:Method"))  #The interaction appears notable along two dimensions
term3$structure  #The first dimension seems a composite of both, while the second one mostly represents the science test
plot(term3)  #The circles still represent each 'group' which in this case is each of the 8 cells of the 2x4 interaction

heplot(manovamod)  #all effects together
#heplot(manovamod, size="eff")  #alternative

#Univariate interaction plots
interaction.plot(Method, Achieve,Test1)
boxplot(Test1~Method:Achieve)
interaction.plot(Method, Achieve,Test2)


#Example of simple effects on the first canonical variate for the interaction
achievemod=lm(Can1~Achieve+Method+Achieve*Method,data=term3$scores)
library(contrast)
contrast(achievemod,a=list(Achieve="Low", Method=levels(Method)), b=list(Achieve="High", Method=levels(Method)))  #Suggests that most of the difference in achievers regards method 4 (maybe the second as well)

#Second canonical variate
achievemod2=lm(Can2~Achieve+Method+Achieve*Method,data=term3$scores)
contrast(achievemod2,a=list(Achieve="Low", Method=levels(Method)), b=list(Achieve="High", Method=levels(Method)))  #Here the difference lies mostly with the first two methods