Normal data (MCAR)
[,1] [,2] [,3]
[1,] 0.0010376 0.01055693 0.1796
[2,] 0.0010274 0.02185791 0.1301
[3,] 0.0000272 0.01532224 0.1451
| | bias | MSE | power |
| complete | 0.002082 | 0.01030371 | 0.174 |
| MLE | 0.007537 | 0.02099279 | 0.116 |
| Auxiliary | 0.002144 | 0.01503080 | 0.138 |
## 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))
Normal data (MAR)
[,1] [,2] [,3]
[1,] -0.0003715 0.01026548 0.1761
[2,] 0.0001807 0.02151804 0.1252
[3,] 0.0011686 0.03087685 0.1051
[,1] [,2] [,3]
[1,] 0.002137 0.01066757 0.189
[2,] 0.005039 0.02162380 0.124
[3,] 0.007761 0.03123450 0.100
## 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))
MNAR data
[,1] [,2] [,3]
[1,] 0.0004883 0.01068175 0.1809
[2,] -0.0018588 0.02183569 0.1250
[3,] -0.0482991 0.01002827 0.1037
[,1] [,2] [,3]
[1,] 0.001203 0.01056922 0.192
[2,] -0.003167 0.02165664 0.126
[3,] -0.051676 0.01061156 0.110
## 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))