## Program to run Mplus for the simulation

## 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(est, sd(res), BC)
}


## run the simulation
## run the bootstrap to get sd, ci for indirect effect
boot.mplus<-function(B, N, datafile, scriptfile, resfile){
  temp.data<-read.table(datafile)
  write.table(temp.data, 'data.txt', row.names=F, col.names=F)
  runmplus<-paste("c:/programs/mplus/mplus ", scriptfile)
  system(runmplus,show.output.on.console=F)
  if (file.exists(resfile)){
    res<-scan(resfile)
    unlink(resfile)
  }else{
    cat('Problem with original data analysis')
    break
  }

  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(runmplus,show.output.on.console=F)
     temp<-scan(resfile)
     boot.res<-rbind(boot.res, temp)
     cat(j, '\n')
  }

  ## Calculate confidence intervals
  res<-res[c(19:21,24,25,27,31,33,36)]
  res<-c(res, res[4]*res[5])

  boot.res<-boot.res[,c(19:21,24,25,27,31,33,36)]
  boot.res<-cbind(boot.res, boot.res[,4]*boot.res[,5])
  
  ## BC in a matrix
  BC.matrix<-NULL

  for (i in 1:10){
    BC.matrix<-rbind(BC.matrix, bc.ci(boot.res[,i], res[i]))  
  }

  colnames(BC.matrix)<-c('Estimate','S.E.', '2.5%', '97.5%')
  rownames(BC.matrix)<-c('iM','iY','iX', 'a','b',"c'",'eM2','eY2','eX2','ab')

  list(res=res, boot.res=boot.res, BC=BC.matrix)
}

## Run bootstrap
system.time(
allres<-boot.mplus(1000,2802,'active.txt','MODEL.ACTIVE.AUX.inp','est.txt')
)