===== Normal data (MCAR) =====
<code>
          [,1]       [,2]   [,3]
[1,] 0.0010376 0.01055693 0.1796
[2,] 0.0010274 0.02185791 0.1301
[3,] 0.0000272 0.01532224 0.1451
</code>

^    |     bias | MSE| power ^
| complete |0.002082 |0.01030371| 0.174 |
| MLE |0.007537 |0.02099279 |0.116 |
| Auxiliary| 0.002144| 0.01503080| 0.138|

<code>## simulation with auxiliary variables

## Generate data
gen.data<-function(){
  N<-100
  x1<-rnorm(N)
  x2<-.1*x1 + rnorm(N)
  z1<-.9*x2 + rnorm(N)
  z2<-.9*x1 + rnorm(N)
  write.table(cbind(x1,x2), 'comp.txt', na ='.', row.names=F, col.names=F)

  ## missing data in x2
  x2[sample(1:N, N*.3)]<-NA

  write.table(cbind(x1,x2,z1,z2), 'data.txt', na ='.', row.names=F, col.names=F)
}

read.mplus<-function(out){
  out<- suppressWarnings(readLines(out))
  out<- out[substr(out,1,5)=='    X']
  out<- strsplit(out, " +")
  b1<-out[[1]][3]
  b1se<-out[[1]][4]  
  c(b1, b1se)
}

## conducting the simulation
res<-NULL
for (i in 1:1000){
## generate data
gen.data()
## run the simulation
system("mplus comp.inp")
## get the results
out<-'comp.out'
temp1<-read.mplus(out)
system("mplus mle.inp")
## get the results
out<-'mle.out'
temp2<-read.mplus(out)
system("mplus regression.inp")
## get the results
out<-'regression.out'
temp3<-read.mplus(out)
temp<-as.numeric(c(temp1,temp2,temp3))
res<-rbind(res, temp)
}

## calculate bias, MSE, and power
bias1<-mean(res[,1])-.1
MSE1<-mean( (res[,1] - .1)^2 )
power1<-mean( (res[,1]-1.96*res[,2] >0)| (res[,1]+1.96*res[,2] <0) )

bias2<-mean(res[,3])-.1
MSE2<-mean( (res[,3] - .1)^2 )
power2<-mean( (res[,3]-1.96*res[,4] >0)| (res[,3]+1.96*res[,4] <0) )

bias3<-mean(res[,5])-.1
MSE3<-mean( (res[,5] - .1)^2 )
power3<-mean( (res[,5]-1.96*res[,6] >0)| (res[,5]+1.96*res[,6] <0) )

rbind(c(bias1, MSE1, power1), c(bias2, MSE2, power2), c(bias3, MSE3, power3))</code>

===== Normal data (MAR) =====


<code>
           [,1]       [,2]   [,3]
[1,] -0.0003715 0.01026548 0.1761
[2,]  0.0001807 0.02151804 0.1252
[3,]  0.0011686 0.03087685 0.1051
</code>
<code>         [,1]       [,2]  [,3]
[1,] 0.002137 0.01066757 0.189
[2,] 0.005039 0.02162380 0.124
[3,] 0.007761 0.03123450 0.100</code>
<code>
## simulation with auxiliary variables

## Generate data
gen.data<-function(){
  N<-100
  x1<-rnorm(N)
  x2<-.1*x1 + rnorm(N)
  z1<-.9*x2 + rnorm(N)
  z2<-.9*x1 + rnorm(N)
  write.table(cbind(x1,x2), 'comp.txt', na ='.', row.names=F, col.names=F)

  ## missing data in x2
  ## x2[sample(1:N, N*.3)]<-NA
  ## Missingness depend on x1
  x2[x1<quantile(x1, .3)]<-NA
  write.table(cbind(x1,x2,z1,z2), 'data.txt', na ='.', row.names=F, col.names=F)
}

read.mplus<-function(out){
  out<- suppressWarnings(readLines(out))
  out<- out[substr(out,1,5)=='    X']
  out<- strsplit(out, " +")
  b1<-out[[1]][3]
  b1se<-out[[1]][4]  
  c(b1, b1se)
}

## conducting the simulation
res<-NULL
for (i in 1:1000){
## generate data
gen.data()
## run the simulation
system("mplus comp.inp")
## get the results
out<-'comp.out'
temp1<-read.mplus(out)
system("mplus mle.inp")
## get the results
out<-'mle.out'
temp2<-read.mplus(out)
system("mplus regression.inp")
## get the results
out<-'regression.out'
temp3<-read.mplus(out)
temp<-as.numeric(c(temp1,temp2,temp3))
res<-rbind(res, temp)
}

## calculate bias, MSE, and power
bias1<-mean(res[,1])-.1
MSE1<-mean( (res[,1] - .1)^2 )
power1<-mean( (res[,1]-1.96*res[,2] >0)| (res[,1]+1.96*res[,2] <0) )

bias2<-mean(res[,3])-.1
MSE2<-mean( (res[,3] - .1)^2 )
power2<-mean( (res[,3]-1.96*res[,4] >0)| (res[,3]+1.96*res[,4] <0) )

bias3<-mean(res[,5])-.1
MSE3<-mean( (res[,5] - .1)^2 )
power3<-mean( (res[,5]-1.96*res[,6] >0)| (res[,5]+1.96*res[,6] <0) )

rbind(c(bias1, MSE1, power1), c(bias2, MSE2, power2), c(bias3, MSE3, power3))
</code>
===== MNAR data =====
<code>
           [,1]       [,2]   [,3]
[1,]  0.0004883 0.01068175 0.1809
[2,] -0.0018588 0.02183569 0.1250
[3,] -0.0482991 0.01002827 0.1037
</code>

<code>
          [,1]       [,2]  [,3]
[1,]  0.001203 0.01056922 0.192
[2,] -0.003167 0.02165664 0.126
[3,] -0.051676 0.01061156 0.110
</code>
<code>
## simulation with auxiliary variables

## Generate data
gen.data<-function(){
  N<-100
  x1<-rnorm(N)
  x2<-.1*x1 + rnorm(N)
  z1<-.9*x2 + rnorm(N)
  z2<-.9*x1 + rnorm(N)
  write.table(cbind(x1,x2), 'comp.txt', na ='.', row.names=F, col.names=F)

  ## missing data in x2
  ## x2[sample(1:N, N*.3)]<-NA
  ## Missingness depend on x1
  ## x2[x1<quantile(x1, .3)]<-NA
  ## Missingness depend on x1
  x2[x2<quantile(x2, .3)]<-NA
  write.table(cbind(x1,x2,z1,z2), 'data.txt', na ='.', row.names=F, col.names=F)
}

read.mplus<-function(out){
  out<- suppressWarnings(readLines(out))
  out<- out[substr(out,1,5)=='    X']
  out<- strsplit(out, " +")
  b1<-out[[1]][3]
  b1se<-out[[1]][4]  
  c(b1, b1se)
}

## conducting the simulation
res<-NULL
for (i in 1:1000){
## generate data
gen.data()
## run the simulation
system("mplus comp.inp")
## get the results
out<-'comp.out'
temp1<-read.mplus(out)
system("mplus mle.inp")
## get the results
out<-'mle.out'
temp2<-read.mplus(out)
system("mplus regression.inp")
## get the results
out<-'regression.out'
temp3<-read.mplus(out)
temp<-as.numeric(c(temp1,temp2,temp3))
res<-rbind(res, temp)
}

## calculate bias, MSE, and power
bias1<-mean(res[,1])-.1
MSE1<-mean( (res[,1] - .1)^2 )
power1<-mean( (res[,1]-1.96*res[,2] >0)| (res[,1]+1.96*res[,2] <0) )

bias2<-mean(res[,3])-.1
MSE2<-mean( (res[,3] - .1)^2 )
power2<-mean( (res[,3]-1.96*res[,4] >0)| (res[,3]+1.96*res[,4] <0) )

bias3<-mean(res[,5])-.1
MSE3<-mean( (res[,5] - .1)^2 )
power3<-mean( (res[,5]-1.96*res[,6] >0)| (res[,5]+1.96*res[,6] <0) )

rbind(c(bias1, MSE1, power1), c(bias2, MSE2, power2), c(bias3, MSE3, power3))
</code>