
B<-100

N<-1000
R<-0.3

runsim<-function(B, N, R){
  
  ## (1) Simulation based on different data sets
  b.mcar.all<-NULL
  b.mar.all<-NULL
  
  for (i in 1:B){

     x<-rnorm(N)
     y<-1+2*x +rnorm(N)

     ## estimation based on MLE
     ## MAR data
     y1<-y
     c.value<-x[rank(x)==N*R]
     y1[x < c.value]<-NA

     b.mar<-coef(lm(y1~x))
     b.mar.all<-rbind(b.mar.all, b.mar)

     ##MCAR data
     x<-rnorm(N)
     y<-1+2*x +rnorm(N)

     ## estimation based on MLE
     #MCAR data
     y2<-y
     y2[sample(1:N, N*R)]<-NA

     b.mcar<-coef(lm(y2~x))

     b.mcar.all<-rbind(b.mcar.all, b.mcar)
  }

  m.mcar<-apply(b.mcar.all, 2, mean)
  sd.mcar<-apply(b.mcar.all, 2, sd)

  m.mar<-apply(b.mar.all, 2, mean)
  sd.mar<-apply(b.mar.all, 2, sd)

  res1<-c(m.mcar-c(1,2), m.mar-c(1,2), sd.mar-sd.mcar)
  ## rbind(m.mcar, sd.mcar, m.mar, sd.mar)


  ## (2) using same datasets

  b.mcar.all<-NULL
  b.mar.all<-NULL


  for (i in 1:B){

     x<-rnorm(N)
     y<-1+2*x +rnorm(N)

     ## estimation based on MLE
     #MAR data
     y1<-y
     c.value<-x[rank(x)==N*R]
     y1[x < c.value]<-NA

     b.mar<-coef(lm(y1~x))
     b.mar.all<-rbind(b.mar.all, b.mar)

     ##MCAR data
     ## estimation based on MLE
     #MAR data
     y2<-y
     y2[sample(1:N, N*R)]<-NA

     b.mcar<-coef(lm(y2~x))

     b.mcar.all<-rbind(b.mcar.all, b.mcar)
  }

  m.mcar<-apply(b.mcar.all, 2, mean)
  sd.mcar<-apply(b.mcar.all, 2, sd)

  m.mar<-apply(b.mar.all, 2, mean)
  sd.mar<-apply(b.mar.all, 2, sd)

  res2<-c(m.mcar-c(1,2), m.mar-c(1,2), sd.mar-sd.mcar)
  ##rbind(m.mcar, sd.mcar, m.mar, sd.mar)
  c(res1,res2)
}

allres<-NULL
for (k in 1:1000){
	temp<-runsim(100,100,.3)
	allres<-cbind(allres, temp)
	}