* For data without auxiliary variables
## Program to run Mplus for the simulation
## write the model
gen.model<-function(datafile){
## 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=', datafile, ';\n',file=mplus.null, append=T);
cat('VARIABLE: \n',file=mplus.null, append=T);
cat(' NAMES ARE X M Y;\n', file=mplus.null, append=T);
cat(' USEVARIABLES ARE X M Y;\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<- output[lnumber[2]] #as.numeric(strsplit(output[lnumber[2]], " ")[[1]][c(16,22,29)])
est
}
## run the simulation
mis.pat<-c('mar', 'mcar', 'mnar')
rate1<-c(1.5, 4, 9)
rate2<-c(.1,.2,.4)
for (pat in mis.pat){
if (pat == 'mar'){
rate <- rate1
}else{
rate <- rate2
}
for (r in rate){
resfile <- paste(pat ,'-',r,'.txt', sep='')
for (i in 1:1000){
datafile <- paste('C:\\zzy\\MNARmediation\\newdata\\', pat ,'\\data-',i,'-',r,'.txt', sep='')
gen.model(datafile)
system("c:/programs/mplus/mplus model.inp",show.output.on.console=F)
temp<-proc.res()
cat(c(i, temp), file=resfile, append=T)
cat("\n", file=resfile,append=T)
}
}
}
* For data with auxiliary variables
## Program to run Mplus for the simulation
## write the model
gen.model<-function(datafile){
## 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=', datafile, ';\n',file=mplus.null, append=T);
cat('VARIABLE: \n',file=mplus.null, append=T);
cat(' NAMES ARE X M Y U V;\n', file=mplus.null, append=T);
cat(' USEVARIABLES ARE X M Y U V;\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(' U with X M Y; \n',file=mplus.null, append=T);
cat(' V with X M Y; \n',file=mplus.null, append=T);
cat(' U with V; \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<- output[lnumber[2]] #as.numeric(strsplit(output[lnumber[2]], " ")[[1]][c(16,22,29)])
est
}
## run the simulation
mis.pat<-c('mcar0','mcar1','mcar2','mar0','mar1','mar2','mnar0','mnar1','mnar2','mnar0-1','mnar1-1','mnar2-1')
rate1<-c(1.5, 4, 9)
rate2<-c(.1,.2,.4)
for (pat in mis.pat){
if (substr(pat,1,3)=="mar"){
rate <- rate1
}else{
rate <- rate2
}
for (r in rate){
resfile <- paste(pat ,'-',r,'.txt', sep='')
for (i in 1:1000){
datafile <- paste('C:\\zzy\\MNARmediation\\newdata\\', pat ,'\\data-',i,'-',r,'.txt', sep='')
gen.model(datafile)
system("c:/programs/mplus/mplus model.inp",show.output.on.console=F)
temp<-proc.res()
cat(c(i, temp), file=resfile, append=T)
cat("\n", file=resfile,append=T)
}
}
}