#HOTELLING'S T2 TEST FOR PAIRED DESIGN:

#LOAD DATA & HYPOTHESIS:
M=read.table("c:/DATA/Multivariate/JW6-1-longform.txt")
M   #DATA TABLE
  
#SEPARATING GROUPS:           

X=M[(M$Lab=='1'),1:2]
Y=M[(M$Lab=='2'),1:2]
X
Y

#=============================================
#FUNCTION FOR HOTELLING'S T2 TEST FOR PAIRED DESIGN:
#X     = dataset 1
#Y     = dataset 2
#mu0   = hypothesis vector
#alpha = alpha of the test

PD.Hotelling.T2 <- function(X,Y,mu0,alpha)
{
  D=X-Y      #subtract matrices
  X=D        #uses OS code for X=D from this point on
  n=nrow(X)
  p=ncol(X)
  Xbar=colMeans(X)
  S=cov(X)

  Tsq=n%*%((Xbar-mu0)%*%solve(S)%*%(Xbar-mu0))
  Tsq

  Sigmahat=((n-1)/n)*S
  XX=as.matrix(sweep(X,MARGIN=2,mu0,FUN="-"))
  Sigma0hat=t(XX)%*%XX
  Wilks.lambda=(det(n*(Sigmahat))/(det(Sigma0hat)))^(n/2)
  Wilks.lambda

  C=(((n-1)*p)/(n-p))*qf(1-alpha,p,n-p)
  C

  F=Tsq*(n-p)/((n-1)*p)
  df1=p
  df2=n-p

  P=1-pf(F,p,n-p)
  P

  DiscriminantCoef=solve(S)%*%Xbar

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

  FF=(p*(n-1))/(n*(n-p))
  CT=sqrt(FF*qf(1-alpha,p,n-p))
  CT
  ct=qt(1-alpha/2,n-1)
  ct
  cb=qt(1-alpha/(2*p),n-1)
  cb

  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/n)*C*lambda[i])
    sigma[i]=S[i,i]


    T2.CI.lower[i]=Xbar[i]-CT*sqrt(sigma[i])
    T2.CI.upper[i]=Xbar[i]+CT*sqrt(sigma[i])

    Univar.CI.lower[i]=Xbar[i]-ct*sqrt(sigma[i]/n)
    Univar.CI.upper[i]=Xbar[i]+ct*sqrt(sigma[i]/n)

    Bonf.CI.lower[i]=Xbar[i]-cb*sqrt(sigma[i]/n)
    Bonf.CI.upper[i]=Xbar[i]+cb*sqrt(sigma[i]/n)

  }

  CI=cbind(T2.lower=T2.CI.lower,B.lower=Bonf.CI.lower,t.lower=Univar.CI.lower,Mean=Xbar,
         t.upper=Univar.CI.upper,B.upper=Bonf.CI.upper,T2.upper=T2.CI.upper)

  result=list(hypothesis.vector=mu0,HotellingT2=Tsq,WilksLambda=Wilks.lambda,CV=C,Prob=P,
           CI=CI,eigenvectors=eigenvec,eigenvalues=lambda,mean.vector=Xbar,CT=CT,D=D)
  class(result) = 'OS.Hotelling.T2'
  result
  cat("\n","Paired Design Hotelling's T2","\n","\n",    
     "Hypothesis Vector:       (",mu0,")","\n",
     "Difference Vector        (",Xbar,")","\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)
  cat("\n")

  return(result)

}
#=============================================

RES=PD.Hotelling.T2(X,Y,mu0=c(0,0),alpha=0.05)

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

#FUNCTION PLOT CI:
#R = return from function OS.Hotelling.T2 as a list
#X = variable to plot along the X axis
#Y = variable to plot along the Y axis

plot.CI <- function(R,X,Y)
{

  mu0=R$hypothesis.vector
  Xbar=R$mean.vector
  
  plot(RES$D[,X],RES$D[,Y],type='n',xlab=X,ylab=Y)

  points(RES$D[,X],RES$D[,Y],pch=19,cex=0.5,col='blue')
  points(mu0[X],mu0[Y],pch=2,col='red')
  points(Xbar[X],Xbar[Y],pch=3,col='blue')

  abline(v=R$CI[X,1],col='red')
  abline(v=R$CI[X,7],col='red')
  abline(v=R$CI[X,2],col='brown',lty=3)
  abline(v=R$CI[X,6],col='brown',lty=3)
  abline(v=R$CI[X,3],col='brown',lty=2)
  abline(v=R$CI[X,5],col='brown',lty=2)

  abline(h=R$CI[Y,1],col='red')
  abline(h=R$CI[Y,7],col='red')
  abline(h=R$CI[Y,2],col='brown',lty=3)
  abline(h=R$CI[Y,6],col='brown',lty=3)
  abline(h=R$CI[Y,3],col='brown',lty=2)
  abline(h=R$CI[Y,5],col='brown',lty=2)

}

plot.CI(RES,1,2)