wants  =  c("car", "mvtnorm","Hotelling")
has    =  wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])
library(mvtnorm)
library(Hotelling)

################################################################################################
# Example 1: simulated data
# Source: http://www.uni-kiel.de/psychologie/rexrepos/posts/multMANOVA.html
################################################################################################
set.seed(123)
p     =  3
n     =  c(15, 25, 20)
Sigma =  matrix(c(16,-2, -2,9), byrow=TRUE, ncol=2)
mu1   =  c(-4,  4)
mu2   =  c( 3,  3)
mu3   =  c( 1, -1)

y1    =  round(rmvnorm(n[1], mean=mu1, sigma=Sigma))
y2    =  round(rmvnorm(n[2], mean=mu2, sigma=Sigma))
y3    =  round(rmvnorm(n[3], mean=mu3, sigma=Sigma))

data1 =  data.frame(y=rbind(y1, y2, y3),IV=factor(rep(1:p,n)))

manova1  =  manova(cbind(y.1,y.2) ~ IV, data=data1)

summary(manova1, test="Wilks")

summary(manova1, test="Roy")

summary(manova1, test="Pillai")

summary(manova1, test="Hotelling-Lawley")


################################################################################################
# Example 2: Example on producing plastic film from Krzanowski (1988, p. 381)
# Source: R documentation for summary.manova
# Krzanowski (1988) Principles of Multivariate Analysis. A User's Perspective. Oxford
################################################################################################

tear  =  c(6.5, 6.2, 5.8, 6.5, 6.5, 6.9, 7.2, 6.9, 6.1, 6.3,
          6.7, 6.6, 7.2, 7.1, 6.8, 7.1, 7.0, 7.2, 7.5, 7.6)
          
gloss  =  c(9.5, 9.9, 9.6, 9.6, 9.2, 9.1, 10.0, 9.9, 9.5, 9.4,
           9.1, 9.3, 8.3, 8.4, 8.5, 9.2, 8.8, 9.7, 10.1, 9.2)
           
opacity  =  c(4.4, 6.4, 3.0, 4.1, 0.8, 5.7, 2.0, 3.9, 1.9, 5.7,
             2.8, 4.1, 3.8, 1.6, 3.4, 8.4, 5.2, 6.9, 2.7, 1.9)
             
Y  =  cbind(tear, gloss, opacity)

rate  =  factor(gl(2,10), labels = c("Low", "High"))

additive  =  factor(gl(2, 5, length = 20), labels = c("Low", "High"))

fit  =  manova(Y ~ rate * additive)

summary.aov(fit)             # univariate ANOVA tables

summary(fit, test = "Wilks") # ANOVA table of Wilks' lambda

summary(fit)                # same F statistics as single-df terms




################################################################################################
# Exemplo3: Notas da aula 2
################################################################################################
p = 2
n = c(10,10)

y1 = matrix(
c(5,3,4.5,3.2,6,3.5,6,4.6,6.2,5.6,6.9,5.2,6.8,6,5.3,5.5,6.6,7.3,7.3,6.5),n[1],p,byrow=TRUE)

y2 = matrix(
c(4.6,4.9,4.9,5.9,4,4.1,3.8,5.4,6.2,6.1,5,7,5.3,4.7,7.1,6.6,5.8,7.8,6.8,8.0),n[2],p,byrow=TRUE)


data2    =  data.frame(y=rbind(y1, y2),IV=factor(rep(1:p,n)))

manova2  =  manova(cbind(y.1,y.2) ~ IV, data=data2)

summary(manova2, test="Hotelling-Lawley")

hotelling2 = hotelling.test(y1,y2)