Mon - Fri : 08:00 AM - 5:00 PM
734.544-8038

# Estimating HLM Models Using R: Part 1

### Note: For a fuller treatment, download our series of lectures Hierarchical Linear Models.

#### The Empty Model

As a first step, R&B begin with an empty model containing no covariates.

$$\tag{1} Y_{ij} = \beta_{0j} + e_{ij}$$

Each school’s intercept, β0j, is then set equal to a grand mean, γ00, and a random error u0j.

$$\tag{2} \beta_{0j} = \gamma_{00} + u_{0j}$$

Substituting (2) into (1) produces

$$\tag{3} Y_{ij} = \gamma_{00} + u_{0j} + e_{ij}$$

There are actually two R packages that will estimate HLM models: nlme and lme4. As of this writing (June 2010), there is only a lme4 binary available for Windows. In addition, the comprehensive text “Mixed Effects Models in S and S-Plus” from Jose Pinheiro and Douglas Bates (2000, Springer Press) covers nlme, so it will be the focus of this tutorial. As a first step, the package must be loaded.

library(nlme)

The following estimates the empty model (assuming the data have already been read into R and stored as an object named hsb).

> model1<-lme(mathach ~ 1, random = ~ 1 | id, data=hsb)

> summary(model1)

The lme function estimates linear mixed effects models. The first argument specifies the structure of the fixed effects portion of the model using the same syntax as one uses to specify traditional linear models in R. The dependent variable is listed first, followed by a tilde, followed by the independent variables. Here the 1 tells R to estimate an intercept only. The random portion of the model is specified with similar syntax using the random = argument, but the equation is one-sided. Again, the 1 tells R to only include an intercept in this portion of the model. A vertical bar follows the listing of random effects and is in turn followed by the grouping variable. The final optional argument tells R that the variable names relate to columns in the hsb data frame. The results are stored in an object named model1 of class lme. The summary method summarizes the results. The output is the following: The output is the following:

> summary(model1)
Linear mixed-effects model fit by REML
Data: hsb
AIC      BIC    logLik
47122.79 47143.43 -23558.40

Random effects:
Formula: ~1 | id
(Intercept) Residual
StdDev:    2.934966 6.256862

Fixed effects: mathach ~ 1
Value Std.Error   DF  t-value p-value
(Intercept) 12.63697 0.2443936 7025 51.70747       0

Standardized Within-Group Residuals:
Min          Q1         Med          Q3         Max
-3.06312473 -0.75387398  0.02670132  0.76062171  2.74262579

Number of Observations: 7185
Number of Groups: 160

These results correspond to Table 4.2 in R&B. Note that the random effects are printed as standard deviations. The VarCorr function returns these estimates on both the variance and standard deviation scales.

> VarCorr(model1)
id = pdLogChol(1)
Variance  StdDev
(Intercept)  8.614025 2.934966
Residual    39.148322 6.256862

The next step is to estimate a means-as-outcomes model.