cranio <- read.csv('DATA/Howell.csv',header=T, na.strings='0')
data <- cranio[cranio$Sex == 'M' & cranio$Population %in% c('AINU', 'N JAPAN', 'PERU'),
               c('Population', 'XFB', 'ZYB', 'ZMB')]
data$Population <- factor(data$Population)

table(data$Population)

colMeans(data[data$Population=='AINU',-1])

colMeans(data[data$Population=='N JAPAN',-1])

colMeans(data[data$Population=='PERU',-1])

var(data[data$Population=='AINU',-1])

var(data[data$Population=='N JAPAN',-1])

var(data[data$Population=='PERU',-1])

par(mfrow=c(1,3))
boxplot(data$XFB ~ data$Population, ylab='frontal breadth (mm)')
boxplot(data$ZYB ~ data$Population, ylab='bizygomatic breadth (mm)')
boxplot(data$ZMB ~ data$Population, ylab='bimaxillary breadth (mm)')

library(MVN)
par(mfrow=c(1,3))
mvn(data, subset='Population', mvnTest = 'mardia', multivariatePlot = 'qq')$multivariateNormality

mvn(data, subset='Population', mvnTest = 'hz')$multivariateNormality

library(biotools)
boxM(data[,-1], grouping=data$Population)

model <- manova(as.matrix(data[,-1]) ~ data$Population)
summary(model,test='Wilks')

summary(model,test='Pillai')

summary(model,test='Hotelling-Lawley')

summary(model,test='Roy')

SSB <- summary(model,test='Wilks')$SS[[1]]
SSE <- summary(model,test='Wilks')$SS[[2]]
SST <- SSB + SSE
n <- nrow(data)
k <- 3 # pocet promennych
r <- 3 # pocet skupin
const <- (n- (k+r)/2 -1)
K <- -const * log(diag(SSE)/diag(SST))
K
## XFB ZYB ZMB
## 24.496114 19.807912 5.739724
kvantil <- qchisq(0.95, df=k*(r-1))
kvantil
## [1] 12.59159

alpha.korig <- 0.05/ choose(r,2)
alpha.korig
## [1] 0.01666667
library(ICSNP)
HotellingsT2(data[data$Population=='AINU',-1], data[data$Population=='N JAPAN',-1])

HotellingsT2(data[data$Population=='AINU',-1], data[data$Population=='PERU',-1])

HotellingsT2(data[data$Population=='N JAPAN',-1], data[data$Population=='PERU',-1])

alpha.korig2 <- 0.05 / (k*r*(r-1)/2)
alpha.korig2
## [1] 0.005555556
# AINU vs N JAPAN
t.test(data[data$Population=='AINU',2], data[data$Population=='N JAPAN',2])

t.test(data[data$Population=='AINU',3], data[data$Population=='N JAPAN',3])

t.test(data[data$Population=='AINU',4], data[data$Population=='N JAPAN',4])

# AINU vs PERU
t.test(data[data$Population=='AINU',2], data[data$Population=='PERU',2])
t.test(data[data$Population=='AINU',3], data[data$Population=='PERU',3])
t.test(data[data$Population=='AINU',4], data[data$Population=='PERU',4])

# N JAPAN vs PERU
t.test(data[data$Population=='N JAPAN',2], data[data$Population=='PERU',2])
t.test(data[data$Population=='N JAPAN',3], data[data$Population=='PERU',3])
t.test(data[data$Population=='N JAPAN',4], data[data$Population=='PERU',4])