
## run BMEM for mediation analysis

## generate batch file

gen.batch<-function(datafile){
 filename<-'batch.txt'
 cat("out.txt\n", file=filename)
 cat(datafile, file=filename, append=T)
 cat("\n", file=filename, append=T)
 cat("99999\n", file=filename, append=T)
 cat(".95\n", file=filename, append=T)  
 cat(runif(1), file=filename, append=T)
 cat("\n", file=filename, append=T)
 cat("1000\n", file=filename, append=T)
}


runone<-function(datafile, resfile){
    for (i in 1001:2000){
    	  file<-paste(datafile,'-',i,'.txt', sep='')
       gen.batch(file)
       system("./bmem < batch.txt")
       temp<-scan('out.txt')
       cat(c(i, temp), file=resfile, append=T)
       cat("\n", file=resfile,append=T)
    }
}




for (ab in c(0, .14, .39, .59)){
	for (N in c(50,100,200)){
		for (miss in c('MCAR','MNAR1', 'MNAR2')){
			for (rate in c(.1,.2,.4) ){
				for (av in c('av','noav')){
				## folder name
				data.name<-paste('/pscratch/zzhang4/Private/mediation/data/ab',ab,'N',N,miss,'M',rate
,'/data-',av,sep='')
				res.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'-em-results.txt',sep='')
				runone(data.name, res.name)              
				}
			}
			}
			
		for (miss in 'MAR'){
			for (rate in c(1.5, 4, 9) ){
			    for (av in c('av','noav')){
				## folder name
				data.name<-paste('/pscratch/zzhang4/Private/mediation/data/ab',ab,'N',N,miss,'M',rate
,'/data-',av,sep='')
				res.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'-em-results.txt',sep='')
				runone(data.name, res.name)              
				}              
				}
			}	
		}
	
	}

==== Process the results ====


runone<-function(filename, par){
    temp<-read.table(filename, fill=T)
    if (dim(temp)[2]==1){
    c(0,0,0)
    }else{
    par<-par*par
    if (par==0){
      bias<-100*(mean(temp[,2],na.rm=T)-par)
    }else{
      bias<-100*(mean(temp[,2],na.rm=T)-par)/par
    }
    n<-dim(temp)[1]
    cvg<-mean((temp[,6]par), na.rm=T)*100
    power<-mean((temp[,6]>0) | (temp[,7]<0), na.rm=T)*100
    c(bias, cvg, power)
    }
}




for (ab in c(0, .14, .39, .59)){
	for (N in c(50,100,200)){
		for (miss in c('MCAR','MNAR1', 'MNAR2')){
			for (rate in c(.1,.2,.4) ){
				for (av in c('av','noav')){
		res.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'-em-results.txt',sep='')
		if (file.exists(res.name)){
				cat(ab, N, miss, rate, av, runone(res.name, ab),'\n')    }          
				}
			}
			}
			
		for (miss in 'MAR'){
			for (rate in c(1.5, 4, 9) ){
			    for (av in c('av','noav')){
		res.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'-em-results.txt',sep='')
				if (file.exists(res.name)){
				cat(ab, N, miss, rate, av, runone(res.name, ab),'\n')    }  
				
				}              
				}
			}	
		}
	
	}