# 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.

\(

\begin{equation} \tag{1}

Y_{ij} = \beta_{0j} + e_{ij}

\end{equation}

\)

Each school’s intercept, *β _{0j}*, is then set equal to a grand mean,

*γ*, and a random error

_{00}*u*.

_{0j}\(

\begin{equation} \tag{2}

\beta_{0j} = \gamma_{00} + u_{0j}

\end{equation}

\)

Substituting (2) into (1) produces

\(

\begin{equation} \tag{3}

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

\end{equation}

\)

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.

Still have questions? Contact us!