All results in Excel table

1. Stars plot for result comparison
2. Simulation codes on 2010-10-25

1. SAS macros for listwise deletion
2. R codes for running the listwise deletion simulation

1. SAS macros for pairwise deletion
2. R codes for running the pairwise deletion simulation

1. Run BMEME in R for simulation
2. BMEM simulation setup using R
3. Process BMEM results for BMEM paper
4. Results from BMEM as of April 06, 2010 (no missing in X)
5. Plot running time for BMEM

• MI results with 100 imputations (no missing in X)

• SAS macros for simulation
• R codes for running MI simulation
• Generate sas and R codes for simulation
• Submit all jobs
• SAS codes with auxiliary variables
• SAS codes without auxiliary variables
• Submit one job
• Remove all files
• R codes to run simulation
• Simulation data for MNAR condition
• Data for imputation numbers of 25 and 50

• R codes for processing results of MI using SAS
• R codes for plot of relative deviance according to number of imputations

• R codes for running Mplus for simulation
• R codes for processing results from Mplus
• Simulation results as of April 07 2010 (no missing in X)
• An R function to run mplus and obtain BC intervals

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

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)