User Tools

Site Tools


lab:zhang:r_codes_for_running_mplus_for_simulation
  • 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)
    }
  }
}

Page Tools