Multilevel Analysis in Stata: A Step-by-Step Guide

1. Introduction

2. Estimating Linear Mixed-effects Models in Stata

3. Useful Resources

We use multilevel or mixed-effects models (also known as hierarchical models) when the data is grouped, structured, or nested in multiple levels. Mixed-effects models consist of fixed effects (coefficients that do not vary by group) and random effects (coefficients that vary by group). 

Examples of grouped or nested data:

• A dataset contains data at patients and doctors levels or patients, doctors, and hospital levels. 
• A dataset contains data at individual and state levels or individual, state, and country levels.
• A dataset contains data at students and classrooms levels or students, classrooms, and schools levels, or students, classrooms, schools, and school-district levels. 

Multilevel models allow:

• Study the effects that vary by entity (or groups)
• Estimate group level averages

Some advantages:

• Regular regression ignores the average variation between entities
• Individual regression may face sample problems and lack of generalization

Pitfalls of ignoring multilevel data structure (Finch et al., 2019, page 29):

• Underestimating standard errors leading to incorrect statistical inference
• Smaller p-values than they really should be
• May miss important relationships involving each level in the data

Stata has the option to estimate varieties of mixed-effects models, including linear mixed-effects models, generalized linear mixed-effects models, survival mixed-effects models, and nonlinear mixed-effects models. See this page for all the available options. 

This tutorial provides step-by-step guides to estimate linear mixed-effects models using Stata.

This section provides a step-by-step guide to conducting multilevel analysis using cross-sectional data*.

First, get the data. Type the following codes on the Stata command window:

use "https://dss.princeton.edu/training/multilevel_example1.dta"

The dataset contains two levels:  individuals (Level 1) and countries (Level 2). Here, individuals at Level 1 are nested within countries at Level 2.

The dataset contains the following variables: 

Identifier variables:

individualID: level 1 identifier.

country: level 2 identifier. The dataset contains data for 20 countries, which indicates we have 20 groups in it.

Level 1 dependent variable:

income: individuals' annual incomes (log, international dollar).

Level 1 predictors:

education: respondents' level of education (binary, assigned 1 if completed secondary education and 0 otherwise).

age:  age of the respondents.

gender: respondents' gender (binary, assigned 1 if female and 0 otherwise).

Level 2 predictors:

HDI: country-level human development index score.

demScore: country-level democracy score.

Let's now estimate the following basic multilevel/mixed effects models. 

 

1. Varying-Intercept Model with no Predictors (Null Model)

In a varying-intercept model, the intercept (constant term) of the regression model is allowed to vary across different groups or clusters. For instance, if we have data on individuals from different countries, a varying-intercept model allows the intercept to vary for each country, capturing the fact that individuals from different countries may have different baseline levels of the outcome variable.

We estimate the varying-intercept model with no predictors (null model) to determine whether our data show evidence of clustering or group effects. If we find evidence that the dataset contains clustering/grouping effects, we should use a multilevel modeling approach. If we do not find any evidence of clustering effects in the dataset, we should estimate single-level models. 

To estimate the null model, type:

mixed income || country:, mle

Here,

mixed: is the Stata command for estimating a multilevel model.

income: is the dependent variable.

|| country: this part specifies the random-effects structure of the model. In this case, country is a grouping variable (here, level 2 identifier variable) for which random intercepts are included.

, mle: specifies that the Maximum Likelihood Estimation (MLE) method should be used for estimating the model parameters.

Stata will give us the following outputs:

. mixed income || country:, mle

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -28299.631  
Iteration 1:  Log likelihood = -28299.631  

Computing standard errors ...

Mixed-effects ML regression                         Number of obs    =  20,125
Group variable: country                             Number of groups =      20
                                                    Obs per group:
                                                                 min =     641
                                                                 avg = 1,006.2
                                                                 max =   2,369
                                                    Wald chi2(0)     =       .
Log likelihood = -28299.631                         Prob > chi2      =       .

------------------------------------------------------------------------------
      income | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |   8.733811   .1817174    48.06   0.000     8.377652    9.089971
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
country: Identity            |
                  var(_cons) |   .6593976   .2088699       .354423    1.226797
-----------------------------+------------------------------------------------
               var(Residual) |   .9685608   .0096603      .9498109    .9876809
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 9408.90       Prob >= chibar2 = 0.0000
 

Interpretation

The above output table includes fixed effects estimates, random effects variance components, and other model diagnostic statistics.

- In the fixed effects part, the coefficient for _cons states that the grand mean of country level intercepts is 8.733811,and this is statistically significant as p-value is 0.00. In other words, the grand mean income of 20 countries in the dataset is 8.73 (log-scaled).

-  In the random effects part, the variance of the random intercepts for the countries (var(_cons)) is 0.6593976 with a standard error of 0.2088699. This suggests substantial variation in income across countries, which is significantly different from zero as the confidence interval does not include zero (0.354423 to 1.226797).

-  The residual variance (var(Residual)) is 0.9685608 with a standard error of 0.0096603, indicating the variability in income within countries that is not explained by the model.

-  The  likelihood ratio test (LR test) tests the null hypothesis of single level model (i.e., H0: var(_cons) = 0) against a multilevel model. The above table gives a LR test chi-square statistic of 9408.90, with a p-value of 0.0000, indicating that the multilevel model provides a significantly better fit to the data than a single level model without random effects.

In sum, the above null mixed-effects model indicates that the average income (intercept) is 8.733811, with significant variability both between and within countries. The inclusion of random intercepts for countries is justified by the significant LR test, highlighting the importance of accounting for country-level effects in the model. 

Although the LR test can suggest whether we should use mixed effects model or a single level model, we can also calculate the intraclass correlation coefficient to determine it. The intraclass correlation coefficient measures the degree of similarity among observations within the same group or cluster.

To obtain the intraclass correlation coefficient, type:

estat icc

Stata will give us the following result:

. estat icc

Intraclass correlation

------------------------------------------------------------------------------
                       Level |        ICC   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
                     country |   .4050457   .0763716      .2678361    .5588903
------------------------------------------------------------------------------

Decision rule:

If the intraclass correlation coefficient (ICC) approaches 0 then the grouping by entities (here counties) are of no use (so you may just run a single level model). If the ICC approaches 1 then there is no variance to explain at the individual level, everybody is the same.

Here the ICC is greater than  0, which is .4050457 suggesting that we should use multilevel model instead of single level model for this dataset. 

Note: the ICC is the ratio of level 2 variance to the sum of level 1 and level 2 variance estimates. We can thus calculate it manually by plugging in the respective values from the null model output table as follows:

ICC = var(_cons)/ var(_cons)+  var(Residual)

    = .6593976 / (.6593976 + .9685608)

    = .4050457

 

2. Varying-intercept Model with Level 1 Predictors

From the LR test above and from the ICC value, we now know that our dataset requires estimating multilevel/mixed effects models. Let's now expand our model by adding all the available level 1 predictors. Run the following codes in Stata:

mixed income education age gender || country:, mle

Stata will give us the following outputs:

. mixed income education age gender || country:, mle

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -27237.919  
Iteration 1:  Log likelihood = -27237.919  (backed up)

Computing standard errors ...

Mixed-effects ML regression                         Number of obs    =  19,914
Group variable: country                             Number of groups =      20
                                                    Obs per group:
                                                                 min =     639
                                                                 avg =   995.7
                                                                 max =   2,348
                                                    Wald chi2(3)     = 1507.06
Log likelihood = -27237.919                         Prob > chi2      =  0.0000

------------------------------------------------------------------------------
      income | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   education |   .5521803    .014848    37.19   0.000     .5230787    .5812818
         age |   .0020658   .0004749     4.35   0.000     .0011351    .0029966
      gender |  -.0905304   .0134999    -6.71   0.000    -.1169896   -.0640712
       _cons |   8.397502   .1773812    47.34   0.000     8.049841    8.745162
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
country: Identity            |
                  var(_cons) |   .6171358   .1954933      .3316968    1.148207
-----------------------------+------------------------------------------------
               var(Residual) |   .8968703   .0089926      .8794173    .9146697
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 9119.20       Prob >= chibar2 = 0.0000

 

Interpretation

The above output table includes fixed effects estimates, random effects variance components, and other model diagnostic statistics.

Fixed Effects Estimates:

• education: Completing secondary education (coded as 1) is associated with an approximate 0.552 unit increase in the logarithm of income compared to those who have not completed secondary education (coded as 0). Given that income is log-transformed, this coefficient corresponds to a 73.7% increase in income, calculated as (e0.552−1)×100%. This indicates that individuals with secondary education earn approximately 73.7% more than their counterparts without this education level, holding other factors constant. This effect is statistically significant at 1% level as the p-value is < 0.0001. This suggests a robust and substantial impact of secondary education on income levels.

• age: Each additional year of age is associated with a 0.0021 unit increase in the logarithm of income, which is statistically significant (p < 0.0001). This indicates that income tends to grow exponentially by about 0.21% per additional year of age, assuming other variables are constant.

• gender: The coefficient for gender is -0.0905, indicating that females (coded as 1) earn less than males (coded as 0), holding all other factors constant. Given that income is log-transformed, this coefficient corresponds to an 8.6% lower income for females compared to males, calculated as (e−0.0905−1)×100%. In other words, females earn approximately 8.6% less than males, holding all other variables constant. This result is statistically significant at the 1% level as p < 0.0001.

Random Effects Estimates:

• var(_cons): The variance of the random intercepts for countries is 0.6171358 with a standard error of  0.1954933 suggests that there is substantial variation in income across countries. This is statistically significant, as the confidence interval ranging from 0.3316968 to 1.148207 does not include zero.

• var(Residual): The residual variance within countries is 0.8969, indicating the variation in income within countries that is not explained by the included variables (education, age, and gender). This high variance suggests that there remains a notable amount of variation in income within each country, even after accounting for three individual-level variables.

Model Diagnostic Statistics:

• Wald Test: The Wald chi-square test tests the null hypothesis that the fixed effects are simultaneously zero. The chi-square statistic of 1507.06 with 3 degrees of freedom and a p-value of 0.0000 suggests that the model's predictors (i.e., education, age, and gender) significantly help in explaining the variability in income.

• Likelihood Ratio Test: The chi-square statistic of 9119.20 with a p-value of 0.0000 when comparing against a single level model indicates that incorporating random effects (variance across countries) significantly improves the model's fit.

 

3. Varying-intercept Model with Level 1 and Level 2 Predictors

Let us now add the level 2 variables (HDI and demScore) in our model. Use the following codes:

mixed income education age gender HDI demScore || country:, mle

Stata will give us the following outputs:

. mixed income education age gender HDI demScore || country:, mle

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -27228.959  
Iteration 1:  Log likelihood = -27228.959  (backed up)

Computing standard errors ...

Mixed-effects ML regression                         Number of obs    =  19,914
Group variable: country                             Number of groups =      20
                                                    Obs per group:
                                                                 min =     639
                                                                 avg =   995.7
                                                                 max =   2,348
                                                    Wald chi2(5)     = 1538.21
Log likelihood = -27228.959                         Prob > chi2      =  0.0000

------------------------------------------------------------------------------
      income | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
   education |   .5514025   .0148477    37.14   0.000     .5223014    .5805035
         age |   .0020463   .0004749     4.31   0.000     .0011155    .0029771
      gender |  -.0906056   .0134998    -6.71   0.000    -.1170647   -.0641464
         HDI |   6.225759   1.163167     5.35   0.000     3.945994    8.505525
    demScore |  -.0771922   .1001299    -0.77   0.441    -.2734431    .1190587
       _cons |   4.562724   .8916805     5.12   0.000     2.815063    6.310386
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
country: Identity            |
                  var(_cons) |   .2513366   .0798002      .1348952    .4682901
-----------------------------+------------------------------------------------
               var(Residual) |   .8968704   .0089926      .8794173    .9146698
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 4486.86       Prob >= chibar2 = 0.0000

 

Interpretation

The above output table includes fixed effects estimates, random effects variance components, and other model diagnostic statistics.

Fixed Effects Estimates:

• education: Completing secondary education (coded as 1) is associated with an approximate 0.5514 unit increase in the logarithm of income compared to those who have not completed secondary education (coded as 0). Given that income is log-transformed, this coefficient corresponds to a 73.57% increase in income, calculated as (e0.5514−1)×100%.This indicates that individuals with secondary education earn approximately 73.57% more than their counterparts without this education level, holding other factors constant. This effect is statistically significant at 1% level as the p-value is < 0.0001. This suggests a robust and substantial impact of secondary education on income levels.

• age: Each additional year of age is associated with a 0.0020 unit increase in the logarithm of income, which is statistically significant (p < 0.0001). This indicates that income tends to grow exponentially by about 0.20% per additional year of age, assuming other variables are constant.

• gender: The coefficient for gender is -0.0906, indicating that females (coded as 1) earn less than males (coded as 0), holding all other factors constant. Given that income is log-transformed, this coefficient corresponds to an 8.66% lower income for females compared to males, calculated as (e−0.0906−1)×100%. In other words, females earn approximately 8.66% less than males, holding all other variables constant. This result is statistically significant at the 1% level as p < 0.0001.

• HDI: Each unit increase in HDI is associated with a substantial increase in the log of income by 6.2258, which is highly significant (p < 0.0001). This indicates a strong positive relationship between country human development status and individual income.

• demScore: The democracy score has a coefficient of -0.0772, but this is not statistically significant (p = 0.441), suggesting that within the context of this model, democracy score does not have a clear impact on income.

Random Effects Estimates:

• var(_cons): The variance of the random intercepts for countries is 0.2513 with a standard error of  0.0798 suggests that there is substantial variation in income across countries. This is statistically significant, as the confidence interval ranging from 0.1348952 to 0.4682901 does not include zero.

• var(Residual): The within-country variance in income after accounting for the fixed effects is 0.8969, indicating substantial unexplained variability within countries.

Model Diagnostic Statistics:

• Wald Test: The Wald chi-square test tests the null hypothesis that the fixed effects are simultaneously zero. The chi-square statistic of 1538.21 with 5 degrees of freedom and a p-value of 0.0000 suggests that the model's predictors (i.e., education, age, gender, HDI, and demScore) significantly help in explaining the variability in income.

• Likelihood Ratio Test:  The likelihood ratio (LR) test against a linear model without random effects yields a chi-square value of 4486.86 with a p-value of 0.0000, indicating that the mixed-effects model fits significantly better than a simpler linear model without random effects.

 

Postestimation: Comparing Models using Likelihood-ratio Test

We can use the likelihood-ratio test (lrtest) to compare the model with level 1 predictors with the model with both level 1 and level 2 predictors. This test compares the “log likelihood” (as shown in the output table below) of two models and tests whether they are significantly different. Use the following Stata codes:

/* estimate the model with level 1 predictors and store the estimates */
mixed income education age gender || country:, mle
estimates store m1 
/* estiamte the model with both level 1 and level 2 predictors and store the estimates */
mixed income education age gender HDI demScore || country:, mle
estimates store m2 
/* run the likelihood-ratio test to compare*/
lrtest m1 m2

Stata will give us the following outputs:

. lrtest m1 m2

Likelihood-ratio test
Assumption: m1 nested within m2

 LR chi2(2) =  17.92
Prob > chi2 = 0.0001

The null hypothesis is that there is no significant difference between the two models. If Prob > chi2 <0.05, then you may reject the null and conclude that there is a statistically significant difference between the models. In the example above, we reject the null and conclude that the  model with both level 1 and level 2 predictors provides a better fit (as it has the lowest log likelihood). 

 

4. Varying-Intercept, Varying-Slope Model 

The Varying-Intercept, Varying-Slope Model (also known as the Varying-Coefficient, Varying-Intercept Model or simply the Random Slopes Model) is an extension of the basic mixed effects model framework. In a varying-intercept, varying-slope model, both intercepts and coefficients vary across groups simultaneously.  If we want to covariate intercepts and slopes for only two of the level 1 predictors (education and age), use the following codes in Stata. Depending on your research questions or conceptual framework, you can select more or fewer predictors to covary (we have selected two predictors for demonstration purposes only). 

mixed income education age gender HDI demScore || country: education age, covariance(unstructured) mle

Here, 

mixed: the command in Stata for fitting mixed effects models.

income: dependent variable.

education age gender HDI demScore: fixed effects predictors (both level 1 and level 2).

|| country: education age : specifies the random effects part of the model, where

• The || syntax denotes that random effects are being specified for the grouping variable.
• country is the grouping variable (level-2 cluster) where random effects are specified.
• education and age are the variables for which random slopes are included. This means that the model allows the relationship between education and income, and age and income to vary across different countries.

covariance(unstructured): specifies the covariance structure for the random effects.

• unstructured means that the model allows for an unconstrained covariance matrix for the random effects. In other words, the variances and covariances among the random effects (for education and age in this case) are estimated freely without any restrictions.

, mle: specifies that the model should be estimated using Maximum Likelihood Estimation (MLE). This is the default estimation method for mixed effects models in Stata.

Running the above codes will give us the following outputs:

. mixed income education age gender HDI demScore || country: education age, covariance(unstructured) mle

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -27114.889  
Iteration 1:  Log likelihood = -27114.889  

Computing standard errors ...

Mixed-effects ML regression                          Number of obs    = 19,914
Group variable: country                              Number of groups =     20
                                                     Obs per group:
                                                                  min =    639
                                                                  avg =  995.7
                                                                  max =  2,348
                                                     Wald chi2(5)     = 192.76
Log likelihood = -27114.889                          Prob > chi2      = 0.0000

--------------------------------------------------------------------------------------
              income | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
---------------------+----------------------------------------------------------------
           education |    .550412   .0547711    10.05   0.000     .4430626    .6577614
                 age |   .0018581   .0010787     1.72   0.085    -.0002561    .0039722
              gender |  -.0914229   .0134394    -6.80   0.000    -.1177637   -.0650822
                 HDI |   6.193788   1.140026     5.43   0.000     3.959379    8.428197
            demScore |  -.0611939   .0979584    -0.62   0.532    -.2531888     .130801
               _cons |   4.507244   .8770653     5.14   0.000     2.788228    6.226261
--------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
country: Unstructured        |
              var(education) |   .0551994   .0189224       .028193    .1080754
                    var(age) |   .0000184   7.30e-06      8.47e-06      .00004
                  var(_cons) |   .3722062   .1227356      .1950287     .710344
          cov(education,age) |    .000543   .0002936     -.0000324    .0011185
        cov(education,_cons) |  -.0746274   .0388116     -.1506966    .0014419
              cov(age,_cons) |  -.0013804   .0007623     -.0028746    .0001138
-----------------------------+------------------------------------------------
               var(Residual) |   .8831506   .0088638      .8659477    .9006954
------------------------------------------------------------------------------
LR test vs. linear model: chi2(6) = 4715.00               Prob > chi2 = 0.0000

 

Interpretation

The above output table includes fixed effects estimates, random effects variance-covariance components, and other model diagnostic statistics.

Fixed Effects Estimates:

• education: Completing secondary education (coded as 1) is associated with an approximate 0.5504 unit increase in the logarithm of income compared to those who have not completed secondary education (coded as 0). Given that income is log-transformed, this coefficient corresponds to a 73.0% increase in income, calculated as (e0.5504−1)×100%. This indicates that individuals with secondary education earn approximately 73.0% more than their counterparts without this education level, holding other factors constant. This effect is statistically significant at a 1% level as the p-value is < 0.0001. This suggests a robust and substantial impact of secondary education on income levels.

• age: Each additional year of age is associated with a 0.0019 unit increase in the logarithm of income, which is not statistically significant at the conventional 5% level, as the p-value is 0.085. 

• gender: The coefficient for gender is -0.0914, indicating that females (coded as 1) earn less than males (coded as 0), holding all other factors constant. Given that income is log-transformed, this coefficient corresponds to an 8.7% lower income for females compared to males, calculated as (e−-0.0914−1)×100%. In other words, females earn approximately 8.7% less than males, holding all other variables constant. This result is statistically significant at the 1% level as p < 0.0001.

• HDI: Each unit increase in HDI is associated with a substantial increase in the log of income by 6.1938, which is highly significant (p < 0.0001). This indicates a strong positive relationship between country human development status and individual income.

• demScore: The democracy score has a coefficient of -0.0612, but this is not statistically significant (p = 0.532), suggesting that within the context of this model, democracy score does not have a clear impact on income.

Random Effects Estimates:

• var(education): The variance of the random slope for education indicates the extent to which the effect of education on income varies across countries. The positive variance estimate (0.0552) suggests significant heterogeneity in how education impacts income between countries. The confidence interval shows that this variability is statistically significant, as it does not include zero.

• var(age): The variance of the random slope for age is very small (0.00002), indicating that there is minimal variability in the effect of age on income across countries. The confidence interval suggests that this variance is statistically significant, but the small magnitude implies that age has a relatively stable effect on income regardless of the country.

• var(_cons): The variance of the random intercept reflects the variability in baseline income levels across countries after accounting for the fixed effects. A significant intercept variance (0.3722) indicates that there are substantial differences in the average income levels across countries, suggesting that country-specific factors influence baseline income levels.

• cov(education,age):The positive covariance estimate (0.00054) suggests a slight positive relationship between the random slopes for education and age. This means that in countries where the effect of education on income is stronger, the effect of age on income is also likely to be stronger. However, the confidence interval includes zero, indicating that this relationship is not statistically significant. 

• cov(education,_cons): The negative covariance (-0.0746) suggests that in countries with higher baseline income levels, the impact of education on income tends to be weaker. However, the confidence interval includes zero, which implies that this relationship is not statistically significant.

• var(Residual): The residual variance represents the variability in income that is not explained by the model, after accounting for both fixed and random effects. A residual variance of 0.8832 indicates that there is still substantial variability in income that is not captured by the predictors and random effects in the model.

Model Diagnostic Statistics:

• Wald Test: The Wald chi-square test tests the null hypothesis that the fixed effects are simultaneously zero. The chi-square statistic of 192.76 with 5 degrees of freedom and a p-value of 0.0000 suggests that the model's predictors (i.e., education, age, gender, HDI, and demScore) significantly help in explaining the variability in income.

• Likelihood Ratio Test:  The likelihood ratio (LR) test against a linear model without random effects yields a chi-square value of 4715.00 with a p-value of 0.0000, indicating that the mixed-effects model fits significantly better than a simpler linear model without random effects.

* The dataset is obtained from Constantelos et al. (2023), available in Harvard Dataverse, and modified for demonstration purposes.