#TWO SAMPLE PROFILE ANALYSIS:

#=============================================
#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  

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

Xbar=colMeans(X)
Xbar
Ybar=colMeans(Y)
Ybar

S1=cov(X)
S1
S2=cov(Y)
S2

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

#=============================================
#PLOT PROFILE:

minX=min(Xbar)
minY=min(Ybar)
min=min(minX,minY)
maxX=max(Xbar)
maxY=max(Ybar)
max=max(maxX,maxY)

maxmin=cbind(maxX,maxY)

index=c(0,p)

plot(index,maxmin,xlim=c(1,p),ylim=c(min,max),type='n',xlab='variable',ylab='mean')

lines(seq(1:4),Xbar,pch=19,col='blue')
lines(seq(1:4),Ybar,pch=19,col='red')

points(seq(1:4),Xbar,pch=19,col='blue')
points(seq(1:4),Ybar,pch=19,col='red')

#=============================================
#CONTRAST MATRIX:
C=matrix(c(-1,1,0,0,0,-1,1,0,0,0,-1,1),nrow=3,ncol=4,byrow=T)
C

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

#FUNCTION FOR TWO SAMPLE PROFILE ANALYSIS:
#X     = dataset 1
#Y     = dataset 2
#alpha = alpha of the test

T2Profile.Analysis <- function(X,Y,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)
  C=matrix(c(-1,1,0,0,0,-1,1,0,0,0,-1,1),nrow=3,ncol=4,byrow=T)

 
  Tsq01=t(C%*%d)%*%solve(((1/n1)+(1/n2))*(C%*%Spooled%*%t(C)))%*%(C%*%d)
  Tsq01
  c01=(((n1+n2-2)*(p-1))/(n1+n2-p))*qf(1-alpha,p-1,n1+n2-p)
  c01
  P01=1-pf(Tsq01,p-1,n1+n2-p)
  P01
 
  l=rep(1,times=p)
  Tsq02=(t(l)%*%d)%*%solve(((1/n1)+(1/n2))*(t(l)%*%Spooled)%*%l)%*%(l%*%d)
  Tsq02
  t02=(t(l)%*%d)/sqrt((t(l)%*%Spooled%*%l)*((1/n1)+(1/n2)))
  t02
  c02T=qf(1-alpha,1,n1+n2-2)
  c02T
  c02t=qt(1-alpha,n1+n2-2)
  c02t
  P02=1-pf(Tsq02,1,n1+n2-2)
  P02

  N=n1+n2
  UD=rbind(X,Y)
  UDbar=colMeans(UD)
  SUD=cov(UD)
  UDbar
  SUD
  Tsq03JW=N*(t(C%*%UDbar)%*%solve(C%*%SUD%*%t(C))%*%(C%*%UDbar))
  Tsq03JW
  Tsq03AR=N*(t(C%*%UDbar)%*%solve(C%*%Spooled%*%t(C))%*%(C%*%UDbar))
  Tsq03AR
  c03=(((N-1)*(p-1))/(N-p+1))*qf(1-alpha,p-1,N-p+1)
  c03
  P03JW=1-pf(Tsq03JW,p-1,N-p+1)
  P03JW
  P03AR=1-pf(Tsq03AR,p-1,N-p+1)
  P03AR

  Tsq03X=n1*(t(C%*%Xbar)%*%solve(C%*%S1%*%t(C))%*%(C%*%Xbar))
  Tsq03X
  Tsq03Y=n1*(t(C%*%Ybar)%*%solve(C%*%S2%*%t(C))%*%(C%*%Ybar))
  Tsq03Y

  c03X=(((n1-1)*(p-1))/(n1-p+1))*qf(1-alpha,p-1,n1-1)
  c03X
  c03Y=(((n2-1)*(p-1))/(n2-p+1))*qf(1-alpha,p-1,n2-1)
  c03Y
  
  P03X=1-pf(Tsq03X,p-1,n1-1)
  P03X
  P03Y=1-pf(Tsq03Y,p-1,n2-1)
  P03Y

  result=list(mean.difference=d,Cov.X=S1,Cov.Y=S2,Cov.Pooled=Spooled,
            Tsq01=Tsq01,c01=c01,P01=P01,
            Tsq02=Tsq02,c02T=c02T,P02=P02,
            t02=t02,c02t=c02t,
            Tsq03JW=Tsq03JW,Tsq03AR=Tsq03AR,c03=c03,
            P03JW=P03JW,P03AR=P03AR,
            Tsq03X=Tsq03X,Tsq03Y=Tsq03Y,
            c03X=c03X,c03Y=c03Y,
            P03X=P03X,P03Y=P03Y)
  class(result)='Hotelling'
            
  cat("\n","Hotelling's T2 Profile Analysis","\n","\n",
      "Test for Parallel Profiles:","\n","\n",
      "Hotelling's T2 = ",Tsq01,"  critical value = ",c01,"  Probability = ",P01,"\n","\n",

      "Test for Coincident Profiles:","\n","\n",
      "Hotelling's T2 = ",Tsq02,"  critical value = ",c02T,"  Probability = ",P02,"\n",
      "Student's t = ",t02,"  critical value = ",c02t,"\n","\n",

      "Test for Combined Level Profile:","\n","\n",
      "Hotelling's T2 (JW) = ",Tsq03JW,"  critical value = ",c03,"  Probability = ",P03JW,"\n",
      "Hotelling's T2 (AR) = ",Tsq03AR,"  critical value = ",c03,"  Probability = ",P03AR,"\n","\n",

      "Test for Separate Level Profiles:","\n","\n",
      "Hotelling's T2 for X = ",Tsq03X,"  critical value = ",c03X,"  Probability = ",P03X,"\n",
      "Hotelling's T2 for Y = ",Tsq03Y,"  critical value = ",c03Y,"  Probability = ",P03Y,"\n","\n")


  return(result)
}     

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

RES=T2Profile.Analysis(X,Y,alpha=0.01)
RES=T2Profile.Analysis(X,Y,alpha=0.005)

#=============================================
#NOTE To compare results with Worksheet MHT 050
# and AR pp. 167-168, use:
#         alpha01=0.01
#         alpha02=0.005
#         alpha03=0.01

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























  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)
  C=matrix(c(-1,1,0,0,0,-1,1,0,0,0,-1,1),nrow=3,ncol=4,byrow=T)

 
  Tsq01=t(C%*%d)%*%solve(((1/n1)+(1/n2))*(C%*%Spooled%*%t(C)))%*%(C%*%d)
  Tsq01
  c01=(((n1+n2-2)*(p-1))/(n1+n2-p))*qf(1-alpha,p-1,n1+n2-p)
  c01
  P01=1-pf(Tsq01,p-1,n1+n2-p)
  P01
 
  l=rep(1,times=p)
  Tsq02=(t(l)%*%d)%*%solve(((1/n1)+(1/n2))*(t(l)%*%Spooled)%*%l)%*%(l%*%d)
  Tsq02
  t02=(t(l)%*%d)/sqrt((t(l)%*%Spooled%*%l)*((1/n1)+(1/n2)))
  t02
  c02T=qf(1-alpha,1,n1+n2-2)
  c02T
  c02t=qt(1-alpha,n1+n2-2)
  c02t
  P02=1-pf(Tsq02,1,n1+n2-2)
  P02

  N=n1+n2
  UD=rbind(X,Y)
  UDbar=colMeans(UD)
  SUD=cov(UD)
  UDbar
  SUD
  Tsq03JW=N*(t(C%*%UDbar)%*%solve(C%*%SUD%*%t(C))%*%(C%*%UDbar))
  Tsq03JW
  Tsq03AR=N*(t(C%*%UDbar)%*%solve(C%*%Spooled%*%t(C))%*%(C%*%UDbar))
  Tsq03AR
  c03=(((N-1)*(p-1))/(N-p+1))*qf(1-alpha,p-1,N-p+1)
  c03
  P03JW=1-pf(Tsq03JW,p-1,N-p+1)
  P03JW
  P03AR=1-pf(Tsq03AR,p-1,N-p+1)
  P03AR

  Tsq03X=n1*(t(C%*%Xbar)%*%solve(C%*%S1%*%t(C))%*%(C%*%Xbar))
  Tsq03X
  Tsq03Y=n1*(t(C%*%Ybar)%*%solve(C%*%S2%*%t(C))%*%(C%*%Ybar))
  Tsq03Y

  c03X=(((n1-1)*(p-1))/(n1-p+1))*qf(1-alpha,p-1,n1-1)
  c03X
  c03Y=(((n2-1)*(p-1))/(n2-p+1))*qf(1-alpha,p-1,n2-1)
  c03Y
  
  P03X=1-pf(Tsq03X,p-1,n1-1)
  P03X
  P03Y=1-pf(Tsq03Y,p-1,n2-1)
  P03Y