##### 6.7.1 ######

setwd("C:/R springer/R CODE 6")

source("USP.r")
source("CSP.r")
source("IC_CSP.r")
source("IC_USP.r")
source("synchro_summary.r")
source("plot_hsu.r")


set.seed(1000)
I <-3
J<-2
n<-2
y<-rnorm((n*I*J))
x1<-factor(rep(seq(1,I),each=(n*J)))
x2<-factor(rep(rep(c(1,2),each=n),I))
x<-data.frame(A=x1,B=x2)
x
y


t=USP(y,x)
str(t)
synchro.summary(t)


 attach(x)
 fit<-aov(y~A*B)
 summary(fit)


 set.seed(10)
 n<-4
 x1<-factor(rep(seq(1,I),each=(n*J)))
 x2<-factor(rep(rep(c(1,2),each=n),I))
 x<-data.frame(A=x1,B=x2)
 y<-rnorm(n*I*J)
 alpha=rep(c(-1,0,1),each=2*n)
 alpha.beta<-c(rep(c(-1,1),each=n),rep(0,each=(2*n)),rep(c(1,-1),each=n))/5
 y<-y+alpha+alpha.beta
 round(y,digits=5)


 t<-CSP(y,x,exact=TRUE)
 synchro.summary(t)

 t<-CSP(y,x,C=5000)
 synchro.summary(t)



  set.seed(101)
  n.comp<-choose(I,2)
  IC.A<-array(0,dim=c(n.comp,2))
  F<-array(0,dim=c(n.comp,2))
  m<-array(0,dim=c(n.comp,2))
  k<-1
  for(i in 1:(I-1)){
  for(s in (i+1):I){
  F[k,1]<-paste("m",i,sep="")
  F[k,2]<-paste("m",s,sep="")

  IC.A[k,]<-IC.USP(y,x,i,s,C = 200)

  m[k,1]<-mean(y[x[,1]==i])
  m[k,2]<-mean(y[x[,1]==s])
  k=k+1
  }
  }

 out.A<-data.frame(F1=F[,1],F2=F[,2],m1=m[,1],m2=m[,2],MSD=m[,1]-m[,2]-IC.A[,1])
 out.A

 colnames(IC.A) = c("lower","upper")
 rownames(IC.A) = c("1-2","1-3","2-3")
 IC.A




 set.seed(1)
 m1<-mean(y[x[,2]==1])
 m2<-mean(y[x[,2]==2])

 IC.B<-IC.USP(y,x,1,2,row.perm=TRUE,C = 200)

 MSD=m1-m2-IC.B[1]
 out.B<-data.frame(F1="m1",F2="m2",m1=m1,m2=m2,MSD=MSD)
 out.B


### The files "hsu-ICA.eps" and "hsu-ICB.eps" will be on your working directory ### 

plot.hsu(out.A,file="ICA",main="A",title="USP ICs",dev="eps",measure=expression(sigma))
plot.hsu(out.B,file="ICB",main="B",title="USP ICs",dev="eps",measure=expression(sigma))



###########example p. 165#######################################

data<-data.frame(y=y,A=x[,1],B=x[,2])
A<-as.factor(data$A)
B<-as.factor(data$B)
fit<-aov(y ~ A*B)
TukeyHSD(fit,"A")


########### example IC.CSP #######################################

IC.A = array(0,dim=c(choose(I,2),2))

cont=1
 for(i in 1:(I-1)){
 for(s in (i+1):I){
IC.A[cont,]<-IC.CSP(y,x,i,s,row.perm=FALSE,conf.level=(1-2/35))
cont=cont+1
}
}


 colnames(IC.A) = c("lower","upper")
 rownames(IC.A) = c("1-2","1-3","2-3")
 IC.A