#===================== 輪入資料 ===========================

a=c(30,7)
b=c(34,8)
c=c(40,6)
d=c(50,5)
e=c(60,4)


a1=c(d,e)
a2=c(c,e)
a3=c(c,d)
a4=c(b,e)
a5=c(b,d)
a6=c(b,c)
a7=c(a,e)
a8=c(a,d)
a9=c(a,b)
a10=c(a,c)

b1=c(a,b,c)
b2=c(a,b,d)
b3=c(a,b,e)
b4=c(a,c,d)
b5=c(a,c,e)
b6=c(a,d,e)
b7=c(b,c,d)
b8=c(b,c,e)
b9=c(c,d,e)
b10=c(b,d,e)

N=10      #C5取2 or C5取3
k=2       #number of response
m=2       #number of treatment1
n=3       #number of treatment2

aa=rbind(a1,a2,a3,a4,a5,a6,a7,a8,a9,a10)
bb=rbind(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10)

#============================ 將資料轉成矩陣 ==============================

A=matrix(array(dim=N*k*m),nrow=N,ncol=k*m)
for (i in 1:m){
A[,i]=aa[,(2*i-1)]
A[,m+i]=aa[,2*i]
}

B=matrix(array(dim=N*k*n),nrow=N,ncol=k*n)
for (i in 1:n){
B[,i]=bb[,(2*i-1)]
B[,n+i]=bb[,2*i]
}


#============================ F值 =========================================

F=function(a,b){
x=matrix(a,nrow=m,ncol=k)   #Treament1

xbar=matrix(array(dim=k),nrow=1,ncol=k)
for (i in 1:k){
xbar[1,i]=sum(x[,i])/m
}

x
xbar

Cx=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
cx=function(i){((x[i,u]-xbar[1,u])*(x[i,v]-xbar[1,v]))/(m-1)}
i=1:m
Cx[u,v]=sum(cx(i))
}
}
Cx          #matrix of covariance on Treatment1



y=matrix(b,nrow=n,ncol=k)   #Treament2

ybar=matrix(array(dim=k),nrow=1,ncol=k)
for (i in 1:k){
ybar[1,i]=sum(y[,i])/n
}
y
ybar

Cy=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
cy=function(i){((y[i,u]-ybar[1,u])*(y[i,v]-ybar[1,v]))/(n-1)}
i=1:n
Cy[u,v]=sum(cy(i))
}
}
Cy          #matrix of covariance on Treatment2

C=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
C[u,v]=((m-1)*Cx[u,v]+(n-1)*Cy[u,v])/(m+n-2)
}
}
C          # k*k matrix of pooled covariance


library(car)
inv(C)    #solve(C)
T2=((m*n)/(m+n))*((xbar-ybar)%*%inv(C)%*%t(xbar-ybar))    #T2-value
f=((m+n-k-1)/((m+n-2)*k))*T2     #F-value
f
}


F.value=array(dim=c(1,N))
for (j in 1:N){
F.value[j]=F(A[j,],B[j,])
}
F.value

for (i in 1:N){
F.p=(1/N)
if (F.value[i]>=F.value[1]) {F.p=F.p+1/N}
}
F.p

#============================ T2值 =========================================

T2=function(a,b){
x=matrix(a,nrow=m,ncol=k)   #Treament1

x1sum=0
x2sum=0
for (i in 1:m){
x1sum=x1sum+x[i,1]
x2sum=x2sum+x[i,2]
}
xbar=matrix(c(x1sum/m,x2sum/m),nrow=1,ncol=2)
x
xbar


Cx=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
cx=function(i){((x[i,u]-xbar[1,u])*(x[i,v]-xbar[1,v]))/(m-1)}
i=1:m
Cx[u,v]=sum(cx(i))
}
}
Cx          #matrix of covariance on Treatment1



y=matrix(b,nrow=n,ncol=k)   #Treament2

y1sum=0
y2sum=0
for (i in 1:n){
y1sum=y1sum+y[i,1]
y2sum=y2sum+y[i,2]
}
ybar=matrix(c(y1sum/n,y2sum/n),nrow=1,ncol=2)
y
ybar

Cy=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
cy=function(i){((y[i,u]-ybar[1,u])*(y[i,v]-ybar[1,v]))/(n-1)}
i=1:n
Cy[u,v]=sum(cy(i))
}
}
Cy          #matrix of covariance on Treatment2

C=matrix(array(dim=k*k),nrow=k,ncol=k)
for (u in 1:k){
for (v in 1:k){
C[u,v]=((m-1)*Cx[u,v]+(n-1)*Cy[u,v])/(m+n-2)
}
}
C          # k*k matrix of pooled covariance

library(car)
inv(C)    #solve(C)
T2=((m*n)/(m+n))*((xbar-ybar)%*%inv(C)%*%t(xbar-ybar))    #T2-value
T2
}

T2.value=array(dim=N)
for (i in 1:N){
T2.value[i]=T2(A[i,],B[i,])
}
T2.value


#============================= t1 ====================================
t1=function(a,b){
x=matrix(a,nrow=m,ncol=k)   #Treament1

x1sum=0
for (i in 1:m){
x1sum=x1sum+x[i,1]
}
x1bar=x1sum/m
x
x1bar

x1var=var(x[,1])


y=matrix(b,nrow=n,ncol=k)   #Treament2

y1sum=0
for (i in 1:n){
y1sum=y1sum+y[i,1]
}
y1bar=y1sum/n
y
y1bar

y1var=var(y[,1])

t1.sp2=((m-1)*x1var+(n-1)*y1var)/((m-1)+(n-1))

t1=(x1bar-y1bar)/sqrt(t1.sp2*(1/m+1/n))

t1
}

t1.value=array(dim=N)
for (i in 1:N){
t1.value[i]=t1(A[i,],B[i,])
}
t1.value


#============================= t2 ====================================

t2=function(a,b){
x=matrix(a,nrow=m,ncol=k)   #Treament1

x2sum=0
for (i in 1:m){
x2sum=x2sum+x[i,2]
}
x2bar=x2sum/m
x
x2bar

x2var=var(x[,2])


y=matrix(b,nrow=n,ncol=k)   #Treament2

y2sum=0
for (i in 1:n){
y2sum=y2sum+y[i,2]
}
y2bar=y2sum/n
y
y2bar

y2var=var(y[,2])

t2.sp2=((m-1)*x2var+(n-1)*y2var)/((m-1)+(n-1))

t2=(x2bar-y2bar)/sqrt(t2.sp2*(1/m+1/n))

t2
}

t2.value=array(dim=N)
for (i in 1:N){
t2.value[i]=t2(A[i,],B[i,])
}
t2.value


#======================= tmaxabs, tmax ==============================

t=matrix(c(t1.value,t2.value),nrow=10,ncol=k)
colnames(t)=c("t1","t2")
t


tmaxabs.value=array(dim=N)
for (i in 1:N){
tmaxabs.value[i]=max(c(abs(t[i,])))
}
tmaxabs.value

for (i in 1:N){
tmaxabs.p=(1/N)
if (tmaxabs.value[i]>=tmaxabs.value[1]) {tmaxabs.p=tmaxabs.p+1/N}
}
tmaxabs.p

tmax.value=array(dim=N)
for (i in 1:N){
tmax.value[i]=max(t[i,])
}
tmax.value

for (i in 1:N){
tmax.p=(1/N)
if (tmax.value[i]>=tmax.value[1]) {tmax.p=tmax.p+1/N}
}
tmax.p

#================  將data T2 F t1 t2 tmaxabs tmax 匯整成表格  ====================

#=====================  treat1(A) treat2(B)  =========================


AB=cbind(aa,bb)
colnames(AB)=c(rep("treat1",m*k),rep("treat2",n*k))
rownames(AB)=c(1:N)


#============================   All   =================================



All=cbind(AB,t(rbind(T2.value,F.value,t1.value,t2.value,tmaxabs.value,tmax.value)))
rownames(All)=c(1:N)
All

library(xtable)
print(xtable(All),type="html")

p.value=matrix(c(F.p,tmaxabs.p,tmax.p),nrow=1,ncol=3)
rownames(p.value)=c("p-value")
colnames(p.value)=c("F","tmaxabs","tmax")

print(xtable(p.value),type="html")



<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">

<html xmlns="http://www.w3.org/1999/xhtml">
<head><title>

</title></head>
<body>
    <form method="post" action="./OpDlFile.aspx?id=69692C979713D698" id="form1">
<div class="aspNetHidden">
<input type="hidden" name="__VIEWSTATE" id="__VIEWSTATE" value="/wEPDwUKMTc5NTg3NTg3NGRkd5gLStbXUfnHwVjiIUsXvedyIWQ=" />
</div>

<div class="aspNetHidden">

	<input type="hidden" name="__VIEWSTATEGENERATOR" id="__VIEWSTATEGENERATOR" value="E17CBA4F" />
</div>
    <h1 id="sysmsg"></h1>
    </form>
</body>
</html>