#### 1 (a) simulate a data set

## function to simulate data and responses for a random intercept logit model 
## or log-linear poisson model
simRI <- function(  
  m = 50, # (approx.) no. of subjects  
  ni = 4, #(avg.) no. of obs. per subject 
  balanced = TRUE, # balanced design? 
  sdInt = 1, # sd(Random Intercept)
  beta =c(1), # fixed effects (multiple effects possible)  
  family=c("binomial", "poisson") # simulate binary("binomial") or count ("poisson") data?
){  
  
  family <- match.arg(family) # defaults to "binomial"
  p <- length(beta)  	
  
  #design matrices: x_i ~ U[-1,1], simulate p variables
  X <- matrix(runif(m*ni*p, -1, 1), nrow=m*ni, ncol=p)
  id <- if(balanced){
    factor(rep(1:m, each=ni)) # balanced 
  }else{
    factor(sample(1:m, m*ni, rep=TRUE)) # unbalanced
  }
  
  # U = incidence matrix for id-factor
  U <- model.matrix(~ 0 + id)
  # draw random intercepts ~N(0, sdInt^2)
  b <- rnorm(nlevels(id), sd = sdInt)
  
  # generate y
  if(family == "binomial"){	
    mu <- plogis(X %*% beta + U %*% b) # plogis(x) is inverse logit function
    y <- rbinom(m*ni, 1, mu)
  }
  if(family == "poisson"){
    mu <- exp(X %*% beta + U %*% b)
    y <- rpois(m*ni, mu)
  }
  
  return(data.frame(y, X=X, id))
}

data <- simRI(m = 50,  ni = 4, balanced = TRUE, 
              sdInt = 1, beta = c(2), family="binomial")


#### 1 (b) estimate the models and look at coefficients

# fit the GLMM via PQL 
require(MASS)
require(nlme)
fitPQL <- try(glmmPQL(y ~ X, random = ~ 1|id, family = "binomial", data = data))
fixef(fitPQL)
VarCorr(fitPQL)

# fit the GLMM via Laplace
require(lme4)
fitLaplace <- try(glmer(y ~ X + (1|id), family = "binomial", data = data))
fixef(fitLaplace)
VarCorr(fitLaplace)

# fit the GLMM via Adaptive Gauss Quadratur
require(lme4)
fitAGQ <- try(glmer(y ~ X + (1|id), family = "binomial", data = data, nAGQ = 11))
fixef(fitAGQ)
VarCorr(fitAGQ)

# true values for fixed effects are 0 and 2
# true variance of random effect is 1


#### 1 (c) example to change m, ni

data <- simRI(m = 100,  ni = 5, balanced = TRUE, 
              sdInt = 1, beta =c(2), family="binomial")

# fit the GLMM via PQL 
require(MASS)
require(nlme)
fitPQL <- try(glmmPQL(y ~ X, random = ~ 1|id, family = "binomial", data = data, verbose = F))
fixef(fitPQL)
VarCorr(fitPQL)

# fit the GLMM via Laplace
require(lme4)
fitLaplace <- try(glmer(y ~ X + (1|id), family = "binomial", data = data))
fixef(fitLaplace)
VarCorr(fitLaplace)

# fit the GLMM via Adaptive Gauss Quadratur
require(lme4)
fitAGQ <- try(glmer(y ~ X + (1|id), family = "binomial", data = data, nAGQ = 11))
fixef(fitAGQ)
VarCorr(fitAGQ)

# true values for fixed effects are 0 and 2
# true variance of random effect is 1