## Program to run Mplus for the simulation

## write the model

gen.model<-function( ){
## generate mplus scripts for data analysis
## Mplus input scripts
mplus.null<-'model.inp'

## for null model
cat('TITLE: Model\n',file=mplus.null); 
cat('DATA:\n',file=mplus.null, append=T); 
cat('  FILE=data.txt;\n',file=mplus.null, append=T); 
cat('VARIABLE: \n',file=mplus.null, append=T); 
cat('    NAMES ARE X M Y U V W;\n', file=mplus.null, append=T);
cat('    USEVARIABLES ARE X M Y;\n', file=mplus.null, append=T);
cat('    AUXILIARY = (m) U V W;\n', file=mplus.null, append=T);
cat('    MISSING = ALL(99999); \n',file=mplus.null, append=T); 

cat('ANALYSIS: \n',file=mplus.null, append=T); 
cat('    COVERAGE=.01; \n',file=mplus.null, append=T); 
## cat('    bootstrap=1000; \n',file=mplus.null, append=T);

cat('MODEL: \n',file=mplus.null, append=T); 
cat('    M ON X ;\n',file=mplus.null, append=T); 
cat('    Y ON M X; \n',file=mplus.null, append=T); 
cat('MODEL INDIRECT: \n',file=mplus.null, append=T); 
cat('    Y IND X; \n',file=mplus.null, append=T); 

## cat('OUTPUT:\n',  file=mplus.null, append=T);
## cat('    CINTERVAL (BCBOOTSTRAP); \n', file=mplus.null, append=T);
}

## function to process output to get the parameter estimates and the bc intervals

proc.res<-function(){
  out.file<-"model.out"
  output <- suppressWarnings(readLines(out.file))
  lnumber <- grep("  Total indirect", output)
  est<- unlist(strsplit(output[lnumber[1]], " "))
  est<-est[!est==""]
  est<-as.numeric(est[3])
  est
}


## function to calculate the BC for mplus
bc.ci<-function(res, est){
  iB<-length(res)
  ## BC interval
  alpha<-c(.025, .975)
  z0<-qnorm(sum(res<est)/iB)
  zalpha<-qnorm(alpha)
  bcalpha<-pnorm(2*z0+zalpha)
  BC<-quantile(res,bcalpha)
 
  c(BC, est)
}


## run the simulation
## run the bootstrap to get sd, ci for indirect effect
boot.once<-function(i, B, N){
  datafile <- paste('C:\\zzy\\MNARmediation\\mnar\\data-',i,'-',r,'.txt', sep='')
  temp.data<-read.table(datafile)
     write.table(temp.data, 'data.txt', row.names=F, col.names=F)
     system("c:/programs/mplus/mplus model.inp",show.output.on.console=F)
     res<-proc.res()
  boot.res<-NULL
  for (j in 1:B){
     index<-sample(1:N, replace=T)
     dset<-temp.data[index, ]
     write.table(dset, 'data.txt', row.names=F, col.names=F)
     system("c:/programs/mplus/mplus model.inp",show.output.on.console=F)
     temp<-proc.res()
     boot.res<-c(boot.res, temp)
  }
  bc.ci(boot.res, res)
}

##mis.pat<-c('mar', 'mcar', 'mnar')
mis.pat<-'mnar'
B<-1000
N<-100
rate1<-c(.4)
rate2<-c(.1,.2,.4)
gen.model()

for (pat in mis.pat){

  if (pat == 'mar'){
    rate <- rate1
  }else{
    rate <- rate2
  }

  for (r in rate){
    resfile <- paste(pat ,'-',r,'.txt', sep='')
    for (k in 1:1000){
       temp<-boot.once(k, B, N)
       cat(c(k, temp), file=resfile, append=T)
       cat("\n", file=resfile,append=T)
    }
  }
}



res<-read.table('mnar-0.4.txt')
sum(res[,2]<.39^2 & res[,3]>.39^2)/dim(res)[1]
sum(res[,2]>0)/dim(res)[1]