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

### Citations

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