Ｒ関数
#############R 言語関数の更新###########
Ｒ言語の関数（2007年 1月 30日更新）
修正されたプログラム(manova2,Tukey.test)の修正箇所には、目印として#20070130が記入されています。
Ｒ言語の巻数（2006年12月17日更新）
関数lqda を修正(修正箇所には、目印として#20061217を記入)
Ｒ言語の関数（2006年 6月24日更新）
修正されたプログラム(Tukey.test)の修正箇所には、目印として#20060624が記入されています。
Ｒ言語の関数（2006年 6月 1日更新）
修正されたプログラム(Tukey.test)の修正箇所には、目印として#20060601が記入されています。
#############R 関数####################
コピーする際は、#start#####からend#############までのみをコピーする
こと。それより外側までコピーするとプログラム実行の際にエラーとなる
ことがあります。
###########cmdscale2##########################################################
cmdscale2<-function(d,k){
library(stats)
out<-cmdscale(d,k,eig=T)
evalue<-out$eig
points2<-out$points
for(j in 1:k) if(points2[1,j]<=0) points2[,j]<--points2[,j]
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.02)
axis1<-range(points2[,1:2])
axis2<-axis1
plot(axis1,axis2,type="n")
text(points2[,1:2],labels=c(1:nrow(d))) #corrected
write.table(c("points"),file="cmdscale2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(points2,3),file="cmdscale2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(evalue=round(evalue,3),points=round(points2,3))
}
##########fscore#########################################################
fscore<-function(x,a){
n<-ncol(x)
m<-ncol(a)
n1<-nrow(x)
n2<-ncol(x)
sd<-sqrt(apply(x,2,var)*(n1-1)/n1)
z<-scale(x,T,F)/matrix(rep(sd,n1),ncol=n2,byrow=T)
a2<-t(a)%*%a
f<-z%*%a%*%solve(a2)
par(mfrow=c(2,2),pty="s",las=1,lwd=2,tck=0.02)
axis1<-range(f[,1])
axis2<-range(f[,2])
plot(axis1,axis2,type="n")
text(a[,1],a[,2])
text(f[,1],f[,2],font=2)
write.table(c("fscore"),file="fscore-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(f,3),file="fscore-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(rawdata=round(x,3),loading=round(a,3),fscore=round(f,3))
}
###########hclust2########################################################
hclust2<-function(d){
library(stats)
d<-as.dist(d)
out<-hclust(d,method="average")
plclust(out)
}
###########isoMDS2##########################################################
isoMDS2<-function(d,y,k){
library(MASS)
out<-isoMDS(d,y,k)
points2<-out$points
stress<-out$stress
for(j in 1:k) if(points2[1,j]<=0) points2[,j]<--points2[,j]
par(mfrow=c(2,2),pty="s",las=1,lwd=2,tck=0.02)
axis1<-range(points2[,1:2])
axis2<-axis1
plot(axis1,axis2,type="n")
text(points2[,1:2])
dst.obtained<-c(0,max(c(dist(points2),as.dist(d))))
dst.given<-dst.obtained
plot(dst.obtained,dst.given,type="n")
points(dist(points2),as.dist(d),pch="o")
write.table(c("points"),file="isomds2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(points2,3),file="isomds2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(points=round(points2,3),stress=round(stress,3))
}
##########kmeans#######################################################
kmeans2<-function(x,center){
library(stats)
iter<-20
n<-nrow(x)
cl <- kmeans(x, center, iter)
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.02)
axis1<-range(c(x[,1:2],cl$center[,1:2]))
axis2<-axis1
plot(axis1,axis2,type="n")
text(x[,1:2])
points(cl$centers[,1:2])
cluster<-rbind(1:n,cl$cluster)
list(center=cl$center,result=cluster)
}
##########hqt1###########################################################
hqt1<-function(x,y){
a<-solve(t(x)%*%x)%*%t(x)%*%y
list(a=round(a,3))
}
##########hqt2###########################################################
hqt2<-function(n,d){
n1<-n[1]
n2<-n[2]
d1<-d[1:n1,]
n3<-n1+1
d2<-d[n3:nrow(d),]
y<-apply(d,2,"sum")
y<-matrix(y,nrow=1,byrow=T)
x1<-apply(d1,2,"sum")
x2<-apply(d2,2,"sum")
x1<-matrix(x1,nrow=1,byrow=T)
x2<-matrix(x2,nrow=1,byrow=T)
t1<-t(d)%*%d-t(y)%*%y/nrow(d)
b<-t(x1)%*%x1/n1+t(x2)%*%x2/n2-t(y)%*%y/nrow(d)
evalue1<-eigen(t1)$values
evector1<-eigen(t1)$vectors
p<-evector1%*%sqrt(diag(evalue1))
a<-solve(p)%*%b%*%solve(t(p))
evalue<-eigen(a)$values
evector2<-eigen(a)$vectors
evector3<-solve(t(p))%*%evector2
xx<-sqrt(diag(t(evector3)%*%evector3))
evector<-evector3/xx
Y<-d%*%evector[,1]
write.table(c("discriminant score"),file="hqt2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(Y,3),file="hqt2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(d=d,evalue=round(evalue[1],3),coef=round(evector[,1],3),score=round(Y,3))
}
##########hqt3##########################################################
hqt3<-function(x){
n<-ncol(x)
d1<-apply(x,1,sum)
d2<-apply(x,2,sum)
d3<-ifelse(d1==0,0,1/d1)
d4<-ifelse(d2==0,0,1/sqrt(d2))
d3<-diag(d3,nrow=length(d3))
d4<-diag(d4,nrow=length(d4))
m1<-t(x)
y<-d4%*%m1%*%d3%*%x%*%d4
evalue<-eigen(y)$values
evector<-eigen(y)$vectors
for ( j in 1:ncol(evector)) if(evector[1,j]<=0) evector[,j]<--evector[,j]
evector<-evector*sqrt(sum(x))
b<-d4%*%evector
evalue2<-1/sqrt(evalue)
evalue2<-diag(evalue2,ncol=length(evalue2))
a<-d3%*%x%*%b%*%evalue2
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.02)
axis2<-range(c(a[,2],a[,3],b[,2],b[,3]))
axis3<-axis2
plot(axis2,axis3,type="n")
text(a[,2],a[,3],font=1)
text(b[,2],b[,3],font=2)
write.table(c("row"),file="hqt3-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(a,3),file="hqt3-table.xls",
col.names=F,append=T,quote=F,sep="\t")
write.table(c("col"),file="hqt3-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(b,3),file="hqt3-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(data=x,evalue=round(evalue,3),row=round(a,3),col=round(b,3))
}
##########hqt4###########################################################
hqt4<-function(x){
nr<-nrow(x)
nc<-ncol(x)
y<-x
e<-matrix(0,nrow=nc,ncol=nc)
x[x==0]<--1
for( i in 1:nc) for ( j in 1:nc) e[i,j]<-sum(abs(x[,i]+x[,j]))/2
e[row(e)==col(e)]<-0
d<-apply(e,1,"sum")
e2<-e-diag(d,nrow=length(d))
evalue<-round(eigen(e2)$values,3)
evector<-round(eigen(e2)$vectors,3)
for ( j in 1:ncol(evector)) if(evector[1,j]<=0) evector[,j]<--evector[,j]
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.03)
axis2<-c(-1,1)
axis3<-axis2
plot(axis2,axis3,type="n")
text(evector[,2],evector[,3])
write.table(c("points"),file="hqt4-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(evector,3),file="hqt4-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(data=y,eij=e2,evalue=evalue,points=evector)
}
##########lqda###########################################################
lqda<-function(y1,y2){
np<-ncol(y1)
ny1<-nrow(y1)
ny2<-nrow(y2)
n<-ny1+ny2
cov1<-var(y1)
cov2<-var(y2)
cov12<-((ny1-1)*cov1+(ny2-1)*cov2)/(ny1+ny2-2)
det1<-det(cov1)
det2<-det(cov2)
det12<-det(cov12)
v<-det1^(ny1/2)*det2^(ny2/2)/det12^(n/2)
chisqvalue<--2*log(v)
df<-np*(np+1)/2
pvalue<-1-pchisq(chisqvalue,df)
n1<-np*ny1
n2<-np*ny1+1
cmeany1<-apply(y1,2,mean)
cmeany2<-apply(y2,2,mean)
dmean<-cmeany1-cmeany2
if(pvalue<0.025) result<-c("different cov") else result<-c("same cov")
if(pvalue<0.025) select<-2 else select<-1
print(cbind(result),quote=F)
print(cbind(chisqvalue,df,pvalue))
if(select==1){
zy1<-y1-matrix(rep(cmeany1,ny1),ncol=np,byrow=T)
zy2<-y2-matrix(rep(cmeany2,ny2),ncol=np,byrow=T)
sij<-(t(zy1)%*%zy1+t(zy2)%*%zy2)/(ny1+ny2-2)
a<-solve(sij)%*%dmean
midmean<-(cmeany1+cmeany2)/2
a0<- -t(a)%*%midmean
z1<-y1%*%a+rep(a0,ny1)
z2<-y2%*%a+rep(a0,ny2)
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.02)
x1<-range(y1[,1],y2[,1])
x2<-range(y1[,2],y2[,2])
plot(x1,x2,type="n")
text(y1[,1],y1[,2],font=1)
text(y2[,1],y2[,2],font=2)
abline(-a0/a[2],-a[1]/a[2])
}
if(select==2){
ny3<-ny1+ny2
zy11<-y1-matrix(rep(cmeany1,ny1),ncol=np,byrow=T)
zy12<-y2-matrix(rep(cmeany1,ny2),ncol=np,byrow=T)
zy21<-y1-matrix(rep(cmeany2,ny1),ncol=np,byrow=T)
zy22<-y2-matrix(rep(cmeany2,ny2),ncol=np,byrow=T)
cov1<-var(y1)
cov2<-var(y2)
d11<-diag((zy11)%*%solve(cov1)%*%t(zy11))
d12<-diag((zy12)%*%solve(cov1)%*%t(zy12))
d21<-diag((zy21)%*%solve(cov2)%*%t(zy21))
d22<-diag((zy22)%*%solve(cov2)%*%t(zy22))
d1<-c(d11,d12)
d2<-c(d21,d22)
dd<-d2-d1
dd2<-dd
dd2[dd2>=0]<-1
dd2[dd2<0]<-2
}
write.table(c("discriminant score1"),file="lqda-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(z1,3),file="lqda-table.xls",
col.names=F,append=T,quote=F,sep="\t")
write.table(c("discriminant score2"),file="lqda-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(z2,3),file="lqda-table.xls",
col.names=F,append=T,quote=F,sep="\t")
if(select==1){
list(group1=y1,group2=y2,cov=round(sij,3),coef=round(c(a0,a),3),
score1=round(z1,3),score2=round(z2,3))
}
else {
list(group1=y1,group2=y2,cov1=round(cov1,3),cov2=round(cov2,3),
dist1=round(d1,3),dist2=round(d2,3),result1=round(dd,3),result2=dd2)
}
}
##########lsfit2##########################################################
lsfit2<-function(x,y)
{
if(is.vector(x)){
n<-length(x)
meanx<-mean(x)
meany<-mean(y)
sx<-sqrt(var(x)*(n-1)/n)
sy<-sqrt(var(y)*(n-1)/n)
sxy<-mean(x*y)-meanx*meany
a1<-sxy/sx^2
a0<-meany-a1*meanx
yh<-a0+a1*x
R<-cor(y,yh)
RR<-R^2
us<-sqrt(sum((y-yh)^2)/(n-2))
v0<-sum(x^2)/(n*sum((x-meanx)^2))
t0<-a0/(us*sqrt(v0))
v1<-1/sum((x-meanx)^2)
t10<-a1/(us*sqrt(v1))
t11<-abs(a1-1)/(us*sqrt(v1))
tr0<-(R*sqrt(n-2))/sqrt(1-R^2)
df<-n-2
pvalue0<-1-pt(abs(t0),df)
pvalue10<-1-pt(abs(t10),df)
pvalue11<-1-pt(abs(t11),df)
pvaluer<-1-pt(abs(tr0),df)
}
if(is.matrix(x)){
nr<-nrow(x)
nc<-ncol(x)
meanx<-apply(x,2,mean)
meany<-mean(y)
var1<-var(cbind(x,y))*(nr-1)/nr
var2<-var1[1:nc,1:nc]
varx<-diag(var2)
covxy<-var1[1:nc,nc+1]
invv<-solve(var2)
a<-invv%*%covxy
a0<-meany-sum(a*meanx)
yh<-a0+x%*%a
R<-cor(y,yh)
Rh<-sqrt(((nr-1)*R^2-nc)/(nr-nc-1))
e1<-y-yh
df<-nr-nc-1
ve<-sum((yh-y)^2)/df
df1<-nc
df2<-nr-nc-1
fvalue<-R^2*df2/((1-R^2)*df1)
if( fvalue >= 1) pvalue2<-1-pf(fvalue,df1,df2)
if( fvalue < 1) pvalue2<-pf(fvalue,df1,df2)
tvalue<-a/sqrt(diag(invv)*ve/nr)
pvalue<-1-pt(abs(tvalue),df)
sumxx<-0
for ( i in 1:nc) for ( j in 1:nc) sumxx<-sumxx+meanx[i]*meanx[j]*invv[i,j]
tvalue0<-a0/sqrt((1+sumxx)*ve/nr)
pvalue0<-1-pt(abs(tvalue0),df)
x1<-matrix(0,nrow=nr,ncol=nc-1)
ruv<-0
tvaluepcor<-0
pvaluepcor<-0
for ( k in 1:nc){
for(j in 1:nc) {
if( j < k) x1[,j]<-x[,j]
if( j > k) x1[,j-1]<-x[,j]
}
xk<-x[,k]
meanx2<-apply(x1,2,mean)
meanxk<-mean(xk)
vary<-var(y)*(nr-1)/nr
varxk<-var(xk)*(nr-1)/nr
var1<-var(cbind(x1,xk))*(nr-1)/nr
var22<-var1[1:nc-1,1:nc-1]
varx<-diag(var22)
covxy2<-var1[1:nc-1,nc]
invv<-solve(var22)
ak<-invv%*%covxy2
ak0<-meanxk-sum(ak*meanx2)
xkh<-ak0+x1%*%ak
out<-lsfit(x1,y)$coef
nc2<-ncol(x1)+1
yh2<-out[1]+x1%*%out[2:nc2]
v<-y-yh2
u<-xk-xkh
ruv[k]<-cor(u,v)
tvaluepcor[k]<-ruv[k]*sqrt(nr-3)/sqrt(1-ruv[k]^2)
dfpcor<-nr-3
pvaluepcor[k]<-1-pt(abs(tvaluepcor[k]),dfpcor)
}
}
if(is.matrix(x)) list(x=t(x),y=y,yh=round(t(yh),3),meanx=round(meanx,3),
meany=round(meany,3),coef=round(c(a0,a),3),
R=round(R,3),Rh=round(Rh,3),pvaluer=round(pvalue2,3),
tvaluea=round(c(tvalue0,tvalue),3),
pvaluea=round(c(pvalue0,pvalue),3),df=df,pcor=round(ruv,3),
tvaluepcor=round(tvaluepcor,3),dfpcor=dfpcor,pvaluepcor=round(pvaluepcor,3))
else list(x=x,y=y,yh=round(yh,2),
meanx=round(meanx,2),meany=round(meany,2),sdx=round(sx,2),sdy=round(sy,2),
sxy=round(sxy,2),a=round(cbind(a0,a1),2),R=round(R,3),RR=round(RR,3),
tvalue0=round(cbind(t0,df,pvalue0),3),tvalue10=round(cbind(t10,df,pvalue10),3),
tvalue11=round(cbind(t11,df,pvalue11),3),tvaluer=round(cbind(tr0,df,pvaluer),3))
}
##########manova2###################################################
manova2<-function(data,ma,mb,mc,n,rep){
na<-ma;nb<-mb;nc<-mc
value<--1
x0<-array(data,c(n,nc,nb,na));x<-aperm(x0,c(4,3,2,1));
x01<-array(0,c(na,nb,nc,n));x01[x!=value]<-1
x[x==value]<-0;n2<-ncol(data);n3<-nc*nb*na;n4<-nc*n;n5<-n3*n;select<-1
if(n2==n3){
nijk<-apply(x01,c(1,2,3),sum);nijk2<-apply(x01,c(1,2,4),sum)
nikk2<-apply(x01,c(1,3,4),sum);njkk2<-apply(x01,c(2,3,4),sum)
ni<-apply(nijk,1,sum);nj<-apply(nijk,2,sum)
nk<-apply(nijk,3,sum);nij<-apply(nijk,c(1,2),sum)
nik<-apply(nijk,c(1,3),sum);njk<-apply(nijk,c(2,3),sum)
na2<-nijk[,1,1];nab<-nijk[,,1];nbc<-nijk[1,,];nac<-nijk[,1,]
hmean<-na*nb*nc/(sum(1/nijk));gmean<-sum(x)/sum(x01)
meanxijk<-apply(x,c(1,2,3),sum)/nijk
meanxijk2<-apply(x,c(1,2,4),sum)/nijk2
meanxikk2<-apply(x,c(1,3,4),sum)/nikk2
meanxjkk2<-apply(x,c(2,3,4),sum)/njkk2
meanxij<-apply(meanxijk,c(1,2),mean)
meanxik<-apply(meanxijk,c(1,3),mean)
meanxik2<-apply(meanxijk2,c(1,3),mean)
meanxjk<-apply(meanxijk,c(2,3),mean)
meanxjk2<-apply(meanxijk2,c(2,3),mean)
meanxkk2<-apply(meanxikk2,c(2,3),mean)
meanxi<-apply(meanxij,1,mean)
meanxj<-apply(meanxij,2,mean)
meanxk<-apply(meanxjk,2,mean)
meanxk2<-apply(meanxik2,2,mean)
SSbetween<-0;SSwithin<-0;SSa<-0
SSb<-0;SSc<-0;SSab<-0;SSac<-0;SSbc<-0
SSabc<-0;SSe<-0;SSe.a<-0
SSe.b<-0;SSe.c<-0
SSe.ab<-0;SSe.ac<-0
SSe.bc<-0;SSe.abc<-0;SSe.w<-0;SSt<-0
abcs<-sum(x^2)
abs1<-sum(apply(x,c(1,2,4),sum)^2/(nc))
a<-sum(apply(x,1,sum)^2/(n*nb*nc))
b<-sum(apply(x,2,sum)^2/(n*na*nc))
c1<-sum(apply(x,3,sum)^2/(n*na*nb))
ab<-sum(apply(x,c(1,2),sum)^2/(nab*nc))
bc<-sum(apply(x,c(2,3),sum)^2/(nbc*na))
ac<-sum(apply(x,c(1,3),sum)^2/(nac*nb) )
abc<-sum(apply(x,c(1,2,3),sum)^2/nijk)
xpp<-sum(x)^2/(n*na*nb*nc)
abcp<-sum(meanxijk^2)
g<-sum(x)
gp<-sum(meanxijk)
xp<-gp^2/(na*nb*nc)
ap<-nb*nc*sum(meanxi^2)
bp<-na*nc*sum(meanxj^2)
abp<-nc*sum(meanxij^2)
cp<-na*nb*sum(meanxk^2)
acp<-nb*sum(meanxik^2)
bcp<-na*sum(meanxjk^2)
if(rep>=2)
{
a<-sum(apply(x,1,sum)^2/(na2*nb*nc))
as<-sum(apply(x,c(1,4),sum)^2/(nb*nc))
ab<-sum(apply(x,c(1,2),sum)^2/(na2*nc))
abs1<-sum(apply(x,c(1,2,4),sum)^2/(nc))
ac<-sum(apply(x,c(1,3),sum)^2/(na2*nb))
acs<-sum(apply(x,c(1,3,4),sum)^2/(nb))
}
if(rep==3)
{
s<-sum(apply(x,4,sum)^2/(na*nb*nc))
bs<-sum(apply(x,c(2,4),sum)^2/(na*nc))
cs<-sum(apply(x,c(3,4),sum)^2/(na*nb))
bcs<-sum(apply(x,c(2,3,4),sum)^2/na)
}
SSa<-hmean*(ap-xp)
SSb<-hmean*(bp-xp)
SSc<-hmean*(cp-xp)
SSab<-hmean*(abp-ap-bp+xp)
SSac<-hmean*(acp-ap-cp+xp)
SSbc<-hmean*(bcp-bp-cp+xp)
SSabc<-hmean*(abcp-abp-acp-bcp+ap+bp+cp-xp)
SSe<-abcs-abc
for(i in 1:na)for(j in 1:nb)for(k in 1:nc)for(k2 in 1:n){
if(x01[i,j,k,k2]!=0) SSt<-SSt+(x[i,j,k,k2]-gmean)^2
}
if(rep==1){
SSe.between<-abs1-ab
SSe.within<-abcs-abc-abs1+ab
}
if(rep==2){
SSe.a<-as-a
SSe.between<-SSe.a
SSe.b<-abs1-ab-as+a
SSe.c<-acs-ac-as+a
SSe.bc<-abcs-abc-abs1-acs+ab+ac+as-a
}
if(rep==3){
SSe.between<-s-xpp
SSe.a<-as-a-s+xpp
SSe.b<-bs-b-s+xpp
SSe.c<-cs-c1-s+xpp
SSe.ab<-abs1-ab-as-bs+a+b+s-xpp
SSe.ac<-acs-ac-as-cs+a+c1+s-xpp
SSe.bc<-bcs-bc-bs-cs+b+c1+s-xpp
SSe.abc<-abcs-abc-abs1-acs-bcs+ab+ac+bc+as+bs+cs-a-b-c1-s+xpp
}
if(na!=1) dfa<-na-1
if(rep==2) {dfe.between<-sum(na2)-na-1
dfe.within<-(sum(na2)-na)*(nb*nc-1)
}
if(rep==3) {
dfe.between<-n-1;dfe.within<-n*(na*nb*nc-1)
}
dfa<-na-1;dfb<-nb-1;dfc<-nc-1;dfab<-(na-1)*(nb-1)
dfac<-(na-1)*(nc-1);dfbc<-(nb-1)*(nc-1)
dfabc<-(na-1)*(nb-1)*(nc-1)
if(rep==0) dfe<-sum(nijk)-na*nb*nc #dfe<-na*nb*nc*(n-1)
if(rep==1) {dfe.within<-(sum(nab)-na*nb)*(nc-1)
dfe.between<-sum(nab)-na*nb
}
if(rep==2) {dfe.a<-sum(na2)-na
dfe.b<-(sum(na2)-na)*(nb-1)
dfe.c<-(sum(na2)-na)*(nc-1)
dfe.bc<-(sum(na2)-na)*(nb-1)*(nc-1)
}
if(rep==3) {dfe.a<-(na-1)*(n-1)
dfe.b<-(nb-1)*(n-1);
dfe.c<-(nc-1)*(n-1)
dfe.ab<-(na-1)*(nb-1)*(n-1)
dfe.ac<-(na-1)*(nc-1)*(n-1)
dfe.bc<-(nb-1)*(nc-1)*(n-1)
dfe.abc<-(na-1)*(nb-1)*(nc-1)*(n-1)
}
dft<-sum(nijk)-1
if(rep==1) dft<-sum(nab)*nc-1
if(rep==2) dft<-sum(na2)*nb*nc-1
MSa<-SSa/dfa;MSb<-SSb/dfb;MSc<-SSc/dfc
MSab<-SSab/dfab;MSac<-SSac/dfac
MSbc<-SSbc/dfbc;MSabc<-SSabc/dfabc;
if(rep==0) MSe<-SSe/dfe
if(rep!=0) MSe.between<-SSe.between/dfe.between
if(rep==1) MSe.within<-SSe.within/dfe.within
if(rep>1) {MSe.a<-SSe.a/dfe.a;MSe.b<-SSe.b/dfe.b
MSe.c<-SSe.c/dfe.c;MSe.bc<-SSe.bc/dfe.bc
}
if(rep==3) {MSe.ab<-SSe.ab/dfe.ab;MSe.ac<-SSe.ac/dfe.ac;
MSe.abc<-SSe.abc/dfe.abc
}
if(rep==0){Fa<-MSa/MSe;Fab<-MSab/MSe;Fac<-MSac/MSe
Fabc<-MSabc/MSe;Fb<-MSb/MSe;Fc<-MSc/MSe;Fbc<-MSbc/MSe
}
if(rep>1) Fa<-MSa/MSe.a
if(rep==1) {Fa<-MSa/MSe.between;Fb<-MSb/MSe.between
Fc<-MSc/MSe.within;Fab<-MSab/MSe.between
Fac<-MSac/MSe.within;Fbc<-MSbc/MSe.within;Fabc<-MSabc/MSe.within
}
if(rep>1) {Fb<-MSb/MSe.b;Fc<-MSc/MSe.c;Fbc<-MSbc/MSe.bc
}
if(rep==2) {Fab<-MSab/MSe.b;Fac<-MSac/MSe.c;Fabc<-MSabc/MSe.bc
}
if(rep==3) {Fab<-MSab/MSe.ab;Fac<-MSac/MSe.ac;Fabc<-MSabc/MSe.abc
}
if(rep==0){
if(na !=1) {Pa<-1-pf(Fa,dfa,dfe);Pab<-1-pf(Fab,dfab,dfe)
Pac<-1-pf(Fac,dfac,dfe);Pabc<-1-pf(Fabc,dfabc,dfe)
}
if(nb !=1) {Pb<-1-pf(Fb,dfb,dfe);Pbc<-1-pf(Fbc,dfbc,dfe)
}
Pc<-1-pf(Fc,dfc,dfe)
}
if(rep==1 && na!=1) Pa<-1-pf(Fa,dfa,dfe.between)
if(rep>1 && na!=1) {Pa<-1-pf(Fa,dfa,dfe.a);Pb<-1-pf(Fb,dfb,dfe.b)
Pc<-1-pf(Fc,dfc,dfe.c);Pbc<-1-pf(Fbc,dfbc,dfe.bc)
}
if(rep>1 && na==1 && nb!=1) {Pb<-1-pf(Fb,dfb,dfe.b)
Pc<-1-pf(Fc,dfc,dfe.c);Pbc<-1-pf(Fbc,dfbc,dfe.bc)
}
# if(rep==1 && nb!=1) {Pb<-1-pf(Fb,dfe.between,dfe.within)
if(rep==1 && nb!=1) {Pb<-1-pf(Fb,dfb,dfe.between) #20070130
Pc<-1-pf(Fc,dfc,dfe.within);Pbc<-1-pf(Fbc,dfbc,dfe.within)
if(na!=1) {Pabc<-1-pf(Fabc,dfabc,dfe.within);
Pab<-1-pf(Fab,dfab,dfe.between);Pac<-1-pf(Fac,dfac,dfe.within)
}
}
if(rep==1 && nb==1) Pc<-1-pf(Fc,dfc,dfe.within)
if(rep==2 && na!=1) {Pab<-1-pf(Fab,dfab,dfe.c)
Pac<-1-pf(Fac,dfac,dfe.c)
Pabc<-1-pf(Fabc,dfabc,dfe.bc)
}
if(rep==3) {Pab<-1-pf(Fab,dfab,dfe.ab);Pac<-1-pf(Fac,dfac,dfe.ac)
Pabc<-1-pf(Fabc,dfabc,dfe.abc)
}
print("data");print(data);
if(nb!=1 && na!=1) {print("meanxijk");print(meanxijk)
print("meanxij");print(meanxij)
print("meanxik");print(meanxik);print("meanxjk")
print(meanxjk);print("meanxi");print(meanxi)
print("meanxj");print(meanxj)
}
if(nb!=1 && na==1) {print("meanxj");print(meanxij)
print("meanxk");print(meanxik);print("meanxjk")
print(meanxjk)
}
if(nb==1 && na==1) {print("meanxk");print(meanxik)
}
if(nb==1 && na==1 && n4!=sum(nijk)){
select<-2
y<-sum(x)^2/sum(x01)
as2<-sum(x^2)
a2<-sum(apply(x[1,,,],1,sum)^2/nijk)
SSc<-a2-y
SSe<-as2-a2
SSe.c<-SSe
dfc<-nc-1
dfe<-sum(x01)-nc
dfe.c<-dfe
MSc<-SSc/dfc
MSe<-SSe/dfe
MSe.c<-MSe
Fc<-MSc/MSe
Pc<-1-pf(Fc,dfc,dfe)
}
print("gmean");print(gmean)
write.table(cbind("Source","SS","df","MS","F","P"),
file="aov-table.xls",row.names=F,col.names=F,
append=F,quote=F,sep="\t")
if(na!=1) {print(cbind(SSa,dfa,MSa,Fa,Pa))
write.table(cbind(SSa,dfa,MSa,Fa,Pa),file="aov-table.xls",
row.names=c("A"),col.names=F,append=T,quote=F,sep="\t")
}
if(na!=1 && nb!=1) {print(cbind(SSb,dfb,MSb,Fb,Pb))
write.table(cbind(SSb,dfb,MSb,Fb,Pb),file="aov-table.xls",
row.names=c("B"),col.names=F,append=T,quote=F,sep="\t")
}
if(na==1 && nb!=1) {print(cbind(SSb,dfb,MSb,Fb,Pb))
write.table(cbind(SSb,dfb,MSb,Fb,Pb),file="aov-table.xls",
row.names=c("B"),col.names=F,append=T,quote=F,sep="\t")
}
if(nb!=1) {print(cbind(SSc,dfc,MSc,Fc,Pc))
write.table(cbind(SSc,dfc,MSc,Fc,Pc),file="aov-table.xls",
row.names=c("C"),col.names=F,append=T,quote=F,sep="\t")
}
if(nb==1) {print(cbind(SSc,dfc,MSc,Fc,Pc))
write.table(cbind(SSc,dfc,MSc,Fc,Pc),file="aov-table.xls",
row.names=c("C"),col.names=F,append=T,quote=F,sep="\t")
}
if(na!=1) {print(cbind(SSab,dfab,MSab,Fab,Pab))
print(cbind(SSac,dfac,MSac,Fac,Pac));
write.table(cbind(SSab,dfab,MSab,Fab,Pab),file="aov-table.xls",
row.names=c("AB"),col.names=F,append=T,quote=F,sep="\t");
write.table(cbind(SSac,dfac,MSac,Fac,Pac),file="aov-table.xls",
row.names=c("AC"),col.names=F,append=T,quote=F,sep="\t")
}
if(nb!=1) {print(cbind(SSbc,dfbc,MSbc,Fbc,Pbc))
write.table(cbind(SSbc,dfbc,MSbc,Fbc,Pbc),file="aov-table.xls",
row.names=c("BC"),col.names=F,append=T,quote=F,sep="\t")
}
if(na!=1) {print(cbind(SSabc,dfabc,MSabc,Fabc,Pabc))
write.table(cbind(SSabc,dfabc,MSabc,Fabc,Pabc),
file="aov-table.xls",row.names=c("ABC"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==0) {print(cbind(SSe,dfe,MSe))
write.table(cbind(SSe,dfe,MSe),file="aov-table.xls",
row.names=c("e.between"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==1 && na!=1 && nb!=1) {print(cbind(SSe.between,dfe.between,MSe.between))
write.table(cbind(SSe.between,dfe.between,MSe.between),file="aov-table.xls",
row.names=c("e.between"),col.names=F,append=T,quote=F,sep="\t");
SSe.c<-SSe.within;dfe.c<-dfe.within;MSe.c<-MSe.within
print(cbind(SSe.c,dfe.c,MSe.c))
write.table(cbind(SSe.c,dfe.c,MSe.c),file="aov-table.xls",
row.names=c("e.c"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==1 && na==1 && nb!=1) {SSe.b<-SSe.between;dfe.b<-dfe.between;
MSe.b<-MSe.between;
SSe.c<-SSe.within;dfe.c<-dfe.within;MSe.c<-MSe.within;
print(cbind(SSe.b,dfe.b,MSe.b));
write.table(cbind(SSe.b,dfe.b,MSe.b),file="aov-table.xls",
row.names=c("e.b"),col.names=F,append=T,quote=F,sep="\t");
print(cbind(SSe.c,dfe.c,MSe.c))
write.table(cbind(SSe.c,dfe.c,MSe.c),file="aov-table.xls",
row.names=c("e.c"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==1 && nb==1) {SSe.c<-SSe.within;dfe.c<-dfe.within;MSe.c<-MSe.within;
print(cbind(SSe.between,dfe.between))
write.table(cbind(SSe.between,dfe.between),file="aov-table.xls",
row.names=c("e.between"),col.names=F,append=T,quote=F,sep="\t")
print(cbind(SSe.c,dfe.c,MSe.c))
write.table(cbind(SSe.c,dfe.c,MSe.c),file="aov-table.xls",
row.names=c("e.c"),col.names=F,append=T,quote=F,sep="\t")
}
if(na!=1 && rep>1) {print(cbind(SSe.a,dfe.a,MSe.a))
write.table(cbind(SSe.a,dfe.a,MSe.a),file="aov-table.xls",
row.names=c("e.a"),col.names=F,append=T,quote=F,sep="\t")
}
if(nb!=1 && rep>1) {print(cbind(SSe.b,dfe.b,MSe.b))
write.table(cbind(SSe.b,dfe.b,MSe.b),file="aov-table.xls",
row.names=c("e.b"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep>1) {print(cbind(SSe.c,dfe.c,MSe.c))
write.table(cbind(SSe.c,dfe.c,MSe.c),file="aov-table.xls",
row.names=c("e.c"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==3) {print(cbind(SSe.ab,dfe.ab,MSe.ab))
write.table(cbind(SSe.ab,dfe.ab,MSe.ab),file="aov-table.xls",
row.names=c("e.ab"),col.names=F,append=T,quote=F,sep="\t");
print(cbind(SSe.ac,dfe.ac,MSe.ac))
write.table(cbind(SSe.ac,dfe.ac,MSe.ac),file="aov-table.xls",
row.names=c("e.ac"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep>1) {print(cbind(SSe.bc,dfe.bc,MSe.bc))
write.table(cbind(SSe.bc,dfe.bc,MSe.bc),file="aov-table.xls",
row.names=c("e.bc"),col.names=F,append=T,quote=F,sep="\t")
}
if(rep==3) {print(cbind(SSe.abc,dfe.abc,MSe.abc))
write.table(cbind(SSe.abc,dfe.abc,MSe.abc),file="aov-table.xls",
row.names=c("e.abc"),col.names=F,append=T,quote=F,sep="\t")
}
if(n5==sum(nijk)){print(cbind(SSt,dft))
write.table(cbind(SSt,dft),file="aov-table.xls",
row.names=c("total"),col.names=F,append=T,quote=F,sep="\t")
}
if(nb==1 && na==1 && n4!=sum(nijk)){print(cbind(SSt,dft))
write.table(cbind(SSt,dft),file="aov-table.xls",
row.names=c("total"),col.names=F,append=T,quote=F,sep="\t")
}
}
if(n2!=n3) print("input error. The column size of data should be equal to na*nb*nc")
}
##########Nonparametric-anova###################################################
nonparametric.anova<-function(data,m,n,opt,cvalue){
x<-data
if(opt==1) {
value<--1
x01<-matrix(0,nrow=n,ncol=m);x01[x!=value]<-1
value2<-max(x)+100
x[x==value]<-value2
n2<-apply(x01,2,sum)
sumn2<-sum(n2)
r<-rank(x)
if(max(n2)!=min(n2)) r[r==max(r)]<-0
r<-matrix(r,ncol=m,byrow=F)
sumr<-apply(r,2,sum)
mr<-sumr/n2
c0<-table(r)
c1<-c0[c0>1]
h<-sum(sumr^2/n2)*12/(sumn2*(sumn2+1))-3*(sumn2+1)
c2<-1-sum(c1^3-c1)/(sumn2^3-sumn2)
H<-h/c2
n<-n2
pvalue<-1-pchisq(H,m-1)
print("Kruskal-Wallis Test Result")
print("rank data")
print(cbind(r))
print(cbind(sumr,mr))
print(cbind(n))
print(cbind(H,pvalue))
if(pvalue < 0.025){
####multiple comparison(Bonferroni-Dunn method)####
alpha<-0.05
alpha2<-alpha/(m*(m-1))
zvalue<--qnorm(alpha2,0,1)
m1<-m-1
print("Multiple Comparison(Bonferroni-Dunn method)")
k<-0
cd<-0
sig005<-0
mrj<-0
mrk<-0
for( j1 in 1:m1) {
j3<-j1+1
for( j2 in j3:m){
k<-k+1
cd[k]<-zvalue*sqrt(sumn2*(sumn2+1)*(1/n2[j1]+1/n2[j2])/12)
if(abs(mr[j1]-mr[j2])>cd[k] ) sig005[k]<-1 else sig005[k]<-0
mrj[k]<-mr[j1]
mrk[k]<-mr[j2]
}
}
print(cbind(mrj,mrk,cd,sig005))
}
}
if(opt==2){
value<--1
x01<-matrix(0,nrow=n,ncol=m);x01[x!=value]<-1
n2<-apply(x01,2,sum)
sumn2<-sum(n2)
x2<-apply(x,2,sort)
print("Jonckheere-Terpstra Test Result")
print("raw data")
print(x2)
uvalue<-0
k<-0
m1<-m-1
for( j1 in 1:m1) {j3<-j1+1
for ( j2 in j3:m){
k<-k+1
sumy<-0
for ( i1 in 1:n) for( i2 in 1:n) {
if(x2[i1,j1]>0 && x2[i2,j2]>0){
if(x2[i1,j1]< x2[i2,j2]) sumy<-sumy+1
if(x2[i1,j1]==x2[i2,j2]) sumy<-sumy+0.5
}
}
uvalue[k]<-sumy
}
}
jvalue<-sum(uvalue)
print(cbind(uvalue))
if(cvalue == 0){
num<-jvalue-(sumn2^2-sum(n2^2))/4
den<-sqrt((sumn2^2*(2*sumn2+3)-sum(n2^2*(2*n2+3)))/72)
zvalue<-num/den
pvalue<-1-pnorm(zvalue,0,1)
print(cbind(jvalue,zvalue,pvalue))
}
if(cvalue !=0 ) {
if(jvalue>=cvalue) print("significant at 5% level ") else print(" not significant at 5% level")
}
if(jvalue>=cvalue || pvalue < 0.025){
value<--1
x01<-matrix(0,nrow=n,ncol=m);x01[x!=value]<-1
value2<-max(x)+100
x[x==value]<-value2
n2<-apply(x01,2,sum)
sumn2<-sum(n2)
r<-rank(x)
if(max(n2)!=min(n2)) r[r==max(r)]<-0
r<-matrix(r,ncol=m,byrow=F)
sumr<-apply(r,2,sum)
mr<-sumr/n2
####multiple comparison(Bonferroni-Dunn method)####
alpha<-0.05
alpha2<-alpha/(m*(m-1))*2
zvalue<--qnorm(alpha2,0,1)
m1<-m-1
print("Multiple Comparison(Bonferroni-Dunn method)")
k<-0
cd<-0
sig005<-0
mrj<-0
mrk<-0
for( j1 in 1:m1) {
j3<-j1+1
for( j2 in j3:m){
k<-k+1
cd[k]<-zvalue*sqrt(sumn2*(sumn2+1)*(1/n2[j1]+1/n2[j2])/12)
if(abs(mr[j1]-mr[j2])>cd[k] ) sig005[k]<-1 else sig005[k]<-0
mrj[k]<-mr[j1]
mrk[k]<-mr[j2]
}
}
print(cbind(mrj,mrk,cd,sig005))
}
}
if(opt==3){
r<-t(apply(x,1,rank))
sumr<-apply(r,2,sum)
mr<-apply(r,2,mean)
fr<-12*(sum(sumr^2))/(n*m*(m+1))-3*n*(m+1)
c0<-apply(r,1,table)
c1<-unlist(c0)
c2<-c1[c1>1]
if(length(c2) > 0) {c3<-1-sum(c2^3-c2)/(n*(m^3-m))
fr2<-fr/c3
} else fr2<-fr
pvalue<-1-pchisq(fr2,m-1)
Q<-fr2
print("Friedman Test Result")
print("rank data")
print(cbind(r))
print(cbind(sumr,mr))
print(cbind(n,Q,pvalue))
if(pvalue < 0.025) {
####multiple comparison(Bonferroni-Dunn method)####
alpha<-0.05
alpha2<-alpha/(m*(m-1))
zvalue<--qnorm(alpha2,0,1)
m1<-m-1
print("Multiple Comparison(Bonferroni-Dunn method)")
k<-0
cd<-0
sig005<-0
mrj<-0
mrk<-0
for( j1 in 1:m1) {
j3<-j1+1
for( j2 in j3:m){
k<-k+1
cd[k]<-zvalue*sqrt(m*(m+1)/(6*n))
if(abs(mr[j1]-mr[j2])>cd[k] ) sig005[k]<-1 else sig005[k]<-0
mrj[k]<-mr[j1]
mrk[k]<-mr[j2]
}
}
print(cbind(mrj,mrk,cd,sig005))
}
}
if(opt==4){
r<-t(apply(x,1,rank))
sumr<-apply(r,2,sum)
sumr2<-sort(sumr)
m2<-1:m
lvalue<-sum(sumr2*m2)
print("Page Test Result")
print(cbind(lvalue))
if(cvalue==0){
num<-lvalue-(n*m*(m+1)^2/4)
den<-sqrt(n*(m^3-m)^2/(144*(m-1)))
zvalue<-num/den
pvalue<-1-pnorm(zvalue,0,1)
print(cbind(zvalue,pvalue))
}
if(cvalue !=0 ) {
if(lvalue>=cvalue) print("significant at 5% level ") else print(" not significant at 5% level")
}
if(lvalue>=cvalue || pvalue < 0.025){
r<-t(apply(x,1,rank))
sumr<-apply(r,2,sum)
mr<-apply(r,2,mean)
####multiple comparison(Bonferroni-Dunn method)####
alpha<-0.05
alpha2<-alpha/(m*(m-1))*2
zvalue<--qnorm(alpha2,0,1)
m1<-m-1
print("Multiple Comparison(Bonferroni-Dunn method)")
k<-0
cd<-0
sig005<-0
mrj<-0
mrk<-0
for( j1 in 1:m1) {
j3<-j1+1
for( j2 in j3:m){
k<-k+1
cd[k]<-zvalue*sqrt(m*(m+1)/(6*n))
if(abs(mr[j1]-mr[j2])>cd[k] ) sig005[k]<-1 else sig005[k]<-0
mrj[k]<-mr[j1]
mrk[k]<-mr[j2]
}
}
cbind(mrj,mrk,cd,sig005)
}
}
}
##########pfa############################################################
pfa<-function(r){
n<-ncol(r)
r[row(r)==col(r)]<-0
commu1<-apply(abs(r),1,max)
r[row(r)==col(r)]<-commu1
r0<-r
dif<-1000
for ( i in 1:100){
r[row(r)==col(r)]<-commu1
evalue<-eigen(r)$values
evector<-eigen(r)$vectors
for ( j in 1:ncol(evector)) if(evector[1,j]<=0) evector[,j]<--evector[,j]
evalue1<-evalue[round(evalue,6)>0]
nc<-length(evalue1)
cont.evalue1<-evalue1/n
d1<-matrix(0,ncol=nc,nrow=nc)
d1[col(d1)>=row(d1)]<-1
cum.evalue<-cont.evalue1%*%d1
lamda<-diag(evalue1,nc)
evector<-evector[,1:nc]
b<-apply(evector^2,2,sum)
for ( j in 1:nc) evector[,j]<-evector[,j]/sqrt(b[j])
a<-evector%*%sqrt(lamda)
commu2<-apply(a^2,1,sum)
dif<-sqrt(mean((commu1-commu2)^2))
commu1<-commu2
if(dif < 0.01) {iter<-i;break}
}
par(mfrow=c(2,2),pty="s",las=1,lwd=2,tck=0.02)
number.axis<-c(1:nc)
plot(number.axis,evalue1,type="b")
plot(number.axis,cum.evalue,type="b")
axis1<-c(-1,1)
axis2<-c(-1,1)
plot(axis1,axis2,type="n")
text(a[,1],a[,2])
write.table(c("loading"),file="pfa-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(a,3),file="pfa-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(r.org=round(r0,3),r.final=round(r,3),iter=iter,
evalue=round(evalue1,3),cont=round(cont.evalue1,3),
cum=round(cum.evalue,3),loading=round(a,3),
communarity=round(commu2,3))
}
##########prcomp2#######################################################
prcomp2<-function(x){
scale2<-function(x)
{
if(is.matrix(x)) {
n<-nrow(x)
meanx<-apply(x,2,mean)
var2<-apply(x,2,var)*(n-1)/n
z<-matrix(0,ncol=ncol(x),nrow=nrow(x))
for( j in 1:ncol(x)) z[,j]<-(x[,j]-meanx[j])/sqrt(var2[j])
}
if(!is.matrix(x)) {
n<-length(x)
meanx<-mean(x)
var2<-var(x)*(n-1)/n
z<-(x-meanx)/sqrt(var2)
}
z
}
n1<-nrow(x)
n2<-ncol(x)
var1<-var(x)*(n1-1)/n1
cor1<-cor(x)
result<-eigen(var1)
evalue<-result$values
evalue3<-round(eigen(cor1)$values,6)
evector3<-eigen(cor1)$vectors
cum.evalue<-cumsum(evalue)
cum.cont<-cum.evalue/sum(evalue)
evector<-result$vectors
for ( j in 1:ncol(evector)) if(evector[1,j]<=0) evector[,j]<--evector[,j]
norm.evector<-sqrt(apply(evector^2,2,sum))
evector2<-evector/norm.evector
evector4<-evector3/sqrt(apply(evector3^2,2,sum))
loading<-evector4%*%diag(sqrt(evalue3))
z<-x%*%evector2
z.norm<-scale2(z)
par(mfrow=c(2,2),pty="s",las=1,lwd=2,tck=0.02)
plot(c(1:length(evalue)),evalue,type="b")
plot(c(1:length(evalue)),cum.cont,type="b")
plot(z.norm[,1],z.norm[,2],type="n")
text(z.norm[,1],z.norm[,2],font=1)
plot(loading[,1],loading[,2],type="n")
text(loading[,1],loading[,2],font=1)
write.table(c("coef"),file="prcomp2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(evector,3),file="prcomp2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
write.table(c("loading"),file="prcomp2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(loading,3),file="prcomp2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
write.table(c("prcompscore1"),file="prcomp2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(z,3),file="prcomp2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
write.table(c("prcompscore2"),file="prcomp2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(z.norm,3),file="prcomp2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(data=round(x,3),cov=round(var1,3),cor=round(cor1,3),coef=round(evector,3),
loading=round(loading,3),evalue=round(evalue,3),cum.cont=round(cum.cont,3),
prcompscore1=round(z,3),prcompscore2=round(z.norm,3))
}
######scale2#############################################################
scale2<-function(x){
if(is.matrix(x)) {
n<-nrow(x)
meanx<-apply(x,2,mean)
var2<-apply(x,2,var)*(n-1)/n
z<-matrix(0,ncol=ncol(x),nrow=nrow(x))
for( j in 1:ncol(x)) z[,j]<-(x[,j]-meanx[j])/sqrt(var2[j])
}
if(!is.matrix(x)) {
n<-length(x)
meanx<-mean(x)
var2<-var(x)*(n-1)/n
z<-(x-meanx)/sqrt(var2)
}
z
}
#####sort.list2.1##########################################################
sort.list2.1<-function(x,y){
n<-length(x)
outx<-sort(x)
outseq<-sort.list(x)
outy<-y
if(is.vector(y)) outy<-y[outseq]
if(is.matrix(y)) {
nc<-ncol(y)
for ( j in 1:nc) outy[,j]<-y[outseq,j]
}
outy
}
##########sreg###########################################################
sreg<-function(x,y){
coef<-lsfit(x,y)$coef
yh<-coef[1]+coef[2]*x
RR<-cor(y,yh)^2
plot(x,y)
abline(coef)
list(coef=coef,yh=yh,RR=RR)
}
##########Tukey.test###################################################
Tukey.test<-function(data,ma,mb,mc,n,rep,opt)
{
na<-ma;nb<-mb;nc<-mc
value<--1
x0<-array(data,c(n,nc,nb,na));x<-aperm(x0,c(4,3,2,1));
x01<-array(0,c(na,nb,nc,n));x01[x!=value]<-1
x[x==value]<-0;n2<-ncol(data);n3<-nc*nb*na;n4<-nc*n;n5<-n3*n;select<-1
if(n2==n3){
nijk<-apply(x01,c(1,2,3),sum);nijk2<-apply(x01,c(1,2,4),sum)
nikk2<-apply(x01,c(1,3,4),sum);njkk2<-apply(x01,c(2,3,4),sum)
ni<-apply(nijk,1,sum);nj<-apply(nijk,2,sum)
nk<-apply(nijk,3,sum);nij<-apply(nijk,c(1,2),sum)
nik<-apply(nijk,c(1,3),sum);njk<-apply(nijk,c(2,3),sum)
na2<-nijk[,1,1];nab<-nijk[,,1];nbc<-nijk[1,,];nac<-nijk[,1,]
hmean<-na*nb*nc/(sum(1/nijk));gmean<-sum(x)/sum(x01)
meanxijk<-apply(x,c(1,2,3),sum)/nijk
meanxijk2<-apply(x,c(1,2,4),sum)/nijk2
meanxikk2<-apply(x,c(1,3,4),sum)/nikk2
meanxjkk2<-apply(x,c(2,3,4),sum)/njkk2
meanxij<-apply(meanxijk,c(1,2),mean)
meanxik<-apply(meanxijk,c(1,3),mean)
meanxik2<-apply(meanxijk2,c(1,3),mean)
meanxjk<-apply(meanxijk,c(2,3),mean)
meanxjk2<-apply(meanxijk2,c(2,3),mean)
meanxkk2<-apply(meanxikk2,c(2,3),mean)
meanxi<-apply(meanxij,1,mean)
meanxj<-apply(meanxij,2,mean)
meanxk<-apply(meanxjk,2,mean)
meanxk2<-apply(meanxik2,2,mean)
SSbetween<-0;SSwithin<-0;SSa<-0
SSb<-0;SSc<-0;SSab<-0;SSac<-0;SSbc<-0
SSabc<-0;SSe<-0;SSe.a<-0
SSe.b<-0;SSe.c<-0
SSe.ab<-0;SSe.ac<-0
SSe.bc<-0;SSe.abc<-0;SSe.w<-0;SSt<-0
abcs<-sum(x^2)
abs1<-sum(apply(x,c(1,2,4),sum)^2/(nc))
a<-sum(apply(x,1,sum)^2/(n*nb*nc))
b<-sum(apply(x,2,sum)^2/(n*na*nc))
c1<-sum(apply(x,3,sum)^2/(n*na*nb))
ab<-sum(apply(x,c(1,2),sum)^2/(nab*nc))
bc<-sum(apply(x,c(2,3),sum)^2/(nbc*na))
ac<-sum(apply(x,c(1,3),sum)^2/(nac*nb) )
abc<-sum(apply(x,c(1,2,3),sum)^2/nijk)
xpp<-sum(x)^2/(n*na*nb*nc)
abcp<-sum(meanxijk^2)
g<-sum(x)
gp<-sum(meanxijk)
xp<-gp^2/(na*nb*nc)
ap<-nb*nc*sum(meanxi^2)
bp<-na*nc*sum(meanxj^2)
abp<-nc*sum(meanxij^2)
cp<-na*nb*sum(meanxk^2)
acp<-nb*sum(meanxik^2)
bcp<-na*sum(meanxjk^2)
if(rep>=2)
{
a<-sum(apply(x,1,sum)^2/(na2*nb*nc))
as<-sum(apply(x,c(1,4),sum)^2/(nb*nc))
ab<-sum(apply(x,c(1,2),sum)^2/(na2*nc))
abs1<-sum(apply(x,c(1,2,4),sum)^2/(nc))
ac<-sum(apply(x,c(1,3),sum)^2/(na2*nb))
acs<-sum(apply(x,c(1,3,4),sum)^2/(nb))
}
if(rep==3)
{
s<-sum(apply(x,4,sum)^2/(na*nb*nc))
bs<-sum(apply(x,c(2,4),sum)^2/(na*nc))
cs<-sum(apply(x,c(3,4),sum)^2/(na*nb))
bcs<-sum(apply(x,c(2,3,4),sum)^2/na)
}
SSa<-hmean*(ap-xp)
SSb<-hmean*(bp-xp)
SSc<-hmean*(cp-xp)
SSab<-hmean*(abp-ap-bp+xp)
SSac<-hmean*(acp-ap-cp+xp)
SSbc<-hmean*(bcp-bp-cp+xp)
SSabc<-hmean*(abcp-abp-acp-bcp+ap+bp+cp-xp)
SSe<-abcs-abc
SSa.bj<-hmean*nc*(apply(meanxij^2,2,sum)-(apply(meanxij,2,sum))^2/na)
SSa.ck<-hmean*nb*(apply(meanxik^2,2,sum)-(apply(meanxik,2,sum))^2/na)
SSb.ai<-hmean*nc*(apply(meanxij^2,1,sum)-(apply(meanxij,1,sum))^2/nb)
SSb.ck<-hmean*na*(apply(meanxjk^2,2,sum)-(apply(meanxjk,2,sum))^2/nb)
SSc.ai<-hmean*nb*(apply(meanxik^2,1,sum)-(apply(meanxik,1,sum))^2/nc)
SSc.bj<-hmean*na*(apply(meanxjk^2,1,sum)-(apply(meanxjk,1,sum))^2/nc)
SSab.ck<-hmean*(apply(meanxijk^2,3,sum)-(apply(meanxijk,3,sum))^2/(na*nb))-SSa.ck-SSb.ck
SSac.bj<-hmean*(apply(meanxijk^2,2,sum)-(apply(meanxijk,2,sum))^2/(na*nc))-SSa.bj-SSc.bj
SSbc.ai<-hmean*(apply(meanxijk^2,1,sum)-(apply(meanxijk,1,sum))^2/(nb*nc))-SSb.ai-SSc.ai
SSa.bjck<-hmean*(apply(meanxijk^2,c(2,3),sum)-(apply(meanxijk,c(2,3),sum))^2/(na))
SSb.aick<-hmean*(apply(meanxijk^2,c(1,3),sum)-(apply(meanxijk,c(1,3),sum))^2/(nb))
SSc.aibj<-hmean*(apply(meanxijk^2,c(1,2),sum)-(apply(meanxijk,c(1,2),sum))^2/(nc))
for(i in 1:na)for(j in 1:nb)for(k in 1:nc)for(k2 in 1:n){
if(x01[i,j,k,k2]!=0) SSt<-SSt+(x[i,j,k,k2]-gmean)^2
}
if(rep==1){
SSe.between<-abs1-ab
SSe.within<-abcs-abc-abs1+ab
}
if(rep==2){
SSe.a<-as-a
SSe.between<-SSe.a
SSe.b<-abs1-ab-as+a
SSe.c<-acs-ac-as+a
SSe.bc<-abcs-abc-abs1-acs+ab+ac+as-a
}
if(rep==3){
SSe.between<-s-xpp
SSe.a<-as-a-s+xpp
SSe.b<-bs-b-s+xpp
SSe.c<-cs-c1-s+xpp
SSe.ab<-abs1-ab-as-bs+a+b+s-xpp
SSe.ac<-acs-ac-as-cs+a+c1+s-xpp
SSe.bc<-bcs-bc-bs-cs+b+c1+s-xpp
SSe.abc<-abcs-abc-abs1-acs-bcs+ab+ac+bc+as+bs+cs-a-b-c1-s+xpp
}
if(na!=1) dfa<-na-1
if(rep==2) {dfe.between<-sum(na2)-na-1
dfe.within<-(sum(na2)-na)*(nb*nc-1)
}
if(rep==3) {
dfe.between<-n-1;dfe.within<-n*(na*nb*nc-1)
}
dfa<-na-1;dfb<-nb-1;dfc<-nc-1;dfab<-(na-1)*(nb-1)
dfac<-(na-1)*(nc-1);dfbc<-(nb-1)*(nc-1)
dfabc<-(na-1)*(nb-1)*(nc-1)
if(rep==0) dfe<-sum(nijk)-na*nb*nc
if(rep==1) {dfe.within<-(sum(nab)-na*nb)*(nc-1)
dfe.between<-sum(nab)-na*nb
}
if(rep==2) {dfe.a<-sum(na2)-na
dfe.b<-(sum(na2)-na)*(nb-1)
dfe.c<-(sum(na2)-na)*(nc-1)
dfe.bc<-(sum(na2)-na)*(nb-1)*(nc-1)
}
if(rep==3) {dfe.a<-(na-1)*(n-1)
dfe.b<-(nb-1)*(n-1)
dfe.c<-(nc-1)*(n-1)
dfe.ab<-(na-1)*(nb-1)*(n-1)
dfe.ac<-(na-1)*(nc-1)*(n-1)
dfe.bc<-(nb-1)*(nc-1)*(n-1)
dfe.abc<-(na-1)*(nb-1)*(nc-1)*(n-1)
}
dft<-sum(nijk)-1
if(rep==1) dft<-sum(nab)*nc-1
if(rep==2) dft<-sum(na2)*nb*nc-1
MSa<-SSa/dfa;MSb<-SSb/dfb;MSc<-SSc/dfc
MSab<-SSab/dfab;MSac<-SSac/dfac
MSbc<-SSbc/dfbc;MSabc<-SSabc/dfabc;
if(rep==0) MSe<-SSe/dfe
if(rep!=0) MSe.between<-SSe.between/dfe.between
if(rep==1) MSe.within<-SSe.within/dfe.within
if(rep>1) {MSe.a<-SSe.a/dfe.a;MSe.b<-SSe.b/dfe.b
MSe.c<-SSe.c/dfe.c;MSe.bc<-SSe.bc/dfe.bc
}
if(rep==3) {MSe.ab<-SSe.ab/dfe.ab;MSe.ac<-SSe.ac/dfe.ac;
MSe.abc<-SSe.abc/dfe.abc
}
if(rep==0){Fa<-MSa/MSe;Fab<-MSab/MSe;Fac<-MSac/MSe
Fabc<-MSabc/MSe;Fb<-MSb/MSe;Fc<-MSc/MSe;Fbc<-MSbc/MSe
}
if(rep>1) Fa<-MSa/MSe.a
if(rep==1) {Fa<-MSa/MSe.between;Fb<-MSb/MSe.between
Fc<-MSc/MSe.within;Fab<-MSab/MSe.between
Fac<-MSac/MSe.within;Fbc<-MSbc/MSe.within;Fabc<-MSabc/MSe.within
}
if(rep>1) {Fb<-MSb/MSe.b;Fc<-MSc/MSe.c;Fbc<-MSbc/MSe.bc
}
if(rep==2) {Fab<-MSab/MSe.b;Fac<-MSac/MSe.c;Fabc<-MSabc/MSe.bc
}
if(rep==3) {Fab<-MSab/MSe.ab;Fac<-MSac/MSe.ac;Fabc<-MSabc/MSe.abc
}
if(rep==0){
if(na !=1) {Pa<-1-pf(Fa,dfa,dfe);Pab<-1-pf(Fab,dfab,dfe)
Pac<-1-pf(Fac,dfac,dfe);Pabc<-1-pf(Fabc,dfabc,dfe)
}
if(nb !=1) {Pb<-1-pf(Fb,dfb,dfe);Pbc<-1-pf(Fbc,dfbc,dfe)
}
Pc<-1-pf(Fc,dfc,dfe)
}
if(rep==1 && na!=1) Pa<-1-pf(Fa,dfa,dfe.between)
if(rep>1 && na!=1) {Pa<-1-pf(Fa,dfa,dfe.a);Pb<-1-pf(Fb,dfb,dfe.b)
Pc<-1-pf(Fc,dfc,dfe.c);Pbc<-1-pf(Fbc,dfbc,dfe.bc)
}
if(rep>1 && na==1 && nb!=1) {Pb<-1-pf(Fb,dfb,dfe.b)
Pc<-1-pf(Fc,dfc,dfe.c);Pbc<-1-pf(Fbc,dfbc,dfe.bc)
}
# if(rep==1 && nb!=1) {Pb<-1-pf(Fb,dfe.between,dfe.within)
if(rep==1 && nb!=1) {Pb<-1-pf(Fb,dfb,dfe.between) #20070130
Pc<-1-pf(Fc,dfc,dfe.within);Pbc<-1-pf(Fbc,dfbc,dfe.within)
print(cbind(Fc,dfc,dfe.within,Pc)) ##20070129
if(na!=1) {Pabc<-1-pf(Fabc,dfabc,dfe.within);
Pab<-1-pf(Fab,dfab,dfe.between);Pac<-1-pf(Fac,dfac,dfe.within)
}
}
if(rep==1 && nb==1) Pc<-1-pf(Fc,dfc,dfe.within)
if(rep==2 && na!=1) {Pab<-1-pf(Fab,dfab,dfe.c)
Pac<-1-pf(Fac,dfac,dfe.c)
Pabc<-1-pf(Fabc,dfabc,dfe.bc)
}
if(rep==3) {Pab<-1-pf(Fab,dfab,dfe.ab);Pac<-1-pf(Fac,dfac,dfe.ac)
Pabc<-1-pf(Fabc,dfabc,dfe.abc)
}
print("data");print(data);
if(nb!=1 && na!=1) {print("meanxijk");print(meanxijk)
print("meanxij");print(meanxij)
print("meanxik");print(meanxik);print("meanxjk")
print(meanxjk);print("meanxi");print(meanxi)
print("meanxj");print(meanxj)
}
if(nb!=1 && na==1) {print("meanxj");print(meanxij)
print("meanxk");print(meanxik);print("meanxjk")
print(meanxjk)
}
if(nb==1 && na==1) {print("meanxk");print(meanxik)
}
if(nb==1 && na==1 && n4!=sum(nijk)){
select<-2
y<-sum(x)^2/sum(x01)
as2<-sum(x^2)
a2<-sum(apply(x[1,,,],1,sum)^2/nijk)
SSc<-a2-y
SSe<-as2-a2
SSe.c<-SSe
dfc<-nc-1
dfe<-sum(x01)-nc
dfe.c<-dfe
MSc<-SSc/dfc
MSe<-SSe/dfe
MSe.c<-MSe
Fc<-MSc/MSe
Pc<-1-pf(Fc,dfc,dfe)
}
###########Tukey Test####
alpha<-0.05
if(na>=1){
if(na>2 && Pa >=alpha) print.noquote("Factor A is not significant at 5% level")
if(na >2 && Pa < alpha){
print.noquote("Factor A is significant at 5% level")
cl<-1-alpha
if(rep == 0) {
df<-dfe
MSw<-MSe
}
if(rep == 1) {
dfe.ab<-dfe.between
MSe.ab<-MSe.between
df<-dfe.ab
MSw<-MSe.ab
}
if(rep >= 2) {df<-dfe.a
MSw<-MSe.a
}
mx<-meanxi
k<-length(mx)
q2<-qtukey(cl,k,df);mx2<-sort(mx);ni2<-ni[sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
q3<-qtukey(cl,m,df);if(opt==2) q3<-q2;qvalue<-(q2+q3)/2
n12<-ni2[i];n22<-ni2[j]
nhm<-2/(nb*nc/n12+nb*nc/n22)
WSD<-qvalue*sqrt(MSw/(nhm*nb*nc))
nhm2<-nhm*nb*nc
if(d > WSD) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm2)
}
}
if(opt==1) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],WSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],dfw=result[,11],
sig005=result[,12])
if(opt==2) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],HSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],
dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison(WSD Test) of factor A");
# write.table(c("Multiple Comparison(WSD Test) of factor A"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison(HSD Test) of factor A")
# write.table(c("Multiple Comparison(HSD Test) of factor A"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(nb > 2 && Pb >=alpha) print.noquote("Factor B is not significant at 5% level")
if(nb > 2 && Pb < alpha){
print.noquote("Factor B is significant at 5% level")
cl<-1-alpha
if(rep == 0) {df<-dfe
MSw<-MSe
}
if(rep == 1) {
dfe.ab<-dfe.between
df<-dfe.ab
MSe.ab<-MSe.between
MSw<-MSe.ab
}
if(rep >= 2) {df<-dfe.b
MSw<-MSe.b
}
mx<-meanxj
k<-length(mx)
q2<-qtukey(cl,k,df);mx2<-sort(mx);nj2<-nj[sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
q3<-qtukey(cl,m,df);if(opt==2) q3<-q2;qvalue<-(q2+q3)/2
n12<-nj2[i];n22<-nj2[j]
nhm<-2/(na*nc/n12+na*nc/n22)
nhm2<-nhm*na*nc
WSD<-qvalue*sqrt(MSw/(nhm*na*nc))
if(d>WSD) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm2)
}
}
if(opt==1) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],WSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],dfw=result[,11],
sig005=result[,12])
if(opt==2) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],HSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],dfw=result[,11],
sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor B")
# write.table(c("Multiple Comparison(WSD Test) of factor B"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison(HSD Test) of factor B")
# write.table(c("Multiple Comparison(HSD Test) of factor B"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(nc > 2 && Pc >=alpha) print.noquote("Factor C is not significant at 5% level")
if(nc > 2 && Pc < alpha){
print.noquote("Factor C is significant at 5% level")
cl<-1-alpha
if(rep == 0) {df<-dfe
MSw<-MSe
}
if(rep==1) {dfe.c<-dfe.within
MSe.c<-MSe.within
df<-dfe.c
MSw<-MSe.c
}
if(rep >= 2) {df<-dfe.c
MSw<-MSe.c
}
mx<-meanxk
k<-length(mx)
q2<-qtukey(cl,k,df);mx2<-sort(mx);nk2<-nk[sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
q3<-qtukey(cl,m,df);if(opt==2) q3<-q2;qvalue<-(q2+q3)/2
n12<-nk2[i];n22<-nk2[j]
nhm<-2/(na*nb/n12+na*nb/n22)
nhm2<-nhm*na*nb
WSD<-qvalue*sqrt(MSw/(nhm*na*nb))
if(d>WSD) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm2)
}
}
if(opt==1) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],WSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],dfw=result[,11],
sig005=result[,12])
if(opt==2) result<-data.frame(meanj=result[,1],meank=result[,2],d=result[,3],HSD=result[,4],m=result[,5],
h=result[,6],qvalue=result[,7],MSw=result[,8],nj=result[,9],nk=result[,10],nhm2=result[,13],dfw=result[,11],
sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor C")
# write.table(c("Multiple Comparison(WSD Test) of factor C"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison(HSD Test) of factor C")
# write.table(c("Multiple Comparison(HSD Test) of factor C"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(na!=1 && Pab >=alpha) print.noquote("Interaction AB is not significant at 5% level")
if(na!=1 && Pab < alpha){
print.noquote("Interaction AB is significant at 5% level")
if(rep==2) {dfpool.a.b<-dfe.a+dfe.b
MSpool.a.b<-(SSe.a+SSe.b)/dfpool.a.b
}
if(rep==3) {dfpool.a.ab<-dfe.a+dfe.ab
dfpool.b.ab<-dfe.b+dfe.ab
MSpool.a.ab<-(SSe.a+SSe.ab)/dfpool.a.ab
MSpool.b.ab<-(SSe.b+SSe.ab)/dfpool.a.ab
}
dfa.bj<-na-1
MSa.bj<-SSa.bj/dfa.bj
if(rep==0) {Fa.bj<-MSa.bj/MSe
Pa.bj<-1-pf(Fa.bj,dfa.bj,dfe)
}
if(rep==1) {
dfe.ab<-dfe.between
MSe.ab<-MSe.between
Fa.bj<-MSa.bj/MSe.ab
Pa.bj<-1-pf(Fa.bj,dfa.bj,dfe.ab)
}
if(rep==2) {Fa.bj<-MSa.bj/MSpool.a.b
Pa.bj<-1-pf(Fa.bj,dfa.bj,dfpool.a.b)
}
if(rep==3) {Fa.bj<-MSa.bj/MSpool.a.ab
Pa.bj<-1-pf(Fa.bj,dfa.bj,dfpool.a.ab)
}
print.noquote("Simple Main Effect of AB Interaction")
# write.table(c("Simple Main Effect of AB Interaction"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
print(data.frame(SSa.bj=SSa.bj,dfa.bj=dfa.bj,MSa.bj=MSa.bj,Fa.bj=Fa.bj,Pa.bj=round(Pa.bj,8)))
# write.table(data.frame(SSa.bj=SSa.bj,dfa.bj=dfa.bj,MSa.bj=MSa.bj,Fa.bj=Fa.bj,
# Pa.bj=round(Pa.bj,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfe.ab=dfe.ab,MSe.ab=MSe.ab))
if(rep==2) print(data.frame(dfpool.a.b=dfpool.a.b,MSpool.a.b=MSpool.a.b))
if(rep==3) print(data.frame(dfpool.a.ab=dfpool.a.ab,MSpool.a.ab=MSpool.a.ab))
dfb.ai<-nb-1
MSb.ai<-SSb.ai/dfb.ai
dfe.b<-(nb-1)*(n-1)
MSe.b<-SSe.b/dfe.b
if(rep==0) {Fb.ai<-MSb.ai/MSe
Pb.ai<-1-pf(Fb.ai,dfb.ai,dfe)
}
if(rep==1) {Fb.ai<-MSb.ai/MSe.ab
Pb.ai<-1-pf(Fb.ai,dfb.ai,dfe.ab)
}
if(rep==2) {Fb.ai<-MSb.ai/MSe.b
Pb.ai<-1-pf(Fb.ai,dfb.ai,dfe.b)
}
if(rep==3) {Fb.ai<-MSb.ai/MSpool.b.ab
Pb.ai<-1-pf(Fb.ai,dfb.ai,dfpool.b.ab)
}
print(data.frame(SSb.ai=SSb.ai,dfb.ai=dfb.ai,MSb.ai=MSb.ai,Fb.ai=Fb.ai,Pb.ai=round(Pb.ai,8)))
# write.table(data.frame(SSb.ai=SSb.ai,dfb.ai=dfb.ai,MSb.ai=MSb.ai,Fb.ai=Fb.ai,
# Pb.ai=round(Pb.ai,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfe.ab=dfe.ab,MSe.ab=MSe.ab))
if(rep==2) print(data.frame(dfe.b=dfe.b,MSe.b=MSe.b))
if(rep==3) print(data.frame(dfpool.b.ab=dfpool.b.ab,MSpool.b.ab=MSpool.b.ab))
if(na > 2) for( jj in 1:nb){
if(Pa.bj[jj] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1) {
df<-dfe.ab
MSw<-MSe.ab
}
if(rep==2) {df<-dfe.a+dfe.b
MSw<-(SSe.a+SSe.b)/df
}
if(rep==3) {
dfe.a.ab<-dfe.a+dfe.ab
df<-dfe.a.ab
MSw<-MSpool.a.ab
}
mx<-meanxij[,jj]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-qtukey(cl,k,dfe.ab)
if(rep==2) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.b)*MSe.b*(nb-1))/(MSe.a+MSe.b*(nb-1))
if(rep==3) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.ab)*MSe.ab*(nb-1))/(MSe.a+MSe.ab*(nb-1))
mx2<-sort(mx);nij2<-nij;nij2[,jj]<-nij[sort.list(mx),jj]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-qtukey(cl,m,dfe.ab)
if(rep==2) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.b)*MSe.b*(nb-1))/(MSe.a+MSe.b*(nb-1))
if(rep==3) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.ab)*MSe.ab*(nb-1))/(MSe.a+MSe.ab*(nb-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nij2[i,jj];n22<-nij2[j,jj]
nhm<-2/(nc/n12+nc/n22)
nhm2<-nhm*nc
WSD<-qvalue*sqrt(MSw/(nhm*nc))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,jj,nhm2)
}
}
if(opt==1) result<-data.frame(Bj=result[,13],meanBu=result[,1],meanBv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nBu=result[,9],nBv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(Bj=result[,13],meanBu=result[,1],meanBv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nBu=result[,9],nBv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor A[Bj]")
# write.table(c("Multiple Comparison (WSD Test) of factor A[Bj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison (HSD Test) of factor A[Bj]")
# write.table(c("Multiple Comparison (HSD Test) of factor A[Bj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
if( nb > 2) for( ii in 1:na){
if(Pb.ai[ii] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1) {
dfe.ab<-dfe.between
MSe.ab<-MSe.between
df<-dfe.ab
MSw<-MSe.ab
}
if(rep==2) {df<-dfe.b
MSw<-SSe.b/df
}
if(rep==3) {
dfe.b.ab<-dfe.b+dfe.ab
df<-dfe.b.ab
MSpool.b.ab<-(SSe.b+SSe.ab)/dfe.b.ab
MSw<-MSpool.b.ab
}
mx<-meanxij[ii,]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-qtukey(cl,k,dfe.ab)
if(rep==2) q2<-qtukey(cl,k,dfe.a)
if(rep==3) q2<-(qtukey(cl,k,dfe.b)*MSe.b+qtukey(cl,k,dfe.ab)*MSe.ab*(na-1))/(MSe.b+MSe.ab*(na-1))
mx2<-sort(mx);nij2<-nij;nij2[ii,]<-nij[ii,sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-qtukey(cl,m,dfe.ab)
if(rep==2) q3<-qtukey(cl,m,dfe.a)
if(rep==3) q3<-(qtukey(cl,m,dfe.b)*MSe.b+qtukey(cl,m,dfe.ab)*MSe.ab*(na-1))/(MSe.b+MSe.ab*(na-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nij2[ii,i];n22<-nij2[ii,j]
nhm<-2/(nc/n12+nc/n22)
nhm2<-nhm*nc
WSD<-qvalue*sqrt(MSw/(nhm*nc))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,ii,nhm2)
}
}
if(opt==1) result<-data.frame(i=result[,13],meanAu=result[,1],meanAv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(i=result[,13],meanAu=result[,1],meanAv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor B[Ai]")
# write.table(c("Multiple Comparison (WSD Test) of factor B[Ai]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison(HSD Test) of factor B[Ai]")
# write.table(c("Multiple Comparison (HSD Test) of factor B[Ai]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
}
if(na!=1 && Pac >=alpha) print.noquote("Interaction AC is not significant at 5% level")
if(na!=1 && Pac < alpha){
print.noquote("Interaction AC is significant at 5% level")
if(rep==1) {
dfe.c<-dfe.within
MSe.c<-MSe.within
dfe.ab<-dfe.between
dfpool.c<-dfe.ab+dfe.c
SSe.between<-abs1-ab
SSe.within<-abcs-abc-abs1+ab
SSe.ab<-SSe.between
SSe.c<-SSe.within
dfpool.c<-dfe.ab+dfe.c
MSpool.c<-(SSe.ab+SSe.c)/dfpool.c
}
if(rep==2) {dfpool.a.c<-dfe.a+dfe.c
MSpool.a.c<-(SSe.a+SSe.c)/dfpool.a.c
}
if(rep==3) {dfpool.a.ac<-dfe.a+dfe.ac
dfpool.c.ac<-dfe.c+dfe.ac
MSpool.a.ac<-(SSe.a+SSe.ac)/dfpool.a.ac
MSpool.c.ac<-(SSe.c+SSe.ac)/dfpool.c.ac
}
dfa.ck<-na-1
MSa.ck<-SSa.ck/dfa.ck
if(rep==0) {Fa.ck<-MSa.ck/MSe
Pa.ck<-1-pf(Fa.ck,dfa.ck,dfe)
}
if(rep==1) {Fa.ck<-MSa.ck/MSpool.c
Pa.ck<-1-pf(Fa.ck,dfa.ck,dfpool.c)
}
if(rep==2) {Fa.ck<-MSa.ck/MSpool.a.c
Pa.ck<-1-pf(Fa.ck,dfa.ck,dfpool.a.c)
}
if(rep==3) {Fa.ck<-MSa.ck/MSpool.a.ac
Pa.ck<-1-pf(Fa.ck,dfa.ck,dfpool.a.ac)
}
print.noquote("Simple Main Effect of AC Interaction")
# write.table(c("Simple Main Effect of AC Interaction"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
print(data.frame(SSa.ck=SSa.ck,dfa.ck=dfa.ck,MSa.ck=MSa.ck,Fa.ck=Fa.ck,Pa.ck=round(Pa.ck,8)))
# write.table(data.frame(SSa.ck=SSa.ck,dfa.ck=dfa.ck,MSa.ck=MSa.ck,Fa.ck=Fa.ck,
# Pa.ck=round(Pa.ck,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfpool.c=dfpool.c,MSpool.c=MSpool.c))
if(rep==2) print(data.frame(dfpool.a.c=dfpool.a.c,MSpool.a.c=MSpool.a.c))
if(rep==3) print(data.frame(dfpool.a.ac=dfpool.a.ac,MSpool.a.ac=MSpool.a.ac))
dfc.ai<-nc-1
MSc.ai<-SSc.ai/dfc.ai
if(rep==0) {Fc.ai<-MSc.ai/MSe
Pc.ai<-1-pf(Fc.ai,dfc.ai,dfe)
}
if(rep==1 || rep==2) {Fc.ai<-MSc.ai/MSe.c
Pc.ai<-1-pf(Fc.ai,dfc.ai,dfe.c)
}
if(rep==3) {Fc.ai<-MSc.ai/MSpool.c.ac
Pc.ai<-1-pf(Fc.ai,dfc.ai,dfpool.c.ac)
}
print(data.frame(SSc.ai=SSc.ai,dfc.ai=dfc.ai,MSc.ai=MSc.ai,Fc.ai=Fc.ai,Pc.ai=round(Pc.ai,8)))
# write.table(data.frame(SSc.ai=SSc.ai,dfc.ai=dfc.ai,MSc.ai=MSc.ai,Fc.ai=Fc.ai,
# Pc.ai=round(Pc.ai,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1 || rep==2) print(data.frame(dfe.c=dfe.c,MSe.c=MSe.c))
if(rep==3) print(data.frame(dfpool.c.ac=dfpool.c.ac,MSpool.c.ac=MSpool.c.ac))
if(na > 2) for( kk in 1:nc){
if(Pa.ck[kk] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1) {df<-dfpool.c
MSw<-MSpool.c
dfe.ab<-dfe.between
MSe.ab<-MSe.between
}
if(rep==2) {df<-dfe.a+dfe.c
MSw<-(SSe.a+SSe.c)/df
}
if(rep==3) {df<-dfpool.a.ac
MSw<-MSpool.a.ac
}
mx<-meanxik[,kk]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-(qtukey(cl,k,dfe.ab)*MSe.ab+qtukey(cl,k,dfe.c)*MSe.c*(nc-1))/(MSe.ab+MSe.c*(nc-1))
if(rep==2) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.c)*MSe.c*(nc-1))/(MSe.a+MSe.c*(nc-1))
if(rep==3) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.ac)*MSe.ac*(nc-1))/(MSe.a+MSe.ac*(nc-1))
mx2<-sort(mx);nik2<-nik;nik2[,kk]<-nik[sort.list(mx),kk]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-(qtukey(cl,m,dfe.ab)*MSe.ab+qtukey(cl,m,dfe.c)*MSe.c*(nc-1))/(MSe.ab+MSe.c*(nc-1))
if(rep==2) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.c)*MSe.c*(nc-1))/(MSe.a+MSe.c*(nc-1))
if(rep==3) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.ac)*MSe.ac*(nc-1))/(MSe.a+MSe.ac*(nc-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nik2[i,kk];n22<-nik2[j,kk]
nhm<-2/(nb/n12+nb/n22)
nhm2<-nhm*nb
WSD<-qvalue*sqrt(MSw/(nhm*nb))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,kk,nhm2)
}
}
if(opt==1) result<-data.frame(Ck=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(Ck=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor A[Ck]")
# write.table(c("Multiple Comparison (WSD Test) of factor A[Ck]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison (HSD Test) of factor A[Ck]")
# write.table(c("Multiple Comparison (HSD Test) of factor A[Ck]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
if( nc > 2) for( ii in 1:na){
if(Pc.ai[ii] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1 || rep==2) {df<-dfe.c
MSw<-MSe.c
}
if(rep==3) {df<-dfe.c+dfe.ac
MSw<-(SSe.c+SSe.ac)/df
}
mx<-meanxik[ii,]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1 || rep==2) q2<-qtukey(cl,k,dfe.c)
if(rep==3) q2<-(qtukey(cl,k,dfe.c)*MSe.c+qtukey(cl,k,dfe.ac)*MSe.ac*(na-1))/(MSe.c+MSe.ac*(na-1))
mx2<-sort(mx);nik2<-nik;nik2[ii,]<-nik[ii,sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1 || rep==2) q3<-qtukey(cl,m,dfe.c)
if(rep==3) q3<-(qtukey(cl,m,dfe.c)*MSe.c+qtukey(cl,m,dfe.ac)*MSe.ac*(na-1))/(MSe.c+MSe.ac*(na-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nik2[ii,i];n22<-nik2[ii,j]
nhm<-2/(nc/n12+nc/n22)
nhm2<-nhm*nc
WSD<-qvalue*sqrt(MSw/(nhm*nc))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,ii,nhm2)
}
}
if(opt==1) result<-data.frame(Ai=result[,13],meanAu=result[,1],meanAv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(Ai=result[,13],meanAu=result[,1],meanAv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor C[Ai]")
# write.table(c("Multiple Comparison (WSD Test) of factor C[Ai]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison (HSD Test) of factor C[Ai]")
# write.table(c("Multiple Comparison (HSD Test) of factor C[Ai]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
}
if(nb!=1 && Pbc >=alpha) print.noquote("Interaction BC is not significant at 5% level")
if(nb!= 1 && Pbc < alpha){
print.noquote("Interaction BC is significant at 5% level")
if(rep==2 || rep==3) {dfpool.b.bc<-dfe.b+dfe.bc
dfpool.c.bc<-dfe.c+dfe.bc
MSpool.b.bc<-(SSe.b+SSe.bc)/dfpool.b.bc
MSpool.c.bc<-(SSe.c+SSe.bc)/dfpool.c.bc
}
dfb.ck<-nb-1
MSb.ck<-SSb.ck/dfb.ck
if(rep==0) {
Fb.ck<-MSb.ck/MSe
Pb.ck<-1-pf(Fb.ck,dfb.ck,dfe)
}
if(rep==1) {
dfe.ab<-dfe.between
dfe.c<-dfe.within
MSe.c<-MSe.within
dfpool.c<-dfe.ab+dfe.c
SSe.between<-abs1-ab
SSe.within<-abcs-abc-abs1+ab
SSe.ab<-SSe.between
SSe.c<-SSe.within
MSpool.c<-(SSe.ab+SSe.c)/dfpool.c
Fb.ck<-MSb.ck/MSpool.c
Pb.ck<-1-pf(Fb.ck,dfb.ck,dfpool.c)
}
if(rep==2 || rep==3) {Fb.ck<-MSb.ck/MSpool.b.bc
Pb.ck<-1-pf(Fb.ck,dfb.ck,dfpool.b.bc)
}
print.noquote("Simple Main Effect of BC Interaction")
# write.table(c("Simple Main Effect of BC Interaction"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
print(data.frame(SSb.ck=SSb.ck,dfb.ck=dfb.ck,MSb.ck=MSb.ck,Fb.ck=Fb.ck,Pb.ck=round(Pb.ck,8)))
# write.table(data.frame(SSb.ck=SSb.ck,dfb.ck=dfb.ck,MSb.ck=MSb.ck,Fb.ck=Fb.ck,
# Pb.ck=round(Pb.ck,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfpool.c=dfpool.c,MSpool.c=MSpool.c))
if(rep==2 || rep==3) print(data.frame(dfpool.b.bc=dfpool.b.bc,MSpool.b.bc=MSpool.b.bc))
dfc.bj<-nc-1
MSc.bj<-SSc.bj/dfc.bj
if(rep==0) {Fc.bj<-MSc.bj/MSe
Pc.bj<-1-pf(Fc.bj,dfc.bj,dfe)
}
if(rep==1) {Fc.bj<-MSc.bj/MSe.c
Pc.bj<-1-pf(Fc.bj,dfc.bj,dfe.c)
}
if(rep==2 || rep==3) {Fc.bj<-MSc.bj/MSpool.c.bc
Pc.bj<-1-pf(Fc.bj,dfc.bj,dfpool.c.bc)
}
print(data.frame(SSc.bj=SSc.bj,dfc.bj=dfc.bj,MSc.bj=MSc.bj,Fc.bj=Fc.bj,Pc.bj=round(Pc.bj,8)))
# write.table(data.frame(SSc.bj=SSc.bj,dfc.bj=dfc.bj,MSc.bj=MSc.bj,Fc.bj=Fc.bj,
# Pc.bj=round(Pc.bj,8)),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfe.c=dfe.c,MSe.c=MSe.c))
if(rep==2 || rep==3) print(data.frame(dfpool.c.bc=dfpool.c.bc,MSpool.c.bc=MSpool.c.bc))
if(nb > 2) for( kk in 1:nc){
if(Pb.ck[kk] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1) {df<-dfpool.c
SSe.ab<-abs1-ab
SSe.c<-abcs-abc-abs1+ab
MSpool.c<-(SSe.ab+SSe.c)/dfpool.c
MSw<-MSpool.c
dfe.ab<-dfe.between
MSe.ab<-MSe.between
}
if(rep>=2) {df<-dfe.b+dfe.c
MSw<-(SSe.b+SSe.c)/df
}
mx<-meanxjk[,kk]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-(qtukey(cl,k,dfe.ab)*MSe.ab+qtukey(cl,k,dfe.c)*MSe.c)/(MSe.ab+MSe.c)
if(rep>=2) q2<-(qtukey(cl,k,dfe.b)*MSe.b+qtukey(cl,k,dfe.bc)*MSe.bc*(nc-1))/(MSe.b+MSe.bc*(nc-1))
mx2<-sort(mx);njk2<-njk;njk2[,kk]<-njk[sort.list(mx),kk]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-(qtukey(cl,m,dfe.ab)*MSe.ab+qtukey(cl,m,dfe.c)*MSe.c)/(MSe.ab+MSe.c)
if(rep>=2) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.b)*MSe.b*(nb-1))/(MSe.a+MSe.b*(nb-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-njk2[i,kk];n22<-njk2[j,kk]
nhm<-2/(na/n12+na/n22)
nhm2<-nhm*na
WSD<-qvalue*sqrt(MSw/(nhm*na))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,kk,nhm2)
}
}
if(opt==1) result<-data.frame(Ck=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(Ck=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor B[Ck]")
# write.table(c("Multiple Comparison (WSD Test) of factor B[Ck]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison (HSD Test) of factor B[Ck]")
# write.table(c("Multiple Comparison (HSD Test) of factor B[Ck]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
if( nc > 2) for( jj in 1:nb){
if(Pc.bj[jj] < alpha) {
cl<-1-alpha
if(rep==0) {df<-dfe
MSw<-MSe
}
if(rep==1) {df<-dfe.c
MSw<-MSe.c
}
if(rep>=2) {df<-dfe.c+dfe.b
MSw<-(SSe.c+SSe.b)/df
}
mx<-meanxjk[jj,]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-qtukey(cl,k,dfe.c)
if(rep>=2) q2<-(qtukey(cl,k,dfe.c)*MSe.c+qtukey(cl,k,dfe.bc)*MSe.bc*(nb-1))/(MSe.c+MSe.bc*(nb-1));
mx2<-sort(mx);njk2<-njk;njk2[jj,]<-njk[jj,sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=14)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-qtukey(cl,m,dfe.c)
if(rep>=2) q3<-(qtukey(cl,m,dfe.c)*MSe.c+qtukey(cl,m,dfe.bc)*MSe.bc*(nb-1))/(MSe.c+MSe.bc*(nb-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-njk2[jj,i];n22<-njk2[jj,j]
nhm<-2/(na/n12+na/n22)
nhm2<-nhm*na
WSD<-qvalue*sqrt(MSw/(nhm*na))
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,jj,nhm2)
}
}
if(opt==1) result<-data.frame(Bj=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==2) result<-data.frame(Bj=result[,13],meanCu=result[,1],meanCv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm2=result[,14],dfw=result[,11],sig005=result[,12])
if(opt==1) {print.noquote("Multiple Comparison (WSD Test) of factor C[Bj]")
# write.table(c("Multiple Comparison (WSD Test) of factor C[Bj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(opt==2) {print.noquote("Multiple Comparison (HSD Test) of factor C[Bj]")
# write.table(c("Multiple Comparison (HSD Test) of factor C[Bj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
print(result)
# write.table(result,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
}
if(na!=1 && Pabc >=alpha) print.noquote("Interaction ABC is not significant at 5% level")
if(na!=1 && Pabc < alpha){
print.noquote("Interaction ABC is significant at 5% level")
print.noquote("Test of Interaction AB[Ck]")
# write.table(c("Test of Interaction AB[Ck]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
dfab.ck<-(na-1)*(nb-1)
MSab.ck<-SSab.ck/dfab.ck
if(rep==0) {
Fab.ck<-MSab.ck/MSe
Pab.ck<-round(1-pf(Fab.ck,dfab.ck,dfe),8)
print(cbind(SSab.ck,dfab.ck,MSab.ck,dfe,MSe,Fab.ck,Pab.ck))
# write.table(cbind(SSab.ck,dfab.ck,MSab.ck,dfe,MSe,Fab.ck,Pab.ck),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==1) {
SSe.ab<-abs1-ab
SSe.c<-abcs-abc-abs1+ab
MSpool.c<-(SSe.ab+SSe.c)/dfpool.c
Fab.ck<-MSab.ck/MSpool.c
Pab.ck<-round(1-pf(Fab.ck,dfab.ck,dfpool.c),8)
print(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.c,MSpool.c,Fab.ck,Pab.ck))
# write.table(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.c,MSpool.c,Fab.ck,Pab.ck),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==2) {dfpool.b.bc<-dfe.b+dfe.bc
MSpool.b.bc<-(SSe.b+SSe.bc)/dfpool.b.bc
Fab.ck<-MSab.ck/MSpool.b.bc
Pab.ck<-round(1-pf(Fab.ck,dfab.ck,dfpool.b.bc),8)
print(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.b.bc,MSpool.b.bc,Fab.ck,Pab.ck))
# write.table(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.b.bc,MSpool.b.bc,Fab.ck,Pab.ck),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==3) {dfpool.ab.abc<-dfe.ab+dfe.abc
MSpool.ab.abc<-(SSe.ab+SSe.abc)/dfpool.ab.abc
Fab.ck<-MSab.ck/MSpool.ab.abc
Pab.ck<-round(1-pf(Fab.ck,dfab.ck,dfpool.ab.abc),8)
print(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.ab.abc,MSpool.ab.abc,Fab.ck,Pab.ck))
# write.table(cbind(SSab.ck,dfab.ck,MSab.ck,dfpool.ab.abc,MSpool.ab.abc,Fab.ck,Pab.ck),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
print.noquote("Test of Interaction AC[Bj]")
# write.table(c("Test of Interaction AC[Bj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
dfac.bj<-(na-1)*(nc-1)
MSac.bj<-SSac.bj/rep(dfac.bj,nb)
if(rep==0) {
Fac.bj<-MSac.bj/MSe
Pac.bj<-round(1-pf(Fac.bj,dfac.bj,dfe),8)
print(cbind(SSac.bj,dfac.bj,MSac.bj,dfe,MSe,Fac.bj,Pac.bj))
# write.table(cbind(SSac.bj,dfac.bj,MSac.bj,dfe,MSe,Fac.bj,Pac.bj),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==1) {
Fac.bj<-MSac.bj/MSe.c
Pac.bj<-round(1-pf(Fac.bj,dfac.bj,dfe.c),8)
print(cbind(SSac.bj,dfac.bj,MSac.bj,dfe.c,MSe.c,Fac.bj,Pac.bj))
# write.table(cbind(SSac.bj,dfac.bj,MSac.bj,dfe.c,MSe.c,Fac.bj,Pac.bj),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==2) {dfpool.c.bc<-dfe.c+dfe.bc
MSpool.c.bc<-(SSe.c+SSe.bc)/dfpool.c.bc
Fac.bj<-MSac.bj/MSpool.c.bc
Pac.bj<-round(1-pf(Fac.bj,dfac.bj,dfpool.c.bc),8)
print(cbind(SSac.bj,dfac.bj,MSac.bj,dfpool.c.bc,MSpool.c.bc,Fac.bj,Pac.bj))
# write.table(cbind(SSac.bj,dfac.bj,MSac.bj,dfpool.c.bc,MSpool.c.bc,Fac.bj,Pac.bj),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==3) {dfpool.ac.abc<-dfe.ac+dfe.abc
MSpool.ac.abc<-(SSe.ac+SSe.abc)/dfpool.ac.abc
Fac.bj<-MSac.bj/MSpool.ac.abc
Pac.bj<-round(1-pf(Fac.bj,dfac.bj,dfpool.ac.abc),8)
print(cbind(SSac.bj,dfac.bj,MSac.bj,dfpool.ac.abc,MSpool.ac.abc,Fac.bj,Pac.bj))
# write.table(cbind(SSac.bj,dfac.bj,MSac.bj,dfpool.ac.abc,MSpool.ac.abc,Fac.bj,Pac.bj),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
print.noquote("Test of Interaction BC[Ai]")
# write.table(c("Test of Interaction BC[Ai]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
dfbc.ai<-(nb-1)*(nc-1)
MSbc.ai<-SSbc.ai/rep(dfbc.ai,na)
if(rep==0) {Fbc.ai<-MSbc.ai/MSe
Pbc.ai<-round(1-pf(Fbc.ai,dfbc.ai,dfe),8)
print(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfe,MSe,Fbc.ai,Pbc.ai))
}
if(rep==1) {Fbc.ai<-MSbc.ai/MSe.c
Pbc.ai<-round(1-pf(Fbc.ai,dfbc.ai,dfe.c),8)
print(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfe.c,MSe.c,Fbc.ai,Pbc.ai))
# write.table(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfe.c,MSe.c,Fbc.ai,Pbc.ai),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==2) {Fbc.ai<-MSbc.ai/MSe.bc
Pbc.ai<-round(1-pf(Fbc.ai,dfbc.ai,dfe.bc),8)
print(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfe.bc,MSe.bc,Fbc.ai,Pbc.ai))
# write.table(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfe.bc,MSe.bc,Fbc.ai,Pbc.ai),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==3) {dfpool.bc.abc<-dfe.bc+dfe.abc
MSpool.bc.abc<-(SSe.bc+SSe.abc)/dfpool.bc.abc
Fbc.ai<-MSbc.ai/MSpool.bc.abc
Pbc.ai<-round(1-pf(Fbc.ai,dfbc.ai,dfpool.bc.abc),8)
print(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfpool.bc.abc,MSpool.bc.abc,Fbc.ai,Pbc.ai))
# write.table(cbind(SSbc.ai,dfbc.ai,MSbc.ai,dfpool.bc.abc,MSpool.bc.abc,Fbc.ai,Pbc.ai),file="mct.xls",
# row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
print.noquote("Simple Simple Main Effect of Factor A[BjCk]")
# write.table(c("Simple Simple Main Effect of Factor A[BjCk]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
for (jj in 1:nb) for(kk in 1:nc) {
if(rep==0){
dfa.bjck<-na-1
MSa.bjck<-SSa.bjck/dfa.bjck
Fa.bjck<-MSa.bjck/MSe
Pa.bjck<-round(1-pf(Fa.bjck,dfa.bjck,dfe),8)
}
if(rep==1){
dfa.bjck<-na-1
MSa.bjck<-SSa.bjck/dfa.bjck
SSe.ab<-abs1-ab
SSe.c<-abcs-abc-abs1+ab
MSpool.c<-(SSe.ab+SSe.c)/dfpool.c
Fa.bjck<-MSa.bjck/MSpool.c
Pa.bjck<-round(1-pf(Fa.bjck,dfa.bjck,dfpool.c),8)
}
if(rep==2){dfpool.a.b.c.bc<-dfe.a+dfe.b+dfe.c+dfe.bc
MSpool.a.b.c.bc<-(SSe.a+SSe.b+SSe.c+SSe.bc)/dfpool.a.b.c.bc
dfa.bjck<-na-1
MSa.bjck<-SSa.bjck/dfa.bjck
Fa.bjck<-MSa.bjck/MSpool.a.b.c.bc
Pa.bjck<-round(1-pf(Fa.bjck,dfa.bjck,dfpool.a.b.c.bc),8)
}
if(rep==3){dfpool.a.ab.ac.abc<-dfe.a+dfe.ab+dfe.ac+dfe.abc
MSpool.a.ab.ac.abc<-(SSe.a+SSe.ab+SSe.ac+SSe.abc)/dfpool.a.ab.ac.abc
dfa.bjck<-na-1
MSa.bjck<-SSa.bjck/dfa.bjck
Fa.bjck<-MSa.bjck/MSpool.a.ab.ac.abc
Pa.bjck<-round(1-pf(Fa.bjck,dfa.bjck,dfpool.a.ab.ac.abc),8)
}
j<-1:nb
print(data.frame(Bj=jj,Ck=kk,SSa.bjck=SSa.bjck[jj,kk],dfa.bjck=dfa.bjck,
MSa.bjck=MSa.bjck[jj,kk],Fa.bjck=Fa.bjck[jj,kk],Pa.bjck=Pa.bjck[jj,kk]))
# write.table(data.frame(Bj=jj,Ck=kk,SSa.bjck=SSa.bjck[jj,kk],dfa.bjck=dfa.bjck,
# MSa.bjck=MSa.bjck[jj,kk],Fa.bjck=Fa.bjck[jj,kk],Pa.bjck=Pa.bjck[jj,kk]),
# file="mct.xls",row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfpool.c=dfpool.c,MSpool.c=MSpool.c))
if(rep==2) print(data.frame(dfpool.a.b.c.bc=dfpool.a.b.c.bc,MSpool.a.b.c.bc=MSpool.a.b.c.bc))
if(rep==3) print(data.frame(dfpool.a.ab.ac.abc=dfpool.a.ab.ac.abc,
MSpool.a.ab.ac.abc=MSpool.a.ab.ac.abc))
print.noquote("Simple Simple Main Effect of Factor B[AiCk]")
# write.table(c("Simple Simple Main Effect of Factor B[AiCk]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
for (ii in 1:na) for(kk in 1:nc) {
if(rep==0) {
dfb.aick<-nb-1
MSb.aick<-SSb.aick/dfb.aick
Fb.aick<-MSb.aick/MSe
Pb.aick<-round(1-pf(Fb.aick,dfb.aick,dfe),8)
}
if(rep==1) {
dfb.aick<-nb-1
MSb.aick<-SSb.aick/dfb.aick
Fb.aick<-MSb.aick/MSpool.c
Pb.aick<-round(1-pf(Fb.aick,dfb.aick,dfpool.c),8)
}
if(rep==2) {dfpool.b.bc<-dfe.b+dfe.bc
MSpool.b.bc<-(SSe.b+SSe.bc)/dfpool.b.bc
dfb.aick<-nb-1
MSb.aick<-SSb.aick/dfb.aick
Fb.aick<-MSb.aick/MSpool.b.bc
Pb.aick<-round(1-pf(Fb.aick,dfb.aick,dfpool.b.bc),8)
}
if(rep==3) {dfpool.b.ab.bc.abc<-dfe.b+dfe.ab+dfe.bc+dfe.abc
MSpool.b.ab.bc.abc<-(SSe.b+SSe.ab+SSe.bc+SSe.abc)/dfpool.b.ab.bc.abc
dfb.aick<-nb-1
MSb.aick<-SSb.aick/dfb.aick
Fb.aick<-MSb.aick/MSpool.b.ab.bc.abc
Pb.aick<-round(1-pf(Fb.aick,dfb.aick,dfpool.b.ab.bc.abc),8)
}
print(data.frame(Ai=ii,Ck=kk,SSb.aick=SSb.aick[ii,kk],dfb.aick=dfb.aick,
MSb.aick=MSb.aick[ii,kk],Fb.aick=Fb.aick[ii,kk],Pb.aick=Pb.aick[ii,kk]))
# write.table(data.frame(Ai=ii,Ck=kk,SSb.aick=SSb.aick[ii,kk],dfb.aick=dfb.aick,
# MSb.aick=MSb.aick[ii,kk],Fb.aick=Fb.aick[ii,kk],Pb.aick=Pb.aick[ii,kk]),
# file="mct.xls",row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfpool.c=dfpool.c,MSpool.c=MSpool.c))
if(rep==2) print(data.frame(dfpool.b.bc=dfpool.b.bc,MSpool.b.bc=MSpool.b.bc))
if(rep==3) print(data.frame(dfpool.b.ab.bc.abc=dfpool.b.ab.bc.abc,
MSpool.b.ab.bc.abc=MSpool.b.ab.bc.abc))
print.noquote("Simple Simple Main Effect of Factor C[AiBj]")
# write.table(c("Simple Simple Main Effect of Factor C[AiBj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
for( ii in 1:na) for(jj in 1:nb) {
if(rep==0) {
dfc.aibj<-nc-1
MSc.aibj<-SSc.aibj/dfc.aibj
Fc.aibj<-MSc.aibj/MSe
Pc.aibj<-round(1-pf(Fc.aibj,dfc.aibj,dfe),8)
}
if(rep==1) {
dfc.aibj<-nc-1
MSc.aibj<-SSc.aibj/dfc.aibj
Fc.aibj<-MSc.aibj/MSe.c
Pc.aibj<-round(1-pf(Fc.aibj,dfc.aibj,dfe.c),8)
}
if(rep==2) {dfpool.c.bc<-dfe.c+dfe.bc
MSpool.c.bc<-(SSe.c+SSe.bc)/dfpool.c.bc
dfc.aibj<-nc-1
MSc.aibj<-SSc.aibj/dfc.aibj
Fc.aibj<-MSc.aibj/MSpool.c.bc
Pc.aibj<-round(1-pf(Fc.aibj,dfc.aibj,dfpool.c.bc),8)
}
if(rep==3) {dfpool.c.ac.bc.abc<-dfe.c+dfe.ac+dfe.bc+dfe.abc
MSpool.c.ac.bc.abc<-(SSe.c+SSe.ac+SSe.bc+SSe.abc)/dfpool.c.ac.bc.abc
dfc.aibj<-nc-1
MSc.aibj<-SSc.aibj/dfc.aibj
Fc.aibj<-MSc.aibj/MSpool.c.ac.bc.abc
Pc.aibj<-round(1-pf(Fc.aibj,dfc.aibj,dfpool.c.ac.bc.abc),8)
}
j<-1:nb
print(data.frame(Ai=ii,Bj=jj,SSc.aibj=SSc.aibj[ii,jj],dfc.aibj=dfc.aibj,
MSc.aibj=MSc.aibj[ii,jj],Fc.aibj=Fc.aibj[ii,jj],Pc.aibj=Pc.aibj[ii,jj]))
# write.table(data.frame(Ai=ii,Bj=jj,SSc.aibj=SSc.aibj[ii,jj],dfc.aibj=dfc.aibj,
# MSc.aibj=MSc.aibj[ii,jj],Fc.aibj=Fc.aibj[ii,jj],Pc.aibj=Pc.aibj[ii,jj]),
# file="mct.xls",row.names=F,col.names=F,append=T,quote=F,sep="\t")
}
if(rep==0) print(data.frame(dfe=dfe,MSe=MSe))
if(rep==1) print(data.frame(dfe.c=dfe.c,MSe.c=MSe.c))
if(rep==2) print(data.frame(dfpool.c.bc=dfpool.c.bc,MSpool.c.bc=MSpool.c.bc))
if(rep==3) print(data.frame(dfpool.c.ac.bc.abc=dfpool.c.ac.bc.abc,
MSpool.c.ac.bc.abc=MSpool.c.ac.bc.abc))
}
if(na > 2) {print.noquote("Test of Multiple Comparison of Factor A[BjCk]")
# write.table(c("Test of Multiple Comparison of Factor A[BjCk]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if(na > 2) for( jj in 1:nb)for(kk in 1:nc){
cl<-1-alpha
if(rep==0) MSw<-MSe
if(rep==1) MSw<-MSpool.c
if(rep>=2) {
dfpool.a.b.c.bc<-dfe.a+dfe.b+dfe.c+dfe.bc
MSpool.a.b.c.bc<-(SSe.a+SSe.b+SSe.c+SSe.bc)/dfpool.a.b.c.bc
MSw<-MSpool.a.b.c.bc
}
mx<-meanxijk[,jj,kk]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-(qtukey(cl,k,dfe.ab)*MSe.ab+qtukey(cl,k,dfe.c)*MSe.c)/(MSe.ab+MSe.c)
if(rep==2) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.b)*MSe.b*(nb-1)+
qtukey(cl,k,dfe.c)*MSe.c*(nc-1)+qtukey(cl,k,dfe.bc)*MSe.bc*(nb-1)*(nc-1))/
(MSe.a+MSe.b*(nb-1)+MSe.c*(nc-1)+MSe.bc*(nb-1)*(nc-1))
if(rep==3) q2<-(qtukey(cl,k,dfe.a)*MSe.a+qtukey(cl,k,dfe.ab)*MSe.ab*(nb-1)+
qtukey(cl,k,dfe.ac)*MSe.ac*(nc-1)+qtukey(cl,k,dfe.abc)*MSe.abc*(nb-1)*(nc-1))/
(MSe.a+MSe.ab*(nb-1)+MSe.ac*(nc-1)+MSe.abc*(nb-1)*(nc-1))
mx2<-sort(mx);nijk2<-nijk;nijk2[,jj,kk]<-nijk[sort.list(mx),jj,kk]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-(qtukey(cl,m,dfe.ab)*MSe.ab+qtukey(cl,m,dfe.c)*MSe.c)/(MSe.ab+MSe.c)
if(rep==2) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.b)*MSe.b*(nb-1)+
qtukey(cl,m,dfe.c)*MSe.c*(nc-1)+qtukey(cl,m,dfe.bc)*MSe.bc*(nb-1)*(nc-1))/
(MSe.a+MSe.b*(nb-1)+MSe.c*(nc-1)+MSe.bc*(nb-1)*(nc-1))
if(rep==3) q3<-(qtukey(cl,m,dfe.a)*MSe.a+qtukey(cl,m,dfe.ab)*MSe.ab*(nb-1)+
qtukey(cl,m,dfe.ac)*MSe.ac*(nc-1)+qtukey(cl,m,dfe.abc)*MSe.abc*(nb-1)*(nc-1))/
(MSe.a+MSe.ab*(nb-1)+MSe.ac*(nc-1)+MSe.abc*(nb-1)*(nc-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nijk2[i,jj,kk];n22<-nijk2[j,jj,kk]
nhm<-2/(1/n12+1/n22)
WSD<-qvalue*sqrt(MSw/nhm)
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm)
}
if(opt==1) result2<-data.frame(Bj=jj,Ck=kk,meanAu=result[,1],meanAv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
if(opt==2) result2<-data.frame(Bj=jj,Ck=kk,meanAu=result[,1],meanAv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nAu=result[,9],nAv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
}
print(result2)
# write.table(result2,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if( na!=1 && nb > 2) {print.noquote("Test of Multiple Comparison of Factor B[AiCk]")
# write.table(c("Test of Multiple Comparison of Factor B[AiCk]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if( na!=1 && nb > 2) for( ii in 1:na)for(kk in 1:nc){
cl<-1-alpha
if(rep==0) MSw<-MSe
if(rep==1) MSw<-MSpool.c
if(rep>=2) MSw<-MSpool.b.bc
mx<-meanxijk[ii,,kk]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-(qtukey(cl,k,dfe.c)*MSe.c+qtukey(cl,k,dfe.ab)*MSe.ab)/(MSe.ab+MSe.c)
if(rep==2) q2<-(qtukey(cl,k,dfe.b)*MSe.b+qtukey(cl,k,dfe.bc)*MSe.bc*(nc-1))/(MSe.b+MSe.bc*(nc-1))
if(rep==3) q2<-(qtukey(cl,k,dfe.b)*MSe.b+qtukey(cl,k,dfe.ab)*MSe.ab*(na-1)+
qtukey(cl,k,dfe.bc)*MSe.bc*(nc-1)+qtukey(cl,k,dfe.abc)*MSe.abc*(na-1)*(nc-1))/
(MSe.b+MSe.ab*(na-1)+MSe.bc*(nc-1)+MSe.abc*(na-1)*(nc-1))
mx2<-sort(mx);nijk2<-nijk;nijk2[ii,,kk]<-nijk[ii,sort.list(mx),kk]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-(qtukey(cl,m,dfe.c)*MSe.c+qtukey(cl,m,dfe.ab)*MSe.ab)/(MSe.ab+MSe.c)
if(rep==2) q3<-(qtukey(cl,m,dfe.b)*MSe.b+qtukey(cl,m,dfe.bc)*MSe.bc*(nc-1))/
(MSe.b+MSe.bc*(nc-1));
if(rep==3) q3<-(qtukey(cl,m,dfe.b)*MSe.b+qtukey(cl,m,dfe.ab)*MSe.ab*(na-1)+
qtukey(cl,m,dfe.bc)*MSe.bc*(nc-1)+qtukey(cl,m,dfe.abc)*MSe.abc*(na-1)*(nc-1))/
(MSe.b+MSe.ab*(na-1)+MSe.bc*(nc-1)+MSe.abc*(na-1)*(nc-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nijk2[ii,i,kk];n22<-nijk2[ii,j,kk]
nhm<-2/(1/n12+1/n22)
WSD<-qvalue*sqrt(MSw/nhm)
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm)
}
if(opt==1) result2<-data.frame(Ai=ii,Ck=kk,meanBu=result[,1],meanBv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nBu=result[,9],nBv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
if(opt==2) result2<-data.frame(Ai=ii,Ck=kk,meanBu=result[,1],meanBv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nBu=result[,9],nBv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
}
print(result2)
# write.table(result2,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if( na!=1 && nc > 2) {print.noquote("Test of Multiple Comparison of Factor C[AiBj]")
# write.table(c("Test of Multiple Comparison of Factor C[AiBj]"),file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
if( na!=1 && nc > 2) for( ii in 1:na)for(jj in 1:nb){
cl<-1-alpha
if(rep==0) MSw<-MSe
if(rep==1) MSw<-MSe.c
if(rep>=2) MSw<-MSpool.c.bc
mx<-meanxijk[ii,jj,]
k<-length(mx)
if(rep==0) q2<-qtukey(cl,k,dfe)
if(rep==1) q2<-qtukey(cl,k,dfe.c)
if(rep==2) q2<-(qtukey(cl,k,dfe.c)*MSe.c+qtukey(cl,k,dfe.bc)*MSe.bc*(nb-1))/(MSe.c+MSe.bc*(nb-1));
if(rep==3) q2<-(qtukey(cl,k,dfe.c)*MSe.c+qtukey(cl,k,dfe.ac)*MSe.ac*(na-1)+
qtukey(cl,k,dfe.bc)*MSe.bc*(nb-1)+qtukey(cl,k,dfe.abc)*MSe.abc*(na-1)*(nb-1))/
(MSe.c+MSe.ac*(na-1)+MSe.bc*(nb-1)+MSe.abc*(na-1)*(nb-1))
mx2<-sort(mx);nijk2<-nijk;nijk2[ii,jj,]<-nijk[ii,jj,sort.list(mx)]
k1<-k-1;cnt<-0;nn<-k*(k-1)/2;result<-matrix(0,nrow=nn,ncol=13)
d2<-0;cnt2<-0;m2<-0
for( i in 1:k1){
i1<-i+1
for( j in i1:k){
cnt2<-cnt2+1
m<-abs(i-j)+1
d<-abs(mx2[i]-mx2[j])
if(rep==0) q3<-qtukey(cl,m,dfe)
if(rep==1) q3<-qtukey(cl,m,dfe.c)
if(rep==2) q3<-(qtukey(cl,m,dfe.c)*MSe.c+qtukey(cl,m,dfe.bc)*MSe.bc*(nb-1))/
(MSe.c+MSe.bc*(nb-1));
if(rep==3) q3<-(qtukey(cl,m,dfe.c)*MSe.c+qtukey(cl,m,dfe.ac)*MSe.ac*(na-1)+
qtukey(cl,m,dfe.bc)*MSe.bc*(nb-1)+qtukey(cl,m,dfe.abc)*MSe.abc*(na-1)*(nb-1))/
(MSe.c+MSe.ac*(na-1)+MSe.bc*(nb-1)+MSe.abc*(na-1)*(nb-1))
if(opt==2) q3<-q2
qvalue<-(q2+q3)/2
n12<-nijk[ii,jj,i];n22<-nijk[ii,jj,j]
nhm<-2/(1/n12+1/n22)
WSD<-qvalue*sqrt(MSw/nhm)
if(d-WSD>0) sig<-1 else sig<-0
result[cnt2,]<-cbind(mx2[i],mx2[j],d,WSD,k,m,qvalue,MSw,n12,n22,df,sig,nhm)
}
if(opt==1) result2<-data.frame(Ai=ii,Bj=jj,meanCu=result[,1],meanCv=result[,2],d=result[,3],
WSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
if(opt==2) result2<-data.frame(Ai=ii,Bj=jj,meanCu=result[,1],meanCv=result[,2],d=result[,3],
HSD=result[,4],m=result[,5],h=result[,6],qvalue=result[,7],MSw=result[,8],
nCu=result[,9],nCv=result[,10],nhm=result[,13],dfw=result[,11],sig005=result[,12])
}
print(result2)
# write.table(result2,file="mct.xls",row.names=F,col.names=F,
# append=T,quote=F,sep="\t")
}
}
}
if(n2!=n3) print("input error. The column size of data should be equal to na*nb*nc")
}
######var2###############################################################
var2<-function(x){
if(is.vector(x)) n<-length(x)
if(is.matrix(x)) n<-nrow(x)
var2<-var(x)*(n-1)/n
var2
}
##########varimax2##########################################################
varimax2<-function(am){
m <- ncol(am)
n <- nrow(am)
commu <- apply(am^2,1,sum)
icomm <- 1/sqrt(commu)
h <- diag(icomm,nrow=length(commu))
aini <- h%*%am
b <- aini
m1 <- m-1
kk <- 1000
for ( k in 1:kk){b2 <- b;
for ( p in 1:m1) {p1<- p+1;
for ( q in p1:m){
bb <- b^2; bpq1 <- bb[,p] - bb[,q];
bpq2 <- b[,p]*b[,q];
a1 <- sum(bpq1); b1 <- 2*sum(bpq2);
c1 <- sum(bpq1^2)-4*sum(bpq2^2);
d1 <- 4*sum(bpq1*bpq2);
v1 <- d1 - 2*a1*b1/n;
v2 <- c1 - (a1^2-b1^2)/n;
if ( v1 > 0 && v2 > 0) theta4 <- atan(v1/v2);
if ( v1 > 0 && v2 < 0) theta4 <- atan(v1/v2) + 3.14;
if ( v1 < 0 && v2 > 0) theta4 <- atan(v1/v2);
if ( v1 < 0 && v2 < 0) theta4 <- atan(v1/v2) - 3.14;
theta <- theta4/4;
costhe <- cos(theta);
sinthe <- sin(theta);
t1 <- diag(m);
t1[p,p] <- costhe;
t1[p,q] <- -sinthe;
t1[q,q] <- costhe;
t1[q,p] <- sinthe;
b <- b%*%t1;
}
}
bidf<- (b2 - b)^2;
dif <- sum(bidf)/(m*n);
if(dif < 0.001) {
kk<- k;
break
}
}
b<-solve(h)%*%b ;
eigen1<- apply(b^2,2,sum);
par(mfrow=c(1,1),pty="s",las=1,lwd=2,tck=0.02)
axis1<-range(c(b[,1],b[,2]))
axis2<-axis1
plot(axis1,axis2,type="n")
text(b[,1],b[,2])
write.table(c("varimax rotaion"),file="varimax2-table.xls",
row.name=F,col.names=F,append=T,quote=F,sep="\t")
write.table(round(b,3),file="varimax2-table.xls",
col.names=F,append=T,quote=F,sep="\t")
list(data=round(am,3),commu=round(commu,3),
loading=round(b,3),evalue=round(eigen1,2),iter=kk)
}
#end###############################################################