Tuesday, March 26, 2013
A simple simulation of Hierarchical Linear Modelling (HLM) using two different R packages - intercept only RE
# HLM is a common tool used to analysize hierarchically structured data.
# Let's see it in action using a couple of the tools programmed up by HLM researchers.
# First let's generate a data set.
# Imagine we have a two level model, students nested within classrooms.
# Let's say we have 20 classrooms
nclass = 20
# And thirty students per classroom
nstud = 30
# Let's imagine that we have a classroom effect that varies randomly and is uncorrelated with the student level effect.
class.effect = rnorm(nclass)*2
# Imagine also we have a unique outcome per student.
# We have nclass*nstud number of students in total.
student.effect = rnorm(nstud*nclass)*3
# Now we have all the data we need in order to generate our analysis.
data.set = data.frame(class.id = rep(1:nclass, each=nstud),
class.effect = rep(class.effect, each=nstud),
student.id = 1:(nstud*nclass),
student.effect = student.effect)
data.set$outcomes = data.set$class.effect + data.set$student.effect
# Looking good. Now let's load our HLM package.
# For quick reference I used:
lme(outcomes ~ 1, random = ~ 1 | class.id, data=data.set)
# We can see that lme is pretty good at estimating a standard deviation of the level 2 effect at around 2 and a level 3 effect around 3.
# Alternatively we can attempt to use the software lme4.
# For quick reference I went to:
lmer(outcomes ~ 1 + (1|class.id), data=data.set)
# Reassuringly we can see that the estimates are indentical though lme4 seems to provide more better organized feedback.