### Create Rows for GEE simulations corresponding to Mancl and DeRouen's method
### First run 3 programs: sim.pois.exchange.fix.ssc, sim.pois.indep.nofix.ssc, sim.preisser.ssc

row1<-rep(NA,8)
nsim<-1000

### GEE Pois, Exchangeable, tau=0
source("out.sim.pois.exchange.fix")
### creates list called out, with results from pois exchange simulations

ftest.md<-function(beta,Vmd,K,p,alpha=.05){
	### create pvalues first
	x<- 1 - pf( (beta**2)/Vmd,1,K-p)
	prop.le.alpha<-length(x[x<=alpha])/length(x)
   prop.le.alpha
}

mean.ci.length.md<-function(Vmd,K,p,alpha=.05){
   ci.length<-2*sqrt(qf(1-alpha,1,K-p))*sqrt(Vmd)
   mean(ci.length)
}





### GEE Pois, Exchangeable, tau=0

index<-rep(T,nsim)
if (!is.null(out$erri.d1)) index[out$erri.d1]<-F 
###number of converging simulated data sets
length(out$bd1[index])
row1[1]<- ftest.md(out$bd1[index],out$vd1[index,3],20,2)
row1[2]<-mean.ci.length.md(out$vd1[index,3],20,2)

### GEE Pois, Exchangeable, tau=1/2

index<-rep(T,nsim)
if (!is.null(out$erri.d4)) index[out$erri.d4]<-F 
###number of converging simulated data sets
length(out$bd4[index])
row1[3]<- ftest.md(out$bd4[index],out$vd4[index,3],20,2)
row1[4]<-mean.ci.length.md(out$vd4[index,3],20,2)


### GEE Pois, Independent, tau=0
source("out.sim.pois.indep.nofix")
out<-out.sim.pois.indep.nofix

index<-rep(T,nsim)
if (!is.null(out$erri.d1)) index[out$erri.d1]<-F 
###number of converging simulated data sets
length(out$bd1[index])
row1[5]<- ftest.md(out$bd1[index],out$vd1[index,3],20,2)
row1[6]<-mean.ci.length.md(out$vd1[index,3],20,2)

### GEE Pois, Independent, tau=1/2

index<-rep(T,nsim)
if (!is.null(out$erri.d4)) index[out$erri.d4]<-F 
###number of converging simulated data sets
length(out$bd4[index])
row1[7]<- ftest.md(out$bd4[index],out$vd4[index,3],20,2)
row1[8]<-mean.ci.length.md(out$vd4[index,3],20,2)



### Preisser data (GUIDE data) simulation results
row2<-rep(NA,10)
source("out.sim.preisser")
out<-out.sim.preisser

index<-rep(T,nsim)
if (!is.null(out$erri2)) index[out$erri2]<-F 
###number of converging simulated data sets
length(out$b2[index])
row2[1]<- ftest.md(out$b2[index],out$v2[index,3],38,6)
row2[2]<-mean.ci.length.md(out$v2[index,3],38,6)

index<-rep(T,nsim)
if (!is.null(out$erri3)) index[out$erri3]<-F 
###number of converging simulated data sets
length(out$b3[index])
row2[3]<- ftest.md(out$b3[index],out$v3[index,3],38,6)
row2[4]<-mean.ci.length.md(out$v3[index,3],38,6)

index<-rep(T,nsim)
if (!is.null(out$erri4)) index[out$erri4]<-F 
###number of converging simulated data sets
length(out$b4[index])
row2[5]<- ftest.md(out$b4[index],out$v4[index,3],38,6)
row2[6]<-mean.ci.length.md(out$v4[index,3],38,6)

index<-rep(T,nsim)
if (!is.null(out$erri5)) index[out$erri5]<-F 
###number of converging simulated data sets
length(out$b5[index])
row2[7]<- ftest.md(out$b5[index],out$v5[index,3],38,6)
row2[8]<-mean.ci.length.md(out$v5[index,3],38,6)

index<-rep(T,nsim)
if (!is.null(out$erri6)) index[out$erri6]<-F 
###number of converging simulated data sets
length(out$b6[index])
row2[9]<- ftest.md(out$b6[index],out$v6[index,3],38,6)
row2[10]<-mean.ci.length.md(out$v6[index,3],38,6)


round(row1,digits=3)
round(row2,digits=3)