Jeremy Albright

Posted on
ANOVA HLM Stata

# Estimating HLM Models Using Stata: Part 1

The data used in this tutorial can be downloaded from (here)[https://github.com/jeralbri/tutorial-data/tree/master/data]. 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:

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.