#########################################################
############  Handout #6 for ST476/576  #################
################  Two-way (MANOVA)  #####################
#########################################################


##################################################################
########################## Example 1 #############################
#### Test for intercation and main effects on Generated Data #####
##################################################################

#### Generate random samples of n1=50, n2=50, n3=50 from 6 populations
 ### N(mu_11, Sigma), N(mu_12, Sigma), N(mu_21, Sigma), 
 ### N(mu_22, Sigma), N(mu_31, Sigma), N(mu_32, Sigma) 

library(MASS)
set.seed(11)
mu <- c(0,0,0)

tau1 <- c(0,0,0)
tau2 <- c(0,-1,2)
tau3 <- c(0,1,-2)

beta1 <- c(0,0,0)  ## no factor 2 effect
beta2 <- c(0,0,0)

mu11 <- mu + tau1 + beta1  ## no ineraction
mu12 <- mu + tau1 + beta2
mu21 <- mu + tau2 + beta1
mu22 <- mu + tau2 + beta2
mu31 <- mu + tau3 + beta1
mu32 <- mu + tau3 + beta2


Sigma <- .5*rbind(c(5,  4,  5),
               c(4, 10,  3),
               c(5,  3, 15))

n11 <- 20
n12 <- 20
n21 <- 20
n22 <- 20
n31 <- 20
n32 <- 20

X11 <- mvrnorm(n11, mu11, Sigma)
X12 <- mvrnorm(n11, mu12, Sigma)
X21 <- mvrnorm(n11, mu21, Sigma)
X22 <- mvrnorm(n11, mu22, Sigma)
X31 <- mvrnorm(n11, mu31, Sigma)
X32 <- mvrnorm(n11, mu32, Sigma)


## Column stack the data and add variables to denote the factors
X <- rbind(X11, X12, X21, X22, X31, X32)
fac1 <- factor(c(rep(1,n11+n12), rep(2,n21+n22), rep(3,n31+n32)))
fac2 <- factor(c(rep(1,n11), rep(2,n12), rep(1,n21), rep(2,n22), rep(1,n31), rep(2,n32))) 



#### MANOVA procedure to test for interaction and main effects ####

manova.obj <- manova(X ~ fac1*fac2)
summary(manova.obj, test="Wilks")


## Perform approximate Chi Sqr Tests by direct calculation 
## (just to verify what R is doing)

N <- nrow(X)
p <- ncol(X)
g <- length(unique(fac1))
b <- length(unique(fac2))

## Get SS matrices
SS_fac1 <- summary(manova.obj)$SS[[1]]
SS_fac2 <- summary(manova.obj)$SS[[2]]
SS_int <- summary(manova.obj)$SS[[3]]
SS_res <- summary(manova.obj)$SS[[4]]

## Test for factor1 effect
Lambda_fac1 <- det(SS_res)/det(SS_fac1+SS_res)
a_fac1 <- -(N-g*b-(p+1-(g-1))/2)*log(Lambda_fac1)
p_val <- 1-pchisq(a_fac1, p*(g-1))
p_val

## Test for factor2 effect
Lambda_fac2 <- det(SS_res)/det(SS_fac2+SS_res)
a_fac2 <- -(N-g*b-(p+1-(b-1))/2)*log(Lambda_fac2)
p_val <- 1-pchisq(a_fac2, p*(b-1))
p_val

## Test for interaction effect
Lambda_int <- det(SS_res)/det(SS_int+SS_res)
a_int <- -(N-g*b-(p+1-(g-1)*(b-1))/2)*log(Lambda_int)
p_val <- 1-pchisq(a_int, p*(g-1)*(b-1))
p_val



## Other possiblilities for a test statistic
summary(manova.obj, test="Pillai")
summary(manova.obj, test="Hotelling-Lawley")
summary(manova.obj, test="Roy")



#### Obtain a univariate ANOVA summary for each component ####

summary.aov(manova.obj)


#### 95% Simulataneous CI's for all pairwise diffs between component means ####


## To make simulataneous CI's just need to make sure the treatment variables 
 ## are called factor1, and factor2 respectively, the multivariate response 
 ## matrix is called X, the SS_res matrix is called W, and alpha is set to 
 ## desired value.

factor1 <- fac1
factor2 <- fac2
X <- X
W <- summary(manova.obj)$SS[[4]]
alpha <- .05


## Now the following commands will calculate the CI's for you

N <- nrow(X)
p <- ncol(X)
factor1 <- as.numeric(factor1)
factor2 <- as.numeric(factor2)
g <- length(unique(factor1))
b <- length(unique(factor2))
conf1 <- 1-alpha/(p*g*(g-1))
conf2 <- 1-alpha/(p*b*(b-1))
CI_fac1 <- data.frame(matrix(nrow=p*choose(g,2), ncol=2))
names(CI_fac1) <- c(paste(" ", 100*(1-alpha), "% Lower", sep=""), paste(" ", 100*(1-alpha), "% Upper", sep=""))
CI_fac2 <- data.frame(matrix(nrow=p*choose(b,2), ncol=2))
names(CI_fac2) <- c(paste(" ", 100*(1-alpha), "% Lower", sep=""), paste(" ", 100*(1-alpha), "% Upper", sep=""))

next.row <- 1
for(i in 1:p){
  for(l in 1:(g-1)){
    for(m in (l+1):g){
      x_bar_li <- colMeans(X[factor1==l,])[i]
      x_bar_mi <- colMeans(X[factor1==m,])[i]
      n_l <- sum(factor1==l)
      n_m <- sum(factor1==m)
      lower <- (x_bar_li - x_bar_mi) - qt(conf1, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      upper <- (x_bar_li - x_bar_mi) + qt(conf1, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      descr <- paste("tau_",l, ",", i, " - ", "tau_",m, ",", i, sep="")
      row.names(CI_fac1)[next.row] <- descr
      CI_fac1[next.row,] <- c(format(round(lower,4),nsmall=4), format(round(upper,4),nsmall=4))
      next.row <- next.row+1
    }
  }
}

next.row <- 1
for(i in 1:p){
  for(l in 1:(b-1)){
    for(m in (l+1):b){
      x_bar_li <- colMeans(X[factor2==l,])[i]
      x_bar_mi <- colMeans(X[factor2==m,])[i]
      n_l <- sum(factor2==l)
      n_m <- sum(factor2==m)
      lower <- (x_bar_li - x_bar_mi) - qt(conf2, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      upper <- (x_bar_li - x_bar_mi) + qt(conf2, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      descr <- paste("beta_",l, ",", i, " - ", "beta_",m, ",", i, sep="")
      row.names(CI_fac2)[next.row] <- descr
      CI_fac2[next.row,] <- c(format(round(lower,4),nsmall=4), format(round(upper,4),nsmall=4))
      next.row <- next.row+1
    }
  }
}

## Print the CI's
CI_fac1
CI_fac2

## Note: the CI's for factor 2 are not meaningful in this case because the
## overall MANOVA test conclusion was to not reject H_0: no factor 2 effect



#### Create Interaction Plots ####

par(mfrow=c(3,1))
interaction.plot(fac1, fac2, X[,1], main="Interaction Plot", ylab="1st Component Means")

interaction.plot(fac1, fac2, X[,2], main="Interaction Plot", ylab="2nd Component Means")

interaction.plot(fac1, fac2, X[,3], main="Interaction Plot", ylab="3rd Component Means")
par(mfrow=c(1,1))



#### Create Main Effects Plots ####

## Only needed for 2nd and 3rd components of factor 1 as indicated by CI's

par(mfrow=c(2,1))
interaction.plot(fac1, rep(1,length(fac1)), X[,2], legend=F, main="Main Effects Plot for Factor 1", ylab="2nd Component Means")

interaction.plot(fac1, rep(1,length(fac1)), X[,3], legend=F, main="Main Effects Plot for Factor 1", ylab="3rd Component Means")
par(mfrow=c(1,1))



## Could also use boxplots
par(mfrow=c(2,1))
plot(fac1, X[,2], main="Boxplot for Factor 1, 2nd component")
plot(fac1, X[,3], main="Boxplot for Factor 1, 3rd component")
par(mfrow=c(1,1))










###########################################################
####################### Example 2 #########################
### Test of MANOVA hyp, and pairwise CI's on Reflectance Data ###
###########################################################

#### Read in data ####
ex.data <- read.table(file="C://jenn/teaching/stat4765762012/data/reflectance_data.txt", header=T)
nrow(ex.data)
ex.data[1:10,]

#### MANOVA for Reflect = Species*Time ####

Reflect <- as.matrix(ex.data[,1:2])

Reflect_t <- cbind(log(ex.data[,1]), ex.data[,2]^(-1))

Species <- factor(ex.data[,3])
Species
Time <- factor(ex.data[,4])
Time

manova.obj <- manova(Reflect_t ~ Species*Time) #same as (Reflect_t ~ Species
#+ Time + Species*Time
summary(manova.obj,test="Wilks")


#### Check MANOVA assumptions ####

## (1) Assess marginal Normality of each components residuals
##    via histograms and QQ-plots 

resid <- manova.obj$residuals

par(mfrow=c(2,2))  

hist(Reflect_t[,1], main=names(ex.data)[1], xlab='X_1')  
qqnorm(resid[,1])  
qqline(resid[,1])  

hist(Reflect_t[,2], main=names(ex.data)[2], xlab='X_2')  
qqnorm(resid[,2])  
qqline(resid[,2])  

par(mfrow=c(1,1))


## (2) Assess bivariate Normaility of each pair
pairs(resid)



#### Obtain a univariate ANOVA summary for each component ####

summary.aov(manova.obj)


#### 90% Simulataneous CI's for all pairwise diffs between component means ####


## To make simulataneous CI's just need to make sure the treatment variables 
 ## are called factor1, and factor2 respectively, the multivariate response 
 ## matrix is called X, the SS_res matrix is called W, and alpha is set to 
 ## desired value.

##note that if interaction is significant, there is no need to do the pairwise comparison of main factors. Instead, we 
##will do pairwise comparisons of the interaction. I leave the coding to you.
factor1 <- Species
factor2 <- Time
X <- Reflect_t
W <- summary(manova.obj)$SS[[4]]
alpha <- .10


## Now the following commands will calculate the CI's for you

N <- nrow(X)
p <- ncol(X)
factor1 <- as.numeric(factor1)
factor2 <- as.numeric(factor2)
g <- length(unique(factor1))
b <- length(unique(factor2))
conf1 <- 1-alpha/(p*g*(g-1))
conf2 <- 1-alpha/(p*b*(b-1))
CI_fac1 <- data.frame(matrix(nrow=p*choose(g,2), ncol=2))
names(CI_fac1) <- c(paste(" ", 100*(1-alpha), "% Lower", sep=""), paste(" ", 100*(1-alpha), "% Upper", sep=""))
CI_fac2 <- data.frame(matrix(nrow=p*choose(b,2), ncol=2))
names(CI_fac2) <- c(paste(" ", 100*(1-alpha), "% Lower", sep=""), paste(" ", 100*(1-alpha), "% Upper", sep=""))

next.row <- 1
for(i in 1:p){
  for(l in 1:(g-1)){
    for(m in (l+1):g){
      x_bar_li <- colMeans(X[factor1==l,])[i]
      x_bar_mi <- colMeans(X[factor1==m,])[i]
      n_l <- sum(factor1==l)
      n_m <- sum(factor1==m)
      lower <- (x_bar_li - x_bar_mi) - qt(conf1, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      upper <- (x_bar_li - x_bar_mi) + qt(conf1, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      descr <- paste("tau_",l, ",", i, " - ", "tau_",m, ",", i, sep="")
      row.names(CI_fac1)[next.row] <- descr
      CI_fac1[next.row,] <- c(format(round(lower,4),nsmall=4), format(round(upper,4),nsmall=4))
      next.row <- next.row+1
    }
  }
}

next.row <- 1
for(i in 1:p){
  for(l in 1:(b-1)){
    for(m in (l+1):b){
      x_bar_li <- colMeans(X[factor2==l,])[i]
      x_bar_mi <- colMeans(X[factor2==m,])[i]
      n_l <- sum(factor2==l)
      n_m <- sum(factor2==m)
      lower <- (x_bar_li - x_bar_mi) - qt(conf2, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      upper <- (x_bar_li - x_bar_mi) + qt(conf2, N-g*b)*sqrt(W[i,i]/(N-g*b)*(1/n_l+1/n_m))
      descr <- paste("beta_",l, ",", i, " - ", "beta_",m, ",", i, sep="")
      row.names(CI_fac2)[next.row] <- descr
      CI_fac2[next.row,] <- c(format(round(lower,4),nsmall=4), format(round(upper,4),nsmall=4))
      next.row <- next.row+1
    }
  }
}

## Print the CI's
CI_fac1
CI_fac2




#### Create Interaction Plots ####

par(mfrow=c(2,1))
interaction.plot(Species, Time, X[,1], main="Interaction Plot", ylab="1st Component Means")

interaction.plot(Species, Time, X[,2], main="Interaction Plot", ylab="2nd Component Means")
par(mfrow=c(1,1))



#### Create Main Effects Plots for Species ####

par(mfrow=c(2,1))
interaction.plot(Species, rep(1,length(Species)), X[,1], legend=F, main="Main Effects Plot for Species", ylab="1st Component Means")

interaction.plot(Species, rep(1,length(Species)), X[,2], legend=F, main="Main Effects Plot for Species", ylab="2nd Component Means")
par(mfrow=c(1,1))


#### Create Main Effects Plots for Time ####

par(mfrow=c(2,1))
interaction.plot(Time, rep(1,length(Time)), X[,1], legend=F, main="Main Effects Plot for Time", ylab="1st Component Means")

interaction.plot(Time, rep(1,length(Time)), X[,2], legend=F, main="Main Effects Plot for Time", ylab="2nd Component Means")
par(mfrow=c(1,1))