###################################################
### chunk number 1: prelim
###################################################
rm(list=ls())
try(detach(package:nlme))
try(detach(package:lme4))
setwd('W:/lehre/longitudinal_sommer08/uebung')
options(width = 80, prompt = "R>")


###################################################
### chunk number 2: 
###################################################
set.seed(123)

######################################################################
# Simulate complete data in wide format from a LME with
# just a random intercept
#
# Params:
#  N - number of individuals
#  n - number of observations per individual
#  dbetween - variance of the random intercept b_i ~ N(0,dbetween)
#  dwithin  - variance of error term
#
# Returns:
#  list with two vectors of length N containing the observation
#  and the corresponding random effect of each individual
######################################################################

sim.basiclme <- function(N=30,n=2, beta=0, dbetween =2.5, dwithin= 0.5) {
  y <- matrix(0, nrow = N, ncol = n)
  b <- rnorm(N, 0, sqrt(dbetween))
  #First simulate as if all people had complete observations (no missing)
  for (i in 1:N) {
    y[i,] <- beta + b[i] + rnorm(n, 0, sqrt(dwithin))
  }
  row.names(y) <- c(1:nrow(y))
  colnames(y) <- paste("t.",0:(n-1),sep="")
  return(list(y=y,b=b))
}

data <- sim.basiclme()
y <- data$y
b <- data$b

sim1 <- y



###################################################
### chunk number 3: 
###################################################

makeMissing <- function(data, beta)
{
pi <- plogis( cbind(1, data) %*% beta )
R  <- rbinom(nrow(data),size=1, prob=pi)
data[R==0,2] <- NA
data <- as.data.frame(cbind(data,R))
}

beta <- matrix(c(0,0,1),3,1)

sim2 <- makeMissing(sim1,beta)


###################################################
### chunk number 4: 
###################################################

makeMissing2 <- function(data, b, alpha)
{
pi <- plogis( cbind(1,b) %*% alpha )
R  <- rbinom(nrow(data),size=1, prob=pi)
data[R==0,2] <- NA
data <- as.data.frame(cbind(data,R))
}
alpha <- matrix(c(-1,1),2,1)
sim3 <- makeMissing2(sim1, b, alpha)


###################################################
### chunk number 5: 
###################################################
#Convert data to long format
sim2long <- function(y) {
  require(doBy)
  y.long <- reshape(y,varying=list(c("t.0","t.1")),v.names="y",direction="long")
  return( groupedData(y~1|id, data=orderBy(~id + t, y.long)))
}

library(nlme)
(m1 <- lme(y ~ 1, random= ~1|id, data=sim2long(as.data.frame(sim1))))
(m2 <- lme(y ~ 1, random= ~1|id, data=sim2long(sim2),na.action=na.exclude))
(m3 <- lme(y ~ 1, random= ~1|id, data=sim2long(sim3),na.action=na.exclude))


###################################################
### chunk number 6: 
###################################################
nsim <- 500

res1<- res2 <- res3 <- matrix(NA,nrow=nsim,ncol=3)

#Get beta_0, var(b) and var(error) from lme-Object
getParams<-function(m)
{
        beta <- fixef(m)
        varb <- getVarCov(m)
        vareps <- m$sigma^2
        return(c(beta, varb, vareps))
}

for(i in 1:nsim)
{
  data <- sim.basiclme()
  y <- data$y
  b <- data$b

  sim1 <- y
  sim2 <- makeMissing(sim1, beta = c(0,0,1))
  sim3 <- makeMissing2(sim1, b, alpha)
  m1 <- lme(y ~ 1, random= ~1|id, data=sim2long(as.data.frame(sim1)))
  res1[i, ] <- getParams(m1)

        m2 <- lme(y ~ 1, random= ~1|id, data=sim2long(sim2), na.action=na.exclude)
  res2[i, ] <- getParams(m2)

        m3 <- lme(y ~ 1, random= ~1|id, data=sim2long(sim3), na.action=na.exclude)
     res3[i, ] <- getParams(m3)
     if(i %% 50 == 0) cat(i,' ')
}

beta <- data.frame(complete=res1[,1],NMAR1=res2[,1],NMAR2=res3[,1])
(bias.beta <- colSums(beta)/nsim)

varb <- data.frame(complete=res1[,2],NMAR1=res2[,2],NMAR2=res3[,2])
(bias.varb <- colSums((varb-2.5))/nsim)


vareps <- data.frame(complete=res1[,3],NMAR1=res2[,3],NMAR2=res3[,3])
(bias.vareps <- colSums((vareps-.5))/nsim)

par(mfrow=c(3,1))
boxplot(beta,main='Intercept',medlwd=.5)
abline(h=0,col='lightgrey')
boxplot(varb,main='Var(b)',medlwd=.5)
abline(h=2.5,col='lightgrey')
boxplot(vareps,main='Var(eps)',medlwd=.5)
abline(h=.5,col='lightgrey')