4  Cross Sectional Multilevel Models

4.1 The Equation

Recall the general model of Equation 1.1, and the syntax outlined in Section 1.3. Below in Equation 4.1, we consider a more substantive example.

\[\text{outcome}_{ij}= \beta_0 + \beta_1 \text{warmth}_{ij} + \tag{4.1}\]

\[\beta_2 \text{physical punishment}_{ij} +\]

\[\beta_3 \text{identity}_{ij} + \beta_4 \text{intervention}_{ij} + \beta_5 \text{HDI}_{ij} +\]

\[u_{0j} + u_{1j} \times \text{warmth}_{ij} + e_{ij}\]

4.2 Correlated and Uncorrelated Random Effects

Consider the covariance matrix of random effects (e.g. \(u_{0j}\) and \(u_{1j}\)). In Equation 4.2 the covariances of the random effects are constrained to be zero.

\[\begin{bmatrix} var(u_{0j}) & 0 \\ 0 & var(u_{1j}) \end{bmatrix} \tag{4.2}\]

As discussed in the Chapter on multilevel models with cross-sectional data, however, one can consider a multilevel model in which the random effects are correlated, as is the case in Equation 4.3.

\[\begin{bmatrix} var(u_{0j}) & cov(u_{0j}, u_{1j}) \\ cov(u_{0j}, u_{1j}) & var(u_{1j}) \end{bmatrix} \tag{4.3}\]

Procedures for estimating models with uncorrelated and correlated random effects are detailed below (Bates et al., 2015; Bates, 2024; StataCorp, 2021).

Table 4.1: Correlated and Uncorrelated Random Effects

Software	Uncorrelated Random Effects	Correlated Random Effects
Stata	default	add option: , cov(uns)
R	separate random effects from grouping variable with ||	separate random effects from grouping variable with |
Julia	separate terms for each random effect e.g. (1 | group) + (0 + x | group)	separate random effects from grouping variable with |.

All models in the examples below are run with uncorrelated random effects, but could just as easily be run with correlated random effects.

4.3 Run Models

Stata

4.3.0.1 Get The Data

use simulated_multilevel_data.dta

4.3.0.2 Run The Model

mixed outcome warmth physical_punishment i.identity i.intervention HDI || country: warmth

Performing EM optimization ...

Performing gradient-based optimization: 
Iteration 0:  Log likelihood = -9626.6279  
Iteration 1:  Log likelihood =  -9626.607  
Iteration 2:  Log likelihood =  -9626.607  

Computing standard errors ...

Mixed-effects ML regression                          Number of obs    =  3,000
Group variable: country                              Number of groups =     30
                                                     Obs per group:
                                                                  min =    100
                                                                  avg =  100.0
                                                                  max =    100
                                                     Wald chi2(5)     = 334.14
Log likelihood =  -9626.607                          Prob > chi2      = 0.0000

-------------------------------------------------------------------------------------
            outcome | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
--------------------+----------------------------------------------------------------
             warmth |   .8345368   .0637213    13.10   0.000     .7096453    .9594282
physical_punishment |  -.9916657   .0797906   -12.43   0.000    -1.148052   -.8352791
         2.identity |  -.3004767   .2170295    -1.38   0.166    -.7258466    .1248933
     1.intervention |   .6396427   .2174519     2.94   0.003     .2134448    1.065841
                HDI |   -.003228   .0199257    -0.16   0.871    -.0422817    .0358256
              _cons |   51.99991   1.371257    37.92   0.000      49.3123    54.68753
-------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects parameters  |   Estimate   Std. err.     [95% conf. interval]
-----------------------------+------------------------------------------------
country: Independent         |
                 var(warmth) |   .0227504   .0257784      .0024689    .2096436
                  var(_cons) |   2.963975   .9737647      1.556777    5.643163
-----------------------------+------------------------------------------------
               var(Residual) |   34.97499   .9097109      33.23668    36.80422
------------------------------------------------------------------------------
LR test vs. linear model: chi2(2) = 205.74                Prob > chi2 = 0.0000

Note: LR test is conservative and provided only for reference.

R

4.3.0.3 Get The Data

library(haven)

df <- read_dta("simulated_multilevel_data.dta")

4.3.0.4 Run The Model

Caution

lme4 does not directly provide p values in results, because of some disagreement over exactly how these p values should be calculated. Therefore, in this Appendix, I also call library lmerTest to provide p values for lme4 results.

Tip

R prefers to use scientific notation when possible. I find that the use of scientific notation can be confusing in reading results. I turn off scientific notation by setting a penalty for its use: options(scipen = 999).

library(lme4) 

library(lmerTest)

options(scipen = 999) 

fit1 <- lmer(outcome ~ warmth + physical_punishment + 
               identity + intervention + HDI +
               (1 + warmth || country),
             data = df)

summary(fit1)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: outcome ~ warmth + physical_punishment + identity + intervention +  
    HDI + (1 + warmth || country)
   Data: df

REML criterion at convergence: 19268.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.9774 -0.6563  0.0187  0.6645  3.6730 

Random effects:
 Groups    Name        Variance Std.Dev.
 country   (Intercept)  3.19056 1.786   
 country.1 warmth       0.02465 0.157   
 Residual              35.01782 5.918   
Number of obs: 3000, groups:  country, 30

Fixed effects:
                       Estimate  Std. Error          df t value
(Intercept)           52.311714    1.446735   33.113738  36.158
warmth                 0.834562    0.064252   41.896966  12.989
physical_punishment   -0.991892    0.079845 2968.010901 -12.423
identity              -0.300350    0.217179 2970.106304  -1.383
intervention           0.639059    0.217603 2971.185215   2.937
HDI                   -0.003395    0.020596   27.598517  -0.165
                                Pr(>|t|)    
(Intercept)         < 0.0000000000000002 ***
warmth              0.000000000000000277 ***
physical_punishment < 0.0000000000000002 ***
identity                         0.16678    
intervention                     0.00334 ** 
HDI                              0.87027    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) warmth physc_ idntty intrvn
warmth      -0.119                            
physcl_pnsh -0.145 -0.003                     
identity    -0.220 -0.012 -0.003              
interventin -0.077  0.034  0.022 -0.018       
HDI         -0.922 -0.006  0.009 -0.001  0.000

Julia

4.3.0.5 Get The Data

using Tables, MixedModels, StatFiles, DataFrames, CategoricalArrays, DataFramesMeta

df = DataFrame(load("simulated_multilevel_data.dta"))

4.3.0.6 Change Country To Categorical

@transform!(df, :country = categorical(:country))

4.3.0.7 Run The Model

m1 = fit(MixedModel, @formula(outcome ~ warmth + physical_punishment + 
               identity + intervention + HDI +
               (1 | country) +
               (0 + warmth | country)), df)

Linear mixed model fit by maximum likelihood
 outcome ~ 1 + warmth + physical_punishment + identity + intervention + HDI + (1 | country) + (0 + warmth | country)
   logLik   -2 logLik     AIC       AICc        BIC    
 -9626.6070 19253.2140 19271.2140 19271.2742 19325.2713

Variance components:
            Column    Variance Std.Dev.   Corr.
country  (Intercept)   2.963849 1.721583
         warmth        0.022756 0.150852   .  
Residual              34.974984 5.913965
 Number of obs: 3000; levels of grouping factors: 30

  Fixed-effects parameters:
─────────────────────────────────────────────────────────────
                          Coef.  Std. Error       z  Pr(>|z|)
─────────────────────────────────────────────────────────────
(Intercept)          52.3004      1.40406     37.25    <1e-99
warmth                0.834537    0.0637228   13.10    <1e-38
physical_punishment  -0.991665    0.0797906  -12.43    <1e-34
identity             -0.300475    0.217029    -1.38    0.1662
intervention          0.639641    0.217452     2.94    0.0033
HDI                  -0.0032286   0.0199255   -0.16    0.8713
─────────────────────────────────────────────────────────────

4.4 Interpretation

Models suggest that parental warmth is associated with increases in the beneficial outcome, while physical punishment is associated with decreases in the beneficial outcome. Membership in the group represented by identity is not associated with the outcome. The intervention is associated with increases in the outcome. The Human Development Index is not associated with the outcome.

Bates, D. (2024). MixedModels.jl Documentation. https://juliastats.org/MixedModels.jl/stable/

Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting linear mixed-effects models using lme4. Journal of Statistical Software, 67(1), 1–48. https://doi.org/10.18637/jss.v067.i01

StataCorp. (2021). Stata 17 multilevel mixed effects reference manual. Stata Press.