lab:projects:15auxiliary_variables [Zhiyong Johnny Zhang]

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))

Table of Contents

Normal data (MCAR)

Normal data (MAR)

MNAR data