#install.packages("ICSNP")

library(ICSNP)
data(pulmonary)


#Hotelling T^2, unknown Sigma
head(pulmonary)

pplot(pulmonary)
cor(pulmonary)

cov(pulmonary)

# Ho -> no difference -> 0
t.test(pulmonary$FVC)
t.test(pulmonary$FEV)
t.test(pulmonary$CC)

# now multivariate
HotellingsT2(pulmonary)

1- pf(df1 = 3, df2 = 9, 3.8231)

mean_vec <- apply(pulmonary, mean, MARGIN = 2)

t2 <- 12*(9)/(3*11)*t(mean_vec) %*% solve(cov(pulmonary)) %*% mean_vec

1- pf(df1 = 3, df2 = 9, 3.8231)


# now with different mu_0

HotellingsT2(pulmonary, mu = c(-0.1,-0.1,0)  )
mu_vec <- c(-0.1,-0.1,0)

ts <- 12*(9)/(3*11)*t(mean_vec - mu_vec) %*% 
        solve(cov(pulmonary)) %*% (mean_vec - mu_vec)

1 - pf(ts, df1 = 3, df2 = 9)

#------------------

# two sample test
X <- rmvnorm(n = 20, c(0,0,0,0), diag(1:4) + 1)
Y <- rmvnorm(30, c(1,1,0.5,0.5), diag(1:4) + 1)

apply(X, mean, MARGIN = 2)
apply(Y, mean, MARGIN = 2)

HotellingsT2(X, Y)

# or 
Z <- rbind(X, Y)

g <- factor(rep(c(1,2),  c(20,30)))
HotellingsT2(Z ~ g)

1-pf(3.4448, df1 = 4, df2 = 45)




### loop multiple testing

B <- 1000
pvals <- numeric(B)

for(i in 1:B)
{
  
  # two sample test
  X <- rmvnorm(20, c(0,0,0,0), diag(1:4) + 1)
  Y <- rmvnorm(30, c(0,0,0,0), diag(1:4) + 2)
  
  pvals[i] <- HotellingsT2(X, Y)$p.value
  
}

table(pvals < 0.05)
hist(pvals, breaks = 20)