===== Simulation results =====
{{:lab:zhang:results_from_bmem.xls|All results in Excel table}}
==== Data simulation ====
  - [[Stars plot for result comparison]]
  - [[Simulation codes on 2010-10-25]]

==== Listwise deletion ====
  - [[SAS macros for listwise deletion]]
  - [[R codes for running the listwise deletion simulation]]
==== Pairwise deletion ====
  - [[SAS macros for pairwise deletion]]
  - [[R codes for running the pairwise deletion simulation]]
==== BMEM ====
  - [[Run BMEME in R for simulation]]
  - [[BMEM simulation setup using R]]
  - [[Process BMEM results for BMEM paper]]
  - {{:lab:zhang:bmem.results.2010april06.zip|Results from BMEM as of April 06, 2010}} (no missing in X)
  - [[Plot running time for BMEM]]
==== Multiple imputation method ====
=== Results ===
  * {{:lab:zhang:mi.results.2010.04.22.zip|MI results with 100 imputations (no missing in X)}}

=== Files for simulation ===
  * [[SAS macros for simulation]]
  * [[R codes for running MI simulation]]
  * {{:lab:zhang:allsim.r|Generate sas and R codes for simulation}}
  * {{:lab:zhang:allsub.sub|Submit all jobs}}
  * {{:lab:zhang:mi.aux.sas|SAS codes with auxiliary variables}}
  * {{:lab:zhang:mi.sas|SAS codes without auxiliary variables}}
  * {{:lab:zhang:qsub.sub|Submit one job}}
  * {{:lab:zhang:rm.sub|Remove all files}}
  * {{:lab:zhang:runsim.r|R codes to run simulation}} 
  * {{:lab:zhang:mnar.zip|Simulation data for MNAR condition}}
  * {{:lab:zhang:newdata-for-xin.rar|Data for imputation numbers of 25 and 50}}
=== Some codes ===
  * [[R codes for processing results of MI using SAS]]
  * [[R codes for plot of relative deviance according to number of imputations]]

==== FIML in Mplus ====
  * [[R codes for running Mplus for simulation]]
  * [[R codes for processing results from Mplus]]
  * {{:lab:zhang:fiml2010april06.rar|Simulation results as of April 07 2010}} (no missing in X)
  * [[An R function to run mplus and obtain BC intervals]]



===== A simple example looking at the influence of MNAR auxiliary variables on regression coefficients =====
<code>
## a simple simulation

proc.res<-function(out.file){
  output <- suppressWarnings(readLines(out.file))
  lnumber <- grep("    X", output)
  est<- output[lnumber[2]]
  sub('X', '', est)
}

res1<-NULL
res2<-NULL

for (i in 1:1000){

N<-500
x<-rnorm(N)
y<-0*x + rnorm(N)

z<-(x+y)+rnorm(N)

cor(cbind(x,y,z))

x[sample(1:N, N*.3)]<-NA
y[sample(1:N, N*.3)]<-NA
z[z>quantile(z, .7)]<-NA

write.table(cbind(x,y,z), 'testreg.txt', row.names=F, col.names=F, na='.')

system("c:/programs/mplus/mplus reg.inp",show.output.on.console=F)
temp1<-proc.res('reg.out')
res1<-cbind(res1, temp1)

system("c:/programs/mplus/mplus reg-aux.inp",show.output.on.console=F)
temp2<-proc.res('reg-aux.out')
res2<-cbind(res2, temp2)
cat(i, '\n')

}

write.table(res1, 'res1.txt', row.names=F, col.names=F, quote=F)
write.table(res2, 'res2.txt', row.names=F, col.names=F, quote=F)

res1<-scan('res1.txt')
res1<-array(res1, dim=c(4,1000))

res2<-scan('res2.txt')
res2<-array(res2, dim=c(4,1000))

apply(res1,1,mean)
apply(res2,1,mean)
</code>
==== Using EM in R ====
<code>
library(norm)
library(MASS)
library(mvtnorm)



N<-500
r<-.9
set.seed(0)
allres<-NULL
for (i in 1:1000){

  sigma <- matrix(c(1,0.3,0.3,1),2,2)
  yx	<-	mvrnorm(N, rep(0,2),sigma)
  yx[sample(1:N, N*r), 1]<-NA
  yx[yx[, 2]>quantile(yx[,2], 1-r), 2]<-NA
  
  m<-prelim.norm(yx)
  est<-em.norm(m)
  
  res<-getparam.norm(m, est)
  
  res<-c(unlist(res), mean(yx[,1], na.rm=T), var(yx[,1], na.rm=T))
  allres<-rbind(allres, res)
}

apply(allres, 2, mean)

## function to get the likelihood function

llik<-function(x, par){
	#par[1]: mu1
	#par[2]: mu2
	#par[3]: s1
	#par[4]: s2
	#par[5]: s12
	
	mu<-par[1:2]
	sigma <- array(c(par[3], par[5], par[5], par[4]), dim=c(2,2))
	
	n<-dim(x)[1]
	
	l<-0
	for (i in 1:n){
		if ((!is.na(x[i,1])) & (!is.na(x[i,2]))){
			l<-l +  dmvnorm(x[i,], mu, sigma, log=T)
			}
		if ((!is.na(x[i,1])) & (is.na(x[i,2]))){
			l<-l + dnorm(x[i,1], par[1], sqrt(par[3]), log=T)
			}
		if ((is.na(x[i,1])) & (!is.na(x[i,2]))){
			l<-l + dnorm(x[i,2], par[2], sqrt(par[4]), log=T)
			}
	}
l
}

nlm(llik, c(0,0,1,1,0), x=yx)
	
</code>