Create fake data for hierarchical models.

fakeData.R

fakeData <- function() {
  n <- 2000
  n.groups <- 100

  x1 <<- rnorm(n); x2 <<- rnorm(n); x3 <<- rnorm(n)
  X <- cbind(1, x1, x2, x3)
  beta0 <<- runif(ncol(X), -1, 1)

  group <<- rep(1:n.groups, length.out = n)
  z1 <<- rnorm(n.groups)[group]; z2 <<- rnorm(n.groups)[group]
  alpha0 <<- rnorm(n.groups)
  gamma1 <<- rnorm(n.groups, 2, 1)
  gamma2 <<- rnorm(n.groups, -2, 1)

  mu <-  alpha0[group] + gamma1[group]*z1 + 
    gamma2[group]*z2 + X%*%beta0
  y <<- rnorm(n, mu)
}

library(arm)
library(blme)

fakeData() # simulate fake data
m1 <- lmer(y ~ x1 + x2 + x3 + z1 + z2 + (1 + z1 + z2| group)) # estimate the model
display(m) # display the results

# We can also look at bits and pieces of the estimates.
fixef(m1) # the estimated fixed effects
beta0 # compare to the true values

# Important: Each vector of random effects is stored
#   in a list. Each element of the list is devoted
#   to one grouping variable. In this case, we have
#   a single grouping variable called "group." Each 
#   element of the list is a matrix, and each column
#   of the matrix contains the random effects for
#   each variable.

# the estimated random effects by group
ranef(m1)$group[,"(Intercept)"]

alpha0 # the true random effects

plot(alpha0, ranef(m1)$group[,"(Intercept)"])