fat <- read.table('DATA/mlrm-fat.txt', header=T)
str(fat)

fat2 <- fat[-36, -3]

colMeans(fat2)

var(fat2)

library(MVN)
mvn(fat2, mvnTest = 'mardia', multivariatePlot = 'qq')$multivariateNormality

mvn(fat2, mvnTest = 'hz')$multivariateNormality

mu0 <- c(60.8, 167.9, 13, 21.5, 22, 25)
library(ICSNP)
## Loading required package: mvtnorm
## Loading required package: ICS
HotellingsT2(fat2, mu=mu0)

alpha.korig <- 0.05 / 6
alpha.korig
## [1] 0.008333333

t.test(fat2$body.W, mu=mu0[1])

t.test(fat2$body.H, mu=mu0[2])

t.test(fat2$rib.F, mu=mu0[3])

t.test(fat2$abdo.F, mu=mu0[4])

t.test(fat2$hip.F, mu=mu0[5])

t.test(fat2$quad.H, mu=mu0[6])

digits <- read.table('DATA/d2d4.txt', header=T)
str(digits)

digits <- digits[,-1]
colMeans(digits[digits$sex=='f', 2:3])

var(digits[digits$sex=='f', 2:3])

colMeans(digits[digits$sex=='m', 2:3])

var(digits[digits$sex=='m', 2:3])

library(MVN)
par(mfrow=c(1,2))
mvn(digits, subset='sex', mvnTest = 'mardia', multivariatePlot = 'qq')$multivariateNormality

mvn(digits, subset='sex', mvnTest = 'hz')$multivariateNormality

library(biotools)
boxM(digits[,2:3], grouping=digits$sex)

library(ICSNP)
HotellingsT2(digits[digits$sex=="f",2:3], digits[digits$sex=="m",2:3])

n1 <- table(digits$sex)[1]
n2 <- table(digits$sex)[2]
n <- n1 + n2
k <- 2 #pocet promennych
mu1 <- colMeans(digits[digits$sex=="f",2:3])
mu2 <- colMeans(digits[digits$sex=="m",2:3])
var1 <- diag(cov(digits[digits$sex=="f",2:3]))
var2 <- diag(cov(digits[digits$sex=="m",2:3]))
var <- ( (n1-1)*var1 + (n2-1)*var2 )/(n-2)
F.stat <- n1*n2*(n-k-1) * (mu1-mu2)^2 /(var*n*k*(n-2))
p.hodnota <- 1-pf(F.stat, k, n-k-1)
kvantil <- qf(0.95, k, n-k-1)
tab <- round(rbind(F.stat,p.hodnota, kvantil),digits=4)
rownames(tab) <- c("F","p-hodnota", "kvantil")
tab