#HOTELLING'S T2 TEST FOR TWO POPULATIONS WITH EQUAL COVARIANCE MATRICES
#AND SMALL SAMPLE SIZES

#=============================================
#LOAD DATA:
M=read.table("c:/DATA/Multivariate/AR-T5-1-longform.txt")
M   #DATA TABLE
X=M[(M$sex=='1'),2:5]
Y=M[(M$sex=='2'),2:5]
X
Y  

d0=vector(length=ncol(X))
alpha=0.05

#=============================================
#FUNCTION FOR HOTELLING'S T2 TEST FOR TWO POPULATIONS EQUAL WITH VARIANCE:
#X     = dataset 1
#Y     = dataset 2
#d0   = hypothesis vector
#alpha = alpha of the test

TwoPop.Hotelling.T2 <- function(X,Y,d0,alpha)
{

  n1=nrow(X)
  n2=nrow(Y)
  p=ncol(X)

  Xbar1=colMeans(X)
  Xbar2=colMeans(Y)
  d=Xbar1-Xbar2

  S1=cov(X)
  S2=cov(Y)
  Spooled=((n1-1)*S1+(n2-1)*S2)/(n1+n2-2)

  Tsq=t(d-d0)%*%solve(((1/n1)+(1/n2))*Spooled)%*%(d-d0)

  Wilks.lambda=1/((Tsq/(n1+n2-2))+1)

  C=(((n1+n2-2)*p)/(n1+n2-p-1))*qf(1-alpha,p,n1+n2-p-1)

  F=Tsq*((n1+n2-p-1)/((n1+n2-2)*p))
  df1=p
  df2=n1+n2-p-1

  P=1-pf(F,p,n1+n2-p-1)

  DiscriminantCoef=solve(Spooled)%*%d

  lambda=eigen(Spooled)$values
  eigenvec=eigen(Spooled)$vectors

  ct=qt(1-alpha/2,n1+n2-2)
  cb=qt(1-alpha/(2*p),n1+n2-2)

  sigma=as.numeric(vector(length=p))
  T2.CI.lower=as.numeric(vector(length=p))
  T2.CI.upper=as.numeric(vector(length=p))
  Univar.CI.lower=as.numeric(vector(length=p))
  Univar.CI.upper=as.numeric(vector(length=p))
  Bonf.CI.lower=as.numeric(vector(length=p))
  Bonf.CI.upper=as.numeric(vector(length=p))
  L=vector(length=p)

  for (i in 1:p)
  {
    L[i]=sqrt(((1/n1)+(1/n2))*C*lambda[i])

    sigma[i]=Spooled[i,i]

    T2.CI.lower[i]=d[i]-sqrt(C*((1/n1)+(1/n2))*sigma[i])
    T2.CI.upper[i]=d[i]+sqrt(C*((1/n1)+(1/n2))*sigma[i])

    Univar.CI.lower[i]=d[i]-ct*sqrt(((1/n1)+(1/n2))*sigma[i])
    Univar.CI.upper[i]=d[i]+ct*sqrt(((1/n1)+(1/n2))*sigma[i])

    Bonf.CI.lower[i]=d[i]-cb*sqrt(((1/n1)+(1/n2))*sigma[i])
    Bonf.CI.upper[i]=d[i]+cb*sqrt(((1/n1)+(1/n2))*sigma[i])
  }
 
  CI=cbind(T2.lower=T2.CI.lower,B.lower=Bonf.CI.lower,t.lower=Univar.CI.lower,Mean=d,
         t.upper=Univar.CI.upper,B.upper=Bonf.CI.upper,T2.upper=T2.CI.upper)

  result=list(hypothesis.vector=d0,HotellingT2=Tsq,WilksLambda=Wilks.lambda,CV=C,Prob=P,
           CI=CI,eigenvectors=eigenvec,eigenvalues=lambda,mean.vector=d,
           Cov.X=S1,Cov.Y=S2,Cov.Pooled=Spooled,L=L)

  class(result) = 'Hotelling'
 
  cat("\n","Two Population Hotelling's T2 with Equal Covariance Matrices","\n",
     "                    for small sample sizes","\n","\n",   
     "Hypothesis Vector:         (",d0,")","\n",
     "Difference vector X-Y:     (",d,")","\n",
     "Discriminant vector  :     (",DiscriminantCoef,")","\n",
     "T2 Ellipsoid half lengths: (",L,")","\n",
     "Hotelling's T2 Statistic:   ",Tsq,"\n",
     "Equivalent F Statistic :    ",F,"\n",
     "F degrees of freedom:      (",df1,df2,")","\n",
     "Wilks's Lambda:             ",Wilks.lambda,"\n",
     "alpha:                      ",alpha,"\n",
     "Critical Value:             ",C,"\n",
     "Probability:                ",P,"\n","\n", 
  "Confidence Intervals: T2 - Bonferroni - t - Mean - t - Bonferroni - T2","\n","\n")
  
  print(CI)
  return(result)
}
#=============================================

RES=TwoPop.Hotelling.T2(X,Y,d0=c(0,0,0,0),alpha=0.05)

RES$Cov.Pooled

#=============================================

library(ICSNP)
HotellingsT2(X,Y,mu=c(0,0,0,0),test="f")

#=============================================

library(rrcov)
T2.test(X,Y,mu=c(0,0,0,0),conf.level=0.95,method="c")

#=============================================
#USING MANOVA:
 
attach(M)
YY=cbind(y1,y2,y3,y4)
detach(M)

FIT=manova(YY~sex,data=M)
summary(FIT,test='Wilks')