setwd("\\\\sttnas/FacultyReds/pszhong/My Documents/MSU-Teaching/STT-843-2018-Spring") ## Please provide the working directory here.
load(file="Alldata.Rdata")

############################################################
## Obtain the response and design matrix
## Trtment 1: alcohol group; Trtment 2: grape juice;
## Trtment 3: red wine; Trtment 4: water.
############################################################

gene1to5<-sapply(Alldata$originaldata[c(1:5),c(3:110)],as.numeric)
sampleIDs<-colnames(gene1to5)
trtment<-(sampleIDs%in%Alldata$trt1)*1+(sampleIDs%in%Alldata$trt2)*2+(sampleIDs%in%Alldata$trt3)*3+(sampleIDs%in%Alldata$trt4)*4
responses<-t(gene1to5)
trtmentFac<-as.factor(trtment)

#############################################################
##
## Perform an MANOVA test using treatment as factor
##
#############################################################

manova1to5<-manova(responses~trtmentFac)
summary(manova1to5)

##############################################################
##
## Comparision among covariance matrices
##
##############################################################

library(biotools)
covcomp<-boxM(responses,trtmentFac)

#############################################################
##
## Perform two-way MANOVA test using treatment and time 
## as factors
##
#############################################################

times<-(sampleIDs%in%Alldata$time_h0)*0+(sampleIDs%in%Alldata$time_h1)*1+(sampleIDs%in%Alldata$time_h2)*2+(sampleIDs%in%Alldata$time_h4)*3+(sampleIDs%in%Alldata$time_h12)*4
timesFac<-as.factor(times)
Twowaymanova<-manova(responses~trtmentFac+timesFac) ## Without interaction
summary(Twowaymanova)
Twowaymanova2<-manova(responses~trtmentFac*timesFac) ## With interaction
summary(Twowaymanova2)