##### 7.7 ######

 setwd("G:/R springer/R CODE 7")

 source("create_design.r")
 source("unreplicated.r")
 source("t2p.r")



 k=3
 X<-create.design(k)
 X

 X<-create.design(k,Yates=TRUE)
 X

 set.seed(101)
 n<-2^k
 b<-rep(0,n)
 e<-rnorm(n)
 y<-X%*%b+e
 y


 source("unreplicated.r")
 source("t2p.r")
 set.seed(101)
 t<-unreplicated(y,X)
 t$beta

 table<-cbind(t(t$p.value),t(t$dec))
 colnames(table)<-c("p.value","Decision")
 table




 set.seed(101)
 U<-unreplicated(y,X,IER=0.2)
 table<-cbind(U$beta[-1],t(U$p.value),t(U$dec))
 colnames(table)=c("beta","p.value","Decision")
 table





 b[2:3]<-1
 y<-X%*%b+rnorm(8)
 set.seed(101) 
 U<-unreplicated(y,X)
 U

 table<-cbind(U$beta[-1],t(U$p.value),t(U$dec))
 colnames(table)=c("beta","p.value","Decision")
 table

 beta<-U$beta[-1]
 SSE<-sum(beta^2)
 T.oss<-beta^2/(SSE-beta^2) 
 names(T.oss)<-colnames(X)[-1]
 round(T.oss,digits=4)



 set.seed(10)
 y.star<-sample(y)
 beta.star<-unreplicated(y.star,X,3)$beta[-1]

 SSE.star<-sum(beta.star^2)
 T.star<-beta.star^2/(SSE.star-beta.star^2) 
 names(T.star)<-colnames(X)[-1]
 round(T.star,digits=4)


############## example p.202 (from Montgomey, 1991) ###########


 d<-read.csv("Mont_911.csv",header=TRUE)
 y<-c(d[,1],d[,2])
 X<-create.design(4,Yates=TRUE)
 A<-factor(rep(X[,2],2))
 B<-factor(rep(X[,3],2))
 C<-factor(rep(X[,4],2))
 D<-factor(rep(X[,5],2))
 summary(aov(y~A*B*C*D))



 set.seed(10)
 X<-create.design(4,Yates=TRUE)
 
 y = apply(d,1,mean)

 U<-unreplicated(y,X,step.up=TRUE,EER=0.003794,B=5000)	## critical p = 0.00379 for 10% EER (7.7.1 & see later on)
 table<-cbind(U$beta[-1],t(U$p.value),t(U$dec))
 colnames(table)=c("beta","p.value","Decision")
 table

 t<-qqnorm(table[,1],ylim=c(-0.2,0.2))
 qx<-t$x
 qy<-t$y
 text(qx,qy+0.01,labels=colnames(X)[-1],cex=.6 )

 dec<-table[,3]
 points(qx[dec==1],qy[dec==1],pch=20)

 ##### Calibration of the T_p test (based on Monte Carlo simulation) ######

 set.seed(101)
 k=4 ; X = create.design(4,Yates=TRUE)
 n=2^k
 G=10000
 F<-array(0,dim=c(G,n-1))
 P<-array(0,dim=c(G,n-1))
 for(cc in 1:G){
 b<-sort(abs(rnorm(n-1,sqrt(1/n))))
 for(j in 2:(n-1)){
 F[cc,j]<-(j-1)*b[j]^2/(sum(b[1:j]^2)-b[j]^2)
 }
 }
 for(j in 2:(n-1)){
 P[,j]<- 1-pf(F[,j],1,(j-1))
 }
 min<-apply(P[,-1],1,min)
 
 alpha<-c(0.01,0.05,0.1,0.2)

 EER <- quantile(min,alpha) 		# EER critical values

# cbind(alpha,EER)

 p<-dim(X)[2]
 alpha <- 0.1
 tilde.alpha<-1-(1-alpha)^(1/(p-2))
 EER.critical<-apply(P,2,function(x){quantile(x,tilde.alpha)})

#critical<-data.frame(step=seq(1,15),critical=critical)  # just for viewing

 EER.critical					

###  EER critical values for each step of the step-up version 
###  (you'll need to modify the body of function 'unreplicated.r')


 p<-dim(X)[2]
 alpha <- 0.2
 IER.critical<-apply(P,2,function(x){quantile(x,alpha)})

 IER.critical					

###  IER critical values for each step of the step-up version 
###  (you'll need to modify the body of function 'unreplicated.r')



#################################################################