# Why Does Standard Error for Individual Predictors Not Increase in Multilevel Modeling?

I have been trying to teach myself multilevel modeling through R and I am relying on the W. Holmes Finch book. According to them, and any resources in general, not doing multilevel modeling causes an underestimation of standard error (since the DV will be correlated for individuals within a group). To test this I used the Demo.twolevel package that comes with lavaan and tested two models using lm and lmer function. I show the code and outputs below but the interesting obervation is that w1 and w2, which are group level predictors, do get a higher standard error once I cluster them in a random intercept model. But the individual level predictors such as x1,x2, x3 have a lower standard error once they are clustered. My question is why is there such a mismatch? Should individual level predictors not also get a higher standard error? Presumably they should be correlated at a group level so once that is accounted for the estimated error should increase.

Multilevel model code: summary(lmer(fw~x1+x2+x3+w1+w2+(1|cluster), data=Demo.twolevel))

Linear mixed model fit by REML ['lmerMod']
Formula: fw ~ x1 + x2 + x3 + w1 + w2 + (1 | cluster)
Data: Demo.twolevel

REML criterion at convergence: 6762.6

Scaled residuals:
Min      1Q  Median      3Q     Max
-3.6210 -0.6593 -0.0160  0.6561  3.0715

Random effects:
Groups   Name        Variance Std.Dev.
cluster  (Intercept) 0.5750   0.7583
Residual             0.7216   0.8494
Number of obs: 2500, groups:  cluster, 200

Fixed effects:
Estimate Std. Error t value
(Intercept) -0.01120    0.05721  -0.196
x1           0.42797    0.01774  24.130
x2           0.34007    0.01750  19.431
x3           0.17029    0.01735   9.818
w1           0.12436    0.06148   2.023
w2           0.11439    0.05964   1.918

Correlation of Fixed Effects:
(Intr) x1     x2     x3     w1
x1  0.001
x2  0.002  0.007
x3 -0.002  0.012  0.002
w1 -0.038 -0.007 -0.013  0.010
w2  0.089 -0.005  0.005  0.012  0.144


summary(lm(fw~x1+x2+x3+w1+w2, data=Demo.twolevel))

Call:
lm(formula = fw ~ x1 + x2 + x3 + w1 + w2, data = Demo.twolevel)

Residuals:
Min      1Q  Median      3Q     Max
-4.4237 -0.7846  0.0190  0.7842  3.4769

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.002431   0.022836   0.106    0.915
x1          0.423292   0.022942  18.451  < 2e-16 ***
x2          0.312961   0.022617  13.837  < 2e-16 ***
x3          0.141685   0.022265   6.364 2.34e-10 ***
w1          0.165201   0.024126   6.847 9.44e-12 ***
w2          0.133400   0.023902   5.581 2.65e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.137 on 2494 degrees of freedom
Multiple R-squared:  0.2064,    Adjusted R-squared:  0.2048
F-statistic: 129.7 on 5 and 2494 DF,  p-value: < 2.2e-16


"Not doing multilevel modeling causes an underestimation of standard error" is an imprecise statement as it doesn't specify which standard error we underestimate if we ignore the hierarchical structure.

In fact, it's an incorrect statement. The residual (so called unexplained) variation is lower when we fit a better model; so we overestimate the residual standard error if we fit a single-level model when we should fit a multi-level model instead.

To understand this better, it's sufficient to consider one individual-level predictor and one group-level predictor. Note: R code to reproduce the analysis is attached at the end.

The single-level model is a linear regression with two predictors, $$x$$ and $$w$$, and three regression coefficients, $$\alpha_0,\beta_x,\beta_w$$.

\begin{aligned} \text{(M1)} \quad\quad \operatorname{y}_i &\sim N \left(\alpha_0 + \beta_x\operatorname{x}_i + \beta_w\operatorname{w}_{c[i]}, \sigma^2 \right) \end{aligned}

I write $$\operatorname{w}_{c[i]}$$ rather than $$\operatorname{w}_i$$ because every individual in cluster $$c$$ has the same $$w$$ value.

# single-level model
M1 <- lm(
y ~ x + w,
data = Demo.twolevel
)
# regression coefficients
#>   term        estimate std.error
#>   (Intercept)   0.0245    0.0329
#>   x             0.495     0.0332
#>   w             0.198     0.0344
# residual standard error
sigma(M1)
#> 1.646398


The multi-level model is a random-intercept model and has two components. Each cluster $$c$$ has an (unobserved) random intercept $$\alpha_c$$ which is modeled with a simple linear regression on the group-level predictor $$w$$. Each observation $$y_i$$ is modeled with a simple linear regression on the random intercept $$\alpha_{c[i]}$$ and the individual-level predictor $$x$$.

\begin{aligned} \text{(M2)} \quad\quad \operatorname{y}_i &\sim N \left(\alpha_{c[i]} + \beta_x\operatorname{x}_i, \sigma^2 \right) \\ \alpha_c &\sim N \left(\alpha_0 + \beta_w\operatorname{w}_c, \sigma^2_{\alpha} \right) \end{aligned}

# multi-level model
M2 <- lmer(
y ~ x + w + (1 | cluster),
data = Demo.twolevel
)
# regression coefficients (fixed effects)
#>   term        estimate std.error
#>   (Intercept)   0.0148    0.0742
#>   x             0.494     0.0277
#>   w             0.161     0.0791
# residual standard error
sigma(M2)
#> 1.326981


We observe:

• The residual standard error is smaller in the hierarchical model.
• The standard error of the individual-level predictor x is lower in the hierarchical model M2.
• The standard error of the group-level predictor w is lower in the single-level model M1.

Check that these statements also hold for the individual-level predictors x1, x2, x3 and the group-level predictors w1, w2 in the bigger models in your original example. There is no loss of generality in considering a simpler model; but it makes plotting the data easier.

Now let's show that the hierarchical model does a better job at explaining the differences between clusters, by computing the "remainder" of $$y$$ after accounting for $$w$$. From the single-level model M1 we compute $$y_i - (\hat{\alpha}_0 + \hat{\beta}_ww_{c[i]})$$; and from the multi-level model M2 we compute $$y_i - \hat{\alpha}_{c[i]}$$.

In the figure below I plot the first 16 clusters in separate panels and add two horizontal lines: For each cluster $$c$$, the "fixed intercept" $$(\hat{\alpha}_0 + \hat{\beta}_ww_c)$$ is the gray solid line; the random intercept $$\hat{\alpha}_{c[i]}$$ is the black dotted line. It's a bit hard to see but the gray lines don't intercept the y axis at exactly the same value. However, there is too little variability to account for the observed differences between clusters. Visually, the dotted lines are closer to the average $$y$$ in each panel.

Once we subtract the "explanation" due to $$w$$, we plot the remainder of $$y$$ as a function of $$x$$ (for 16 clusters only). There is the same number of points in the two side-by-side plots and they have the same $$x$$ values. However, the range and variability of the $$y$$ remainder is smaller in the hierarchical model. This shows right away that M2 has smaller residual variance $$\hat{\sigma}^2$$ than M1. Since the standard errors of regression coefficients are a function of $$\hat{\sigma}^2$$, it also shows why the individual-level coefficient $$\beta_x$$ has smaller standard error in M2 even though $$x$$ is the same in M1 and M2.

On the other hand, the group-level coefficient $$\beta_w$$ has bigger standard error in the hierarchical model M2. Intuitively, this happens because M2 estimates $$\beta_w$$ from a smaller regression of the random intercepts $$\alpha_c$$ on $$w$$ while M1 estimates $$\beta_w$$ from the bigger regression of the observations $$y_i$$ on $$w$$ (and $$x$$ of course). There are fewer clusters than observations, so M2 has more uncertainty about $$\beta_w$$. M1 is too optimistic about its estimate (it underestimates the variance of $$\beta_w$$) as it assumes — incorrectly — that there are $$n$$ independent observations.

The R code in all its gory details:


library("broom.mixed")
library("lme4")
library("tidyverse")

set.seed(1234)

data(Demo.twolevel, package = "lavaan")

# data for analysis
Demo.twolevel <- Demo.twolevel %>%
as_tibble() %>%
select(
y = y1,
x = x1,
w = w1,
cluster
) %>%
mutate(
across(cluster, as.character)
)
# subset for scatterplot
Demo.twolevel.c16 <- Demo.twolevel %>%
filter(
cluster %in% as.character(seq(16))
)

# single-level model
m1 <- lm(
y ~ x + w,
data = Demo.twolevel
)
# regression coefficients
tidy(m1) %>% select(term, estimate, std.error)
# residual standard error
sigma(m1)

# multi-level model
m2 <- lmer(
y ~ x + w + (1 | cluster),
data = Demo.twolevel
)
# regression coefficients (fixed effects)
tidy(m2, effects = "fixed") %>% select(term, estimate, std.error)
# residual standard error
sigma(m2)

demo_scatterplot <- function(points, xvar, yvar) {
points %>%
ggplot(
aes({{ xvar }}, {{ yvar }})
) +
geom_point(
shape = 1
) +
scale_x_continuous(
limits = c(-2.5, 3)
) +
scale_y_continuous(
limits = c(-5, 5.5)
)
}

# Compute (alpha0 + beta1 * w) from the single-level model M1
Demo.twolevel.c16 <-
Demo.twolevel.c16 %>%
mutate(
fixed_intercept = coef(m1)["(Intercept)"] + coef(m1)["w"] * w
)
# Extract alpha_c from the multi-level model M2
Demo.twolevel.c16 <-
tidy(m2, effects = "ran_vals") %>%
select(
cluster = level,
random_intercept = estimate
) %>%
inner_join(
Demo.twolevel.c16,
by = c("cluster")
)

p <-
Demo.twolevel.c16 %>%
demo_scatterplot(x, y - fixed_intercept)
p +
labs(
y = expression(paste(
"y - (", β[0], "-", β[1], "w)"
)),
title = expression(paste(
"Remainder ", italic(y), " variance in M1"
))
)

p <-
Demo.twolevel.c16 %>%
demo_scatterplot(x, y)
p +
geom_hline(
aes(yintercept = fixed_intercept),
linetype = 1,
color = "gray"
) +
geom_hline(
aes(yintercept = random_intercept),
linetype = 3
) +
facet_wrap(
~cluster
) +
labs(
title = expression(paste(
"Variance in ", italic(y), " within clusters 1 to 16, ",
"after accounting for ", italic(w)
))
)

p <-
Demo.twolevel.c16 %>%
demo_scatterplot(x, y - random_intercept)
p +
labs(
y = expression(paste("y - ", α[c])),
title = expression(paste(
"Remainder ", italic(y), " variance in M2"
))
)