now<-proc.time()
.Random.seed<-c(11,22,10,33,13,3,61,59,10,2,25,1)

source("clogist.calc.ssc")
source("clogist.ssc")
source("saws2final.ssc")



nsim<-1000





make.fayetal.data<-function(design,case,K=40,m=30){
	if (design=="ABC"){
			 z1<-rep(c(.5,.5,0,1,0,1,0,1),K/4)
			 z2<-rep(c(0,1,.5,.5,0,1,1,0),K/4)
			}
	else if (design=="D"){
			 z1<-c(rep(c(.5,.5),37),.25,.75,.25,.75,.25,.75)
			 z2<-c(rep(c(0,1),37),.5,.5,0,1,1,0)
			 K<-40
			}
	if (case==1) { alpha<-0 ; siga<-1; beta1<-0;sig1<-0;beta2<-0;sig2<-0 }
	else if (case==2) { alpha<--3.3 ; siga<-1; beta1<-0;sig1<-0;beta2<-0;sig2<-0 }
	else if (case==3) { alpha<-0 ; siga<-1; beta1<-0;sig1<-.5;beta2<-0;sig2<-.5 }
	else if (case==4) { alpha<-0 ; siga<-1; beta1<-0;sig1<-1;beta2<-0;sig2<-2 }
	else if (case==5) { alpha<--3.3 ; siga<-1; beta1<-4.2;sig1<-.5;beta2<-2.4;sig2<-.5 }
	else if (case==6) { alpha<--3.3 ; siga<-1; beta1<-4.2;sig1<-1;beta2<-2.4;sig2<-2 }

   a<-rnorm(K,mean=alpha,sd=sqrt(siga))
      
 if (sig1==0) b1<-rep(beta1,K)
    else b1<-rnorm(K,mean=beta1,sd=sqrt(sig1))
    if (sig2==0) b2<-rep(beta2,K)
    else b2<-rnorm(K,mean=beta2,sd=sqrt(sig2))

    expit<-function(eta) (1+exp(-eta))^(-1)
   y<-rep(0,2*K)
   strata<-rep(0,2*K)

	for (i in 1:K){
		       strata[(2*i-1):(2*i)]<-c(i,i)
				y[(2*i-1):(2*i)]<- rbinom(2,c(m,m),
				expit(a[i] + z1[(2*i-1):(2*i)]*b1[i] + z2[(2*i-1):(2*i)]*b2[i] ) )
	}


data.frame(z1=z1,z2=z2,Y=y/m,strata=strata,N=rep(m,2*K))

}


analyze.fayetal.data<-function(data){
			Y<-data[,"Y"]
		strata<-data[,"strata"]
		z1<-data[,"z1"]				
		z2<-data[,"z2"]
		N<-data[,"N"]
		x<-matrix(c(z1,z2),ncol=2)
	mout<-clogist.calc(N,round(Y*N),x,strata)
 sout<-saws2(mout$beta,mout$u,mout$omega,test=1)
     sout

}





### Do Design A Cases 1 and 4
### Design D Cases 1 and 4

NT<-8
 pa1<-pa4<-pd1<-pd4<-matrix(NA,nsim,NT)
 cia1<-cia4<-cid1<-cid4<-matrix(NA,nsim,NT)
 nrejecta1<-nrejecta4<-nrejectd1<-nrejectd4<-rep(0,NT)
 
 for (i in 1:nsim){


print(paste("i=",i))
	
	### Design A Case 1
	
dat<-make.fayetal.data("ABC",1,K=20)
sout<-analyze.fayetal.data(dat)
	
   pa1[i,]<-sout$twosided.pvalues
   cia1[i,]<-sout$ci.length
   add<-rep(0,NT)
   add[pa1[i,]<=0.05]<-1
   nrejecta1<-nrejecta1+ add
print(round(nrejecta1/i,digits=4))


	### Design A Case 4
	
dat<-make.fayetal.data("ABC",4,K=20)
sout<-analyze.fayetal.data(dat)
	
   pa4[i,]<-sout$twosided.pvalues
   cia4[i,]<-sout$ci.length
   add<-rep(0,NT)
   add[pa4[i,]<=0.05]<-1
   nrejecta4<-nrejecta4+ add
print(round(nrejecta4/i,digits=4))


	### Design D Case 1
	
dat<-make.fayetal.data("D",1,K=40)
sout<-analyze.fayetal.data(dat)
	
   pd1[i,]<-sout$twosided.pvalues
   cid1[i,]<-sout$ci.length
   add<-rep(0,NT)
   add[pd1[i,]<=0.05]<-1
   nrejectd1<-nrejectd1+ add
print(round(nrejectd1/i,digits=4))


	### Design D Case 4
	
dat<-make.fayetal.data("D",4,K=40)
sout<-analyze.fayetal.data(dat)
	
   pd4[i,]<-sout$twosided.pvalues
   cid4[i,]<-sout$ci.length
   add<-rep(0,NT)
   add[pd4[i,]<=0.05]<-1
   nrejectd4<-nrejectd4+ add
print(round(nrejectd4/i,digits=4))




}

out<-list(
   pa1=pa1,pa4=pa4,pd1=pd1,pd4=pd4,
   cia1=cia1,cia4=cia4,cid1=cid1,cid4=cid4)


print.sim<-function(outp,outci){

test.names<-c("Vm,Chi2","Vs,Chi2","Vs,F,K-p","Vs,F,d.tilde.no.H","Vs,F,d.hat.no.H","Va,Chi2","Va,F,d.tilde.with.H",
"Va,F,d.hat.with.H")

nsim<-dim(outp)[[1]]
index<-rep(T,nsim)

prop.le.alpha<-function(x,alpha=.05) 
 length(x[x<=alpha])/length(x)

print("P-Values")
pval<-apply(outp[index,],2,prop.le.alpha)
names(pval)<-test.names
print(pval)
print("Mean ci")
meanci<-apply(outci[index,],2,mean)
names(meanci)<-test.names
print(meanci)
invisible()
}




print.sim(out$pa1,out$cia1)
print.sim(out$pa4,out$cia4)
print.sim(out$pd1,out$cid1)
print.sim(out$pd4,out$cid4)