## perform the simulations of Lin and Wei, 1989, JASA 1074-1078
## to start just consider models 10 (model-based does poor) and 11 (robust poor)
source("mcoxph.ssc")
source("saws2final.ssc")

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

nsim<-2
n<-50	
 bound<-function(x,lower,upper){
	x[x<lower]<-lower
	x[x>upper]<-upper
	x
}	
numM<-8
 p10<-p11<-p0<-matrix(NA,nsim,numM)
 ci10<-ci11<-ci0<-matrix(NA,nsim,numM)
 nreject10<-nreject11<-nreject0<-rep(0,numM)
 
 for (i in 1:nsim){
	z1<-rnorm(n)
	z1<-bound(z1,-5,5)
	z2<-rnorm(n)
	z2<-bound(z2,-5,5)
	phi<-rnorm(n,sd=.5)
	epsilon<-rexp(n)
	t0<- exp(z2)*epsilon
	t11<- exp(-.5*z2) + epsilon
	t10<- exp( -.5*z2 -z1^2 + phi )
	status<-rep(1,n)
	x<-matrix(c(z1,z2),ncol=2)	
print(paste("i=",i))
	cout<-mcoxph(Surv(t0,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
#print(sout)
   p0[i,]<-sout$twosided.pvalues
   ci0[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p0[i,]<=0.05]<-1
   nreject0<-nreject0+ add
print(round(nreject0/i,digits=4))

	
	cout<-mcoxph(Surv(t10,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
	#print(sout)
   p10[i,]<-sout$twosided.pvalues
   ci10[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p10[i,]<=0.05]<-1

   nreject10<-nreject10+ add

print(round(nreject10/i,digits=4))

	cout<-mcoxph(Surv(t11,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
	#print(sout)
   p11[i,]<-sout$twosided.pvalues
   ci11[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p11[i,]<=0.05]<-1
   nreject11<-nreject11+ add
print(round(nreject11/i,digits=4))


}


out50<-list(p0=p0,p10=p10,p11=p11,
    ci0=ci0,ci10=ci10,ci11=ci11)
proc.time()-now









n<-20	
 bound<-function(x,lower,upper){
	x[x<lower]<-lower
	x[x>upper]<-upper
	x
}	
 p10<-p11<-p0<-matrix(NA,nsim,numM)
 ci10<-ci11<-ci0<-matrix(NA,nsim,numM)
 nreject10<-nreject11<-nreject0<-rep(0,numM)
 
 for (i in 1:nsim){
	z1<-rnorm(n)
	z1<-bound(z1,-5,5)
	z2<-rnorm(n)
	z2<-bound(z2,-5,5)
	phi<-rnorm(n,sd=.5)
	epsilon<-rexp(n)
	t0<- exp(z2)*epsilon
	t11<- exp(-.5*z2) + epsilon
	t10<- exp( -.5*z2 -z1^2 + phi )
	status<-rep(1,n)
	x<-matrix(c(z1,z2),ncol=2)	
print(paste("i=",i))
	cout<-mcoxph(Surv(t0,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
	#print(sout)
   p0[i,]<-sout$twosided.pvalues
   ci0[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p0[i,]<=0.05]<-1
   nreject0<-nreject0+ add
print(round(nreject0/i,digits=4))

	
	cout<-mcoxph(Surv(t10,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
#print(sout)  
 p10[i,]<-sout$twosided.pvalues
   ci10[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p10[i,]<=0.05]<-1

   nreject10<-nreject10+ add

print(round(nreject10/i,digits=4))

	cout<-mcoxph(Surv(t11,status)~x)
	sout<-saws2(cout$beta,cout$u,cout$omega,test=1)
#print(sout) 
  p11[i,]<-sout$twosided.pvalues
   ci11[i,]<-sout$ci.length
   add<-rep(0,numM)
   add[p11[i,]<=0.05]<-1
   nreject11<-nreject11+ add
print(round(nreject11/i,digits=4))


}


out20<-list(p0=p0,p10=p10,p11=p11,
    ci0=ci0,ci10=ci10,ci11=ci11)
proc.time()-now



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(out50$p10,out50$ci10)
print.sim(out50$p11,out50$ci11)
print.sim(out50$p0,out50$ci0)

print.sim(out20$p10,out20$ci10)
print.sim(out20$p11,out20$ci11)
print.sim(out20$p0,out20$ci0)