.Random.seed<-c(47,42,43,4,43,3,12,38,15,33,6,2)

nsim<-1000


library(yags)
source("myags.ssc")
source("myags.pmat.fit.mancl.ssc")
source("saws2final.ssc")
source("LZ.exchalp.ssc")
CORSTR<-"independence"
#CORSTR.fit<-"exchangeable"
CORSTR.fit<-"independence"
#CORSTR<-"exchangeable"



n<-137
pr<-read.table("preisser.txt",header=T)
pr[1,]
y<-pr$bothered
female<-pr$female
age<-pr$age
dayacc<-pr$dayacc
severe<-pr$severe
toilet<-pr$toilet
pract.id<-pr$pract.id
x<-matrix(c(rep(1,n),female,age,dayacc,severe,toilet),n,6)

expit<-function(eta){ exp(eta)/(1+exp(eta)) }

#mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr="exchangeable")
#sout<-saws2(mout$coef,mout$u,mout$omega)
#sout

## fit the model without toilet=6th column of design matrix
beta2<-yags(y~age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR.fit)$coef
fit2<-expit(x[,-2] %*% matrix(beta2,5,1) )
beta3<-yags(y~female+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR.fit)$coef
fit3<-expit(x[,-3] %*% matrix(beta3,5,1) )
beta4<-yags(y~female+age+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR.fit)$coef
fit4<-expit(x[,-4] %*% matrix(beta4,5,1) )
beta5<-yags(y~female+age+dayacc+toilet,id=pract.id,family=binomial,corstr=CORSTR.fit)$coef
fit5<-expit(x[,-5] %*% matrix(beta5,5,1) )
beta6<-yags(y~female+age+dayacc+severe,id=pract.id,family=binomial,corstr=CORSTR.fit)$coef
fit6<-expit(x[,-6] %*% matrix(beta6,5,1) )


now<-proc.time()
 numM<-8

 p2<-p3<-p4<-p5<-p6<-matrix(NA,nsim,numM)
 ci2<-ci3<-ci4<-ci5<-ci6<-matrix(NA,nsim,numM)
 nreject2<-nreject3<-nreject4<-nreject5<-nreject6<-rep(0,numM)
 errors2<-errors3<-errors4<-errors5<-errors6<-NULL
 erri2<-erri3<-erri4<-erri5<-erri6<-NULL 
 cnt2<-cnt3<-cnt4<-cnt5<-cnt6<-0
 uf2<-uf3<-uf4<-uf5<-uf6<-rep(NA,nsim)
 b2<-b3<-b4<-b5<-b6<-rep(NA,nsim)
 w2<-w3<-w4<-w5<-w6<-rep(NA,nsim)
 v2<-v3<-v4<-v5<-v6<-matrix(NA,nsim,4,dimnames=list(NULL,c("Vm","Vs","VsMD","Va")))

 for (i in 1:nsim){


print(paste("i=",i))
	

y<-rbinom(n,rep(1,n),fit2)
mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR)
b2[i]<-mout$coefficients[2]
v2[i,]<-c(mout$naive.parmvar[2,2],mout$robust.parmvar[2,2],mout$robust.parmvar.md[2,2],NA)	
	if (is.null(mout$errorcode)){	
sout<-saws2(mout$coef,mout$u,mout$omega,test=2)
   v2[i,4]<-sout$arse1^2
   w2[i]<-sout$ssew
	p2[i,]<-sout$twosided.pvalues[1:numM]
   ci2[i,]<-sout$ci.length[1:numM]
   add<-rep(0,numM)
   add[p2[i,]<=0.05]<-1
   nreject2<-nreject2+ add
   cnt2<-cnt2+1
   uf2[i]<-mout$ufactor
print(round(nreject2/cnt2,digits=4))
    }
else{ errors2<-c(errors2,mout$errorcode) 
	   erri2<-c(erri2,i) 
	   uf2[i]<-mout$ufactor

}


y<-rbinom(n,rep(1,n),fit3)
dat1.3<<-data.frame(ysim.fit3=y,female=female,age=age,dayacc=dayacc,severe=severe,toilet=toilet,pract.id=pract.id)

mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR)
b3[i]<-mout$coefficients[3]
v3[i,]<-c(mout$naive.parmvar[3,3],mout$robust.parmvar[3,3],mout$robust.parmvar.md[3,3],NA)	

	if (is.null(mout$errorcode)){	
sout<-saws2(mout$coef,mout$u,mout$omega,test=3)
   v3[i,4]<-sout$arse1^2
   w3[i]<-sout$ssew
	p3[i,]<-sout$twosided.pvalues[1:numM]
   ci3[i,]<-sout$ci.length[1:numM]
   add<-rep(0,numM)
   add[p3[i,]<=0.05]<-1
   nreject3<-nreject3+ add
   cnt3<-cnt3+1
   uf3[i]<-mout$ufactor

print(round(nreject3/cnt3,digits=4))
    }
else{ errors3<-c(errors3,mout$errorcode) 
	   erri3<-c(erri3,i) 
	   uf3[i]<-mout$ufactor

}


y<-rbinom(n,rep(1,n),fit4)

mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR)
b4[i]<-mout$coefficients[4]
v4[i,]<-c(mout$naive.parmvar[4,4],mout$robust.parmvar[4,4],mout$robust.parmvar.md[4,4],NA)	

	if (is.null(mout$errorcode)){	
sout<-saws2(mout$coef,mout$u,mout$omega,test=4)
   v4[i,4]<-sout$arse1^2
   w4[i]<-sout$ssew
	p4[i,]<-sout$twosided.pvalues[1:numM]
   ci4[i,]<-sout$ci.length[1:numM]
   add<-rep(0,numM)
   add[p4[i,]<=0.05]<-1
   nreject4<-nreject4+ add
   cnt4<-cnt4+1
   uf4[i]<-mout$ufactor
print(round(nreject4/cnt4,digits=4))
    }
else{ errors4<-c(errors4,mout$errorcode) 
	   erri4<-c(erri4,i) 
   uf4[i]<-mout$ufactor
}


y<-rbinom(n,rep(1,n),fit5)

mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR)
b5[i]<-mout$coefficients[5]
v5[i,]<-c(mout$naive.parmvar[5,5],mout$robust.parmvar[5,5],mout$robust.parmvar.md[5,5],NA)	

	if (is.null(mout$errorcode)){	
sout<-saws2(mout$coef,mout$u,mout$omega,test=5)
   v5[i,4]<-sout$arse1^2
   w5[i]<-sout$ssew
	p5[i,]<-sout$twosided.pvalues[1:numM]
   ci5[i,]<-sout$ci.length[1:numM]
   add<-rep(0,numM)
   add[p5[i,]<=0.05]<-1
   nreject5<-nreject5+ add
   cnt5<-cnt5+1
   uf5[i]<-mout$ufactor
print(round(nreject5/cnt5,digits=4))
    }
else{ errors5<-c(errors5,mout$errorcode) 
	   erri5<-c(erri5,i) 
   uf5[i]<-mout$ufactor
}


y<-rbinom(n,rep(1,n),fit6)

mout<-myags(y~female+age+dayacc+severe+toilet,id=pract.id,family=binomial,corstr=CORSTR)
b6[i]<-mout$coefficients[6]
v6[i,]<-c(mout$naive.parmvar[6,6],mout$robust.parmvar[6,6],mout$robust.parmvar.md[6,6],NA)	

	if (is.null(mout$errorcode)){	
sout<-saws2(mout$coef,mout$u,mout$omega,test=6)
   v6[i,4]<-sout$arse1^2
   w6[i]<-sout$ssew
	p6[i,]<-sout$twosided.pvalues[1:numM]
   ci6[i,]<-sout$ci.length[1:numM]
   add<-rep(0,numM)
   add[p6[i,]<=0.05]<-1
   nreject6<-nreject6+ add
   cnt6<-cnt6+1
   uf6[i]<-mout$ufactor
print(round(nreject6/cnt6,digits=4))
    }
else{ errors6<-c(errors6,mout$errorcode) 
	   erri6<-c(erri6,i) 
   uf6[i]<-mout$ufactor
}




}

out<-list(
p2=p2,ci2=ci2,
v2=v2,b2=b2,w2=w2,
errors2=errors2,erri2=erri2,uf2=uf2,
p3=p3,ci3=ci3,
v3=v3,b3=b3,w3=w3,
errors3=errors3,erri3=erri3,uf3=uf3,
p4=p4,ci4=ci4,
v4=v4,b4=b4,w4=w4,
errors4=errors4,erri4=erri4,uf4=uf4,
p5=p5,ci5=ci5,
v5=v5,b5=b5,w5=w5,
errors5=errors5,erri5=erri5,uf5=uf5,
p6=p6,ci6=ci6,
v6=v6,b6=b6,w6=w6,
errors6=errors6,erri6=erri6,uf6=uf6)

#dump("out","sim.EP1.out")
#dump("out","sim.EP0.out")
#dump("out","sim.preis8.17.00.out")
out
now-proc.time()



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

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)
if (!is.null(outerri)) index[outerri]<-F

print(paste("Number Converged=",length((1:nsim)[index]) ) )

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()
}

out<-out.sim.preisser

quantiles(out$w2)
print.sim(out$p2,out$ci2,out$erri2)
quantiles(out$w3)
print.sim(out$p3,out$ci3,out$erri3)
quantiles(out$w4)
print.sim(out$p4,out$ci4,out$erri4)
quantiles(out$w5)
print.sim(out$p5,out$ci5,out$erri5)
quantiles(out$w6)
print.sim(out$p6,out$ci6,out$erri6)


mult<-1
quantiles(out$w2)
quantiles(mult*out$v2[,1])
quantiles(mult*out$v2[,2])
quantiles(mult*out$v2[,3])
quantiles(mult*out$v2[,4])
mult*var(out$b2)

mult<-1
quantiles(mult*out$v3[,1])
quantiles(mult*out$v3[,2])
quantiles(mult*out$v3[,3])
quantiles(mult*out$v3[,4])
mult*var(out$b3)

mult<-1
quantiles(mult*out$v4[,1])
quantiles(mult*out$v4[,2])
quantiles(mult*out$v4[,3])
quantiles(mult*out$v4[,4])
mult*var(out$b4)

mult<-1
quantiles(mult*out$v5[,1])
quantiles(mult*out$v5[,2])
quantiles(mult*out$v5[,3])
quantiles(mult*out$v5[,4])
mult*var(out$b5)

mult<-1
quantiles(mult*out$v6[,1])
quantiles(mult*out$v6[,2])
quantiles(mult*out$v6[,3])
quantiles(mult*out$v6[,4])
mult*var(out$b6)

out.sim.preisser<-out
dump("out.sim.preisser","out.sim.preisser")
names(out)
proc.time()-now