User Tools

Site Tools


lab:zhang:r_codes_for_running_mi_simulation
## R codes for running simulation of pairwise deletion and listwise deletion
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
				folder.name<-paste('ab',ab,'N',N,miss,'M',rate,sep='')
				
				## change the folder name in SAS
				repNS<-paste("sed 's/simucond/",folder.name,"/g' mi.sas > temp.sas",sep='')
              system(repNS)
              repNS<-paste("sed 's/noav/",av,"/g' temp.sas > temp1.sas",sep='')
              system(repNS)
			  
					if (av=='av'){
					  repNS<-paste("sed 's/var00list/x m y u v/g' temp1.sas > temp2.sas",sep='')
					  system(repNS)
					}else{
					  repNS<-paste("sed 's/var00list/x m y/g' temp1.sas > temp2.sas",sep='')
					  system(repNS)
					}
              
              ## Copy mi.sas to a folder
              folder.name<-paste('ab',ab,'N',N,miss,'M',rate,av,sep='')
              system(paste('mkdir ', folder.name))
              file.copy('temp2.sas', paste(folder.name,'/mi.sas',sep=''))
              unlink('temp.sas')
					unlink('temp1.sas')
              unlink('temp2.sas')
              ## Copy the submision file
              file.copy('qsub.sh', paste(folder.name,'/', folder.name, '.sh',sep=''))
           
				}
			}
			}
			
		for (miss in 'MAR'){
			for (rate in c(1.5, 4, 9) ){
			for (av in c('av','noav')){
				## folder name
				folder.name<-paste('ab',ab,'N',N,miss,'M',rate,sep='')
				
				## change the folder name in SAS
				repNS<-paste("sed 's/simucond/",folder.name,"/g' mi.sas > temp.sas",sep='')
              system(repNS)
              repNS<-paste("sed 's/noav/",av,"/g' temp.sas > temp1.sas",sep='')
              system(repNS)
				
				if (av=='av'){
					repNS<-paste("sed 's/var00list/x m y u v/g' temp1.sas > temp2.sas",sep='')
					system(repNS)
				}else{
					repNS<-paste("sed 's/var00list/x m y/g' temp1.sas > temp2.sas",sep='')
					system(repNS)
				}
              
              ## Copy mi.sas to a folder
              folder.name<-paste('ab',ab,'N',N,miss,'M',rate,av,sep='')
              system(paste('mkdir ', folder.name))
              file.copy('temp2.sas', paste(folder.name,'/mi.sas',sep=''))
              unlink('temp.sas')
				unlink('temp1.sas')
              unlink('temp2.sas')
              ## Copy the submision file
              file.copy('qsub.sh', paste(folder.name,'/', folder.name, '.sh',sep=''))           
				}          
				}
			}	
		}
	
	}
	
	
## Generate submission script

subname<-'submitall.sh'

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<-paste('cd /pscratch/zzhang4/Private/mediation/analysis/mi/ab',ab,'N',N,miss,'M',rate,av,'/\n qsub *.sh \n',sep='')
              cat(folder.name, file=subname, append=T)
           
				}
			}
			}
			
		for (miss in 'MAR'){
			for (rate in c(1.5, 4, 9) ){
			for (av in c('av','noav')){
				folder.name<-paste('cd /pscratch/zzhang4/Private/mediation/analysis/mi/ab',ab,'N',N,miss,'M',rate,av,'/\n qsub *.sh \n',sep='')
              cat(folder.name, file=subname, append=T)         
				}          
				}
			}	
		}
	
	}
a<-seq(958614,958900, by=1)
for (i in a){
  system(paste('qdel ', i))	
}



Process the results

runone<-function(filename, par){
    par<-par*par
    temp<-read.table(filename, fill=T,skip=1)
    
    if (par==0){
      bias<-100*(mean(temp[,1],na.rm=T)-par)
    }else{
      bias<-100*(mean(temp[,1],na.rm=T)-par)/par
    }
    n<-dim(temp)[1]
    cvg<-sum((temp[,3]<par) & (temp[,4]>par), na.rm=T)/sum(!is.na((temp[,3]<par) & (temp[,4]>par)))*100
    power<-sum((temp[,3]>0) | (temp[,4]<0), na.rm=T)/sum(!is.na((temp[,3]<par) & (temp[,4]>par)))*100
    c(bias, cvg, power, n)
}


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
				newname<-paste('ab',ab,'N',N,miss,'M',rate,av,'/all.txt', sep='')
				
				if (file.exists(newname)){
				cat(ab, N, miss, rate, runone(newname, ab),'\n')    } 
           
				}
			}
			}
			
		for (miss in 'MAR'){
			for (rate in c(1.5, 4, 9) ){
			for (av in c('av','noav')){
				newname<-paste('ab',ab,'N',N,miss,'M',rate,av,'/all.txt', sep='')
				
				if (file.exists(newname)){
				cat(ab, N, miss, rate, runone(newname, ab),'\n')    } 
				}          
				}
			}	
		}
	
	}
	
for (ab in c(0, .14, .39, .59)){
	for (N in c(50,100,200)){
		for (miss in c('MCAR','MAR','MNAR1', 'MNAR2')){
			for (rate in c(.1,.2,.4) ){
				for (av in c('av','noav')){
              folder.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'/all.txt',sep='')
              if (file.exists(folder.name)) file.remove(folder.name)
           folder.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'/data.txt',sep='')
              if (file.exists(folder.name)) file.remove(folder.name)
folder.name<-paste('ab',ab,'N',N,miss,'M',rate,av,'/mi.log',sep='')
              if (file.exists(folder.name)) file.remove(folder.name)
				}
			}
			}
			
	
			}	
		}
	
	}

Page Tools