Estimating HLM Models Using Stata: Part 1

Jeremy Albright

Posted on

The data used in this tutorial can be downloaded from here. What follows replicates the results from Raudenbush and Bryk’s (2002, herafter R&B) cannonical text on hierarchical linear models (see especially chapter 4).

The Empty Model

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

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

Each school’s intercept, \(\beta_{0j}\), is then set equal to a grand mean, \(\gamma_{00}\), and a random error \(u_{0j}\).

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

Substituting (2) into (1) produces:

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

We can fit this “empty” model in Stata as follows:

mixed mathach || schid: , reml
Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log restricted-likelihood = -23558.397  
Iteration 1:   log restricted-likelihood = -23558.397  

Computing standard errors:

Mixed-effects REML regression                   Number of obs     =      7,185
Group variable: schid                           Number of groups  =        160

                                                Obs per group:
                                                              min =         14
                                                              avg =       44.9
                                                              max =         67

                                                Wald chi2(0)      =          .
Log restricted-likelihood = -23558.397          Prob > chi2       =          .

     mathach |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
       _cons |   12.63697   .2443943    51.71   0.000     12.15797    13.11598

  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
schid: Identity              |
                  var(_cons) |   8.614081   1.078813      6.739162    11.01062
               var(Residual) |   39.14832   .6606445      37.87466    40.46481
LR test vs. linear model: chibar2(01) = 986.12        Prob >= chibar2 = 0.0000

The command to fit multilevel models is mixed. The dependent variable is listed first followed by fixed effects (none in this model), followed by a double pipe (||), followed by the random effects specification. The schid: syntax specifies that schid is the grouping variable. As with all Stata commands, any modeling options follow a comma (,) after specifying the model variables. Here the reml option specifies that the model will be fit via restricted maximum likelihood rather than the default of maximum likelihood. The choice of reml is used here to be consistent with the defaults in other software.

Typically, the results from the empty model are used to determine the amount of variance in the outcome that occurs at level-2 (school) versus level-1 (student), which is known as the intraclass correlation coefficient (ICC). The output tells us that the variance component for schid, our school grouping variable, is 8.614. The level-1 variance component is 39.148.

To get the ICC manually, we would calculate:

\[ ICC = \frac{8.614}{8.614 + 39.148} = .1804 \]

The postestimation function estat icc will get this number for us.

estat icc
Intraclass correlation

                       Level |        ICC   Std. Err.     [95% Conf. Interval]
                       schid |   .1803528   .0187219      .1465168    .2199886

This tells us that 18.04% of total variability exists at the school level.

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

Still have questions? Contact us!


Stephen W. Raudenbush and Anthony S. Bryk (2002). Hierarchical Linear Models: Applications and Data Analysis Methods. Thousand Oaks: Sage.