### Site Tools

lab:projects:15auxiliary_variables

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

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'
system("mplus mle.inp")
## get the results
out<-'mle.out'
system("mplus regression.inp")
## get the results
out<-'regression.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)
}

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'
system("mplus mle.inp")
## get the results
out<-'mle.out'
system("mplus regression.inp")
## get the results
out<-'regression.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)
}

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'
system("mplus mle.inp")
## get the results
out<-'mle.out'
system("mplus regression.inp")
## get the results
out<-'regression.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))