Stata’s new mixed-models estimation makes it easy to specify and to fit two-way, multilevel, and hierarchical random-effects models.To fit a model of SAT scores with fixed coefficient on x1 and random coefficient on x2 at the school level, and with random intercepts at both the school and class-within-school level, you type

Estimation is either by ML or REML, and standard errors and confidence intervals (in fact, the full covariance matrix) are estimated for all variance components.. xtmixed sat x1 x2 || school: x2 || class:

xtmixedprovides four random-effects variance structures—identity, independent, exchangeable, and unstructured—and you can combine them to form even more complex block-diagonal structures.After estimation, you can obtain best linear unbiased predictions (BLUPs) of random effects and of conditional means (fitted values).

For details, see[XT] xtmixed.

After estimation withxtmixed,

estatrecovariancereports the estimated variance–covariance matrix of the random effects for each level.

estatgroupsummarizes the composition of the nested groups, providing minimum, average, and maximum group size for each level in the model.

predictafterxtmixedcan compute best linear unbiased predictions (BLUPs) for each random effect. It can also compute the linear predictor, the standard error of the linear predictor, the fitted values (linear predictor plus contributions of random effects), the residuals, and the standardized residuals.

Here are some examples from thextmixedmanual entry.

Example 1: One-level model

Consider a longitudinal dataset used by both Ruppert, Wand, and Carroll (2003) and Diggle et al. (2002), consisting of weight measurements of 48 pigs on nine successive weeks. Pigs are identified by variableid. Below is a plot of the growth curves for the first ten pigs.

. use http://www.statapress.com/data/r9/pig (Longitudinal analysis of pig weights) . twoway scatter weight week if id<=10, connect(L)

It seems clear that each pig experiences a linear trend in growth, and that overall weight measurements vary from pig to pig. Because we are not really interested in these particular 48 pigsper se, we instead treat them as a random sample from a larger population and model the between-pig variability as a random effect or, equivalently, as a random-intercept term at the pig level. We thus wish to fit the following modelweight_forij=b_0 +b_1 week_ij+u_i+e_ij(3)i= 1,...,48 andj= 1,...,9. The fixed portion of the model,b_0 +b_1*week_ij, simply states that we want one overall regression line representing the population average. The random effectu_iserves to shift this regression line up or down for each pig. Since the random effects occur at the pig-level, we fit the model by typing

. xtmixed weight week || id: Performing EM optimization: Performing gradient-based optimization: Iteration 0: log restricted-likelihood = -1016.8984 Iteration 1: log restricted-likelihood = -1016.8984 Computing standard errors: Mixed-effects REML regression Number of obs = 432 Group variable: id Number of groups = 48 Obs per group: min = 9 avg = 9.0 max = 9 Wald chi2(1) = 25271.50 Log restricted-likelihood = -1016.8984 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ weight | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- week | 6.209896 .0390633 158.97 0.000 6.133333 6.286458 _cons | 19.35561 .603139 32.09 0.000 18.17348 20.53774 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ id: Identity | sd(_cons) | 3.891253 .4143198 3.158334 4.794252 -----------------------------+------------------------------------------------ sd(Residual) | 2.096356 .0757444 1.953034 2.250195 ------------------------------------------------------------------------------ LR test vs. linear regression: chibar2(01) = 473.15 Prob >= chibar2 = 0.0000 . estimates store randint

At this point, a guided tour of the model specification and output is in order:

- By typing "
weight week" we specified the response,weight, and the fixed portion of the model in the same way we would if we were usingregressor any other estimation command. Our fixed effects are a coefficient onweekand a constant term.

- When we added "
|| id:", we specified random effects at the level identified by group variableid, i.e., the pig-level. Since we only wanted a random intercept, that is all we had to type.

- The output title informs us that our model was fitted using REML, the default. For ML estimates, use option
mle. Since this model is a simple random-intercept model, specifying optionmlewould be equivalent to usingxtreg, also with optionmle.

- The first estimation table is for the fixed effects. We estimate
b_0 = 19.36 andb_1 = 6.21.

- The second estimation table shows the estimated variance components. Since we only have one random effect at this level,
xtmixedknew that theIdentitycovariance structure is the only one that makes sense. In any case,sigma_u is estimated as 3.89 with standard error 0.414.

If you prefer variance estimates, e.g.,sigma^2_u, to standard deviation estimates, e.g.,sigma_u, specify optionvariance, either at estimation or upon replay.

- The row labeled "sd(Residual)" displays the estimated standard deviation of the overall error term, i.e.,
sigma_epsilon= 2.10.

- Finally, a likelihood-ratio test comparing the model to ordinary linear regression is provided and is highly significant for these data.

Example 2: One-level model

Extending (3) to also allow for a random slope onweekyields the modelweight_fitted usingij=b_0 +b_1 week_ij+u_oi+u_1iweek_ij+e_ij(4)xtmixed

. xtmixed weight week || id: week Performing EM optimization: Performing gradient-based optimization: Iteration 0: log restricted-likelihood = -870.51473 Iteration 1: log restricted-likelihood = -870.51473 Computing standard errors: Mixed-effects REML regression Number of obs = 432 Group variable: id Number of groups = 48 Obs per group: min = 9 avg = 9.0 max = 9 Wald chi2(1) = 4592.10 Log restricted-likelihood = -870.51473 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ weight | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- week | 6.209896 .0916386 67.77 0.000 6.030287 6.389504 _cons | 19.35561 .4021142 48.13 0.000 18.56748 20.14374 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ id: Independent | sd(week) | .6135471 .0673971 .4947035 .7609409 sd(_cons) | 2.630132 .302883 2.098719 3.296105 -----------------------------+------------------------------------------------ sd(Residual) | 1.26443 .0487971 1.172317 1.363781 ------------------------------------------------------------------------------ LR test vs. linear regression: chi2(2) = 765.92 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference . estimates store randslope

Since we didn’t specify a covariance structure for the random effects(u_0i, u_1i),xtmixedused the defaultIndependentstructure and estimatedsigma_u0 = 2.63 andsigma_u1 = 0.61. Our point estimates of the fixed effects are, for all intents and purposes, identical, but this does not hold generally. Given the 95% confidence interval forsigma_u1, it would seem that the random slope is significant, and we can uselrtestto verify this fact:

. lrtest randslope randint Likelihood-ratio test LR chibar2(01) = 292.77 (Assumption: randint nested in randslope) Prob > chibar2 = 0.0000 Note: LR tests based on REML are valid only when the fixed-effects specification is identical for both models

Example 3: Two-level model

Baltagi, Song, and Jung (2001) estimate a Cobb–Douglas production function examining the productivity of public capital in each state’s private output. Originally provided by Munnell (1990), the data were recorded over the period 1970–1986 for 48 states grouped into nine regions.

. webuse productivity (Public Capital Productivity) . describe Contains data from http://localpress.stata.com/data/r9/productivity.dta obs: 816 Public Capital Productivity vars: 11 29 Mar 2005 10:57 size: 32,640 (96.9% of memory free) (_dta has notes) ------------------------------------------------------------------------------- storage display value variable name type format label variable label ------------------------------------------------------------------------------- state byte %9.0g states 1-48 region byte %9.0g regions 1-9 year int %9.0g years 1970-1986 public float %9.0g public capital stock hwy float %9.0g log(highway component of public) water float %9.0g log(water component of public) other float %9.0g log(bldg/other component of public) private float %9.0g log(private capital stock) gsp float %9.0g log(gross state product) emp float %9.0g log(non-agriculture payrolls) unemp float %9.0g state unemployment rate ------------------------------------------------------------------------------- Sorted by:

Since the states are nested within regions, we fit a two-level mixed model with random intercepts at both the region and at the state-within-region levels.

. xtmixed gsp private emp hwy water other unemp || region: || state: Performing EM optimization: Performing gradient-based optimization: Iteration 0: log restricted-likelihood = 1404.7101 Iteration 1: log restricted-likelihood = 1404.7101 Computing standard errors: Mixed-effects REML regression Number of obs = 816 ----------------------------------------------------------- | No. of Observations per Group Group Variable | Groups Minimum Average Maximum ----------------+------------------------------------------ region | 9 51 90.7 136 state | 48 17 17.0 17 ----------------------------------------------------------- Wald chi2(6) = 18382.39 Log restricted-likelihood = 1404.7101 Prob > chi2 = 0.0000 ------------------------------------------------------------------------------ gsp | Coef. Std. Err. z P>|z| [95% Conf. Interval] -------------+---------------------------------------------------------------- private | .2660308 .0215471 12.35 0.000 .2237993 .3082624 emp | .7555059 .0264556 28.56 0.000 .7036539 .8073579 hwy | .0718857 .0233478 3.08 0.002 .0261249 .1176464 water | .0761552 .0139952 5.44 0.000 .0487251 .1035853 other | -.1005396 .0170173 -5.91 0.000 -.1338929 -.0671862 unemp | -.0058815 .0009093 -6.47 0.000 -.0076636 -.0040994 _cons | 2.126995 .1574864 13.51 0.000 1.818327 2.435663 ------------------------------------------------------------------------------ ------------------------------------------------------------------------------ Random-effects Parameters | Estimate Std. Err. [95% Conf. Interval] -----------------------------+------------------------------------------------ region: Identity | sd(_cons) | .0435471 .0186292 .0188287 .1007161 -----------------------------+------------------------------------------------ state: Identity | sd(_cons) | .0802737 .0095512 .0635762 .1013567 -----------------------------+------------------------------------------------ sd(Residual) | .0368008 .0009442 .034996 .0386987 ------------------------------------------------------------------------------ LR test vs. linear regression: chi2(2) = 1162.40 Prob > chi2 = 0.0000 Note: LR test is conservative and provided only for reference

Some items of note:

- Our model now has two random-effects equations, separated by
||. The first is a random intercept (constant-only) at theregionlevel, the second a random intercept at thestatelevel. The order in which these are specified (from left to right) is significant—xtmixedassumes thatstateis nested withinregion.

- The variance component estimates are now organized and labeled according to level. After adjusting for the nested-level error structure, we find that the highway and water components of public capital had significant positive effects upon private output, while the other public buildings component had a negative effect.

## References

- Baltagi, B. H., S. H. Song, and B. C. Jung. 2001.
Journal of Econometrics101: 357–381.

- Munnell, A. 1990.
- Why has productivity growth declined? Productivity and public investment.
New England Economic Review. Jan. / Feb.:3–22.

- Diggle, P. J., P. Heagerty, K.-Y. Liang, and S. L. Zeger. 2002.
Analysis of Longitudinal Data. 2nd edition. Oxford: Oxford University Press.

- Ruppert, D., M. P. Wand, and R. J. Carroll. 2003.
Semiparametric Regression. Cambridge: Cambridge University Press.