Are covariates in a mixed model estimated between-group, within-group, or somewhere in-between? And why?

Question

Consider a linear (perhaps mixed) model where we have one factor $p$, one covariate $x$, and no interaction. I want to understand how the coefficient estimate $\hat\beta$ of $x$, and its standard error, are affected by the decision to make $p$ random or fixed. That is, how the factor's type affects the inference on the covariate.

Intuitively, it feels to me that a model with a fixed $p$ should give the same $\hat\beta$ as one with random $p$. My reasoning for this is that with no interaction, we are declaring that $\beta$ should be the same across all levels of $p$, so it shouldn't matter how many levels $p$ has, nor whether it is random or fixed.

Playing with toy examples, I have found this to be the confirmed in all but cases with small datasets and correlation between $x$ and $p$.

More specifically: for large datasets, random and fixed $p$ give the same $\hat\beta$. For tiny datasets, a model with random $p$ gives the same $\hat\beta$ as one where $p$ is left out altogether (call this the "pooled model"). For some boundary cases in-between, $\hat\beta$ cones out part-way between the two.

Mathematically, I'm totally happy with why these exceptions happen: because there is a likelihood cost to magnifying the effects of the random factor, and not for the fixed factor. When we don't have enough evidence to support expanding the random factor, it is as if it weren't there. But intuitively, as to how to understand $\hat\beta$ in these cases, I'm stumped. In some cases $\hat\beta$ seems to represent the within-group slope, in others the between-group slope, and in others it is somewhere in-between. It has echoes of partial pooling, but I don't think the same logic applies here. This realisation has undermined my confidence in interpreting a specific real-world example I'm dealing with.

The following code and graph gives an example of one of these boundary cases, which I find particularly hard to justify. Lines represent $\hat\beta$ in the fixed (magenta), pooled (blue) and mixed (black) models, and the mixed model gives us a $\hat\beta$ in no-man's land in-between the other two:

---

Code to generate Graph

library(lme4)
set.seed(1)  

beta = -1
sigma = 5

nppt = 5
nperppt = 10

generatedataforppt = function(ppt,n){
  x = rnorm(n)+ppt
  y = beta*x + 3*ppt + rnorm(n,sd=sigma)
  cbind(x=x,y=y,p=ppt)
}

data = as.data.frame(do.call(rbind,lapply(1:nppt,generatedataforppt,nperppt)))
plot(y~x,data=data,pch=data$p)
data$p = factor(data$p)

pooled.model = lm(  y~x,       data=data)
fixed.model  = lm(  y~x+p-1,   data=data)
mixed.model  = lmer(y~x+(1|p), data=data)

abline(coefficients(pooled.model)[1],    coefficients(pooled.model)[2],    col="blue",    lwd=3)
abline(coefficients(fixed.model)[2],     coefficients(fixed.model)[1],     col="magenta", lwd=3)
abline(coefficients(mixed.model)$p[1,1], coefficients(mixed.model)$p[1,2], col="black",   lwd=3)

Answer 1

This is a very nice demonstration of a well-known issue in multilevel modeling, which is that an uncentered predictor that itself has variance across level 1 and level 2 (i.e., has both within and between cluster variation), is some (very often) uninterpretable blend of the within and between effect of that predictor on the outcome. Accordingly, many people who use multilevel models prefer to use some sort of technique to center level 1 variables to make it more clear what the coefficient represents. One of the best papers I know on this issue is by Enders and Tofighi (2007). I strongly suggest you read this paper. Another very nice paper on this topic is a recent psychological methods paper by McNeish and Kelley (2018).

To get a better sense of what is going on, I took your code and created a few new variables.

First, look at the intraclass correlation coefficient for your predictor, x.

library(sjstats)

var.comp.x  = lmer(x ~ 1 + (1|p), data=data)
icc(var.comp.x) # ICC = .7629

Fully 76% of the variance in x sits at the between level, here p. That is a lot in my experience.

When you ran your mixed model with the uncentered x as the predictor, the output was as follows:

> summary(mixed.model)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | p)
   Data: data

REML criterion at convergence: 302.1

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.63592 -0.74081  0.02644  0.66501  1.78631 

Random effects:
 Groups   Name        Variance Std.Dev.
 p        (Intercept) 25.84    5.083   
 Residual             21.71    4.660   
Number of obs: 50, groups:  p, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   8.4194     2.9711   2.834
x            -1.1126     0.6242  -1.782

Correlation of Fixed Effects:
  (Intr)
x -0.604

A one unit change in x is associated with a -1.11 unit change in y. But where is this change happening - moreso at the within or between level?

One way to focus purely on the within level effect is to run a fixed effects model. This tells you the within-cluster effect of a 1-unit change in x:

> summary(fixed.model)

Call:
lm(formula = y ~ x + p - 1, data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.9012 -3.5012 -0.0116  2.9083  8.7211 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
x    -1.683      0.688  -2.446 0.018509 *  
p1    2.604      1.601   1.626 0.111020    
p2    6.329      1.765   3.585 0.000840 ***
p3    8.798      2.545   3.457 0.001224 ** 
p4   15.727      2.778   5.662 1.06e-06 ***
p5   16.844      4.097   4.111 0.000169 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.633 on 44 degrees of freedom
Multiple R-squared:  0.6772,    Adjusted R-squared:  0.6331 
F-statistic: 15.38 on 6 and 44 DF,  p-value: 2.081e-09

So within clusters, a 1-unit difference in x is associated with -1.68 in y. But what about at the between-cluster level? For this, you are best off centering your x variable. If you center each individual's x value around their group's mean for x, which I label CWC (centering within clusters), then you get the exact same value for x as you get in the fixed effects approach. If you add to the model the group mean for x, you get the pure between cluster association between x and y. First we have to create a couple of variables:

library(dplyr)
data <- data %>%
  group_by(p) %>% 
  mutate(group_mn_x = mean(x))
data <- data %>% 
  mutate(CWC_x = x-group_mn_x)

Now let's run a mixed model regressing y on the group-mean centered x and the group mean of x:

> mixed.model.cwc  = lmer(y ~ CWC_x + group_mn_x + (1|p), data=data)
> summary(mixed.model.cwc)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ CWC_x + group_mn_x + (1 | p)
   Data: data

REML criterion at convergence: 294

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.57439 -0.75814 -0.08622  0.68231  1.91484 

Random effects:
 Groups   Name        Variance Std.Dev.
 p        (Intercept)  5.715   2.391   
 Residual             21.468   4.633   
Number of obs: 50, groups:  p, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)   1.2423     2.5371   0.490
CWC_x        -1.6830     0.6880  -2.446
group_mn_x    1.3820     0.7666   1.803

Correlation of Fixed Effects:
           (Intr) CWC_x 
CWC_x       0.000       
group_mn_x -0.869  0.000

The value for CWC_x is the same as the fixed effects, as expected, but the value for the between effect of x on y is 1.382. At the group level, a 1-unit change in x is associated with a 1.382 increase in the group's average value of y.

A logical next question is whether the group effect of x is different from the individual effect of x. I'm sure it's possible to test this after you run the model in R, but I am not as well-versed in R. So instead, I'll run a different form of the model in which I leave x uncentered and add in the group mean for x as a second predictor. It turns out that when you do so, the coefficient on group mean x is a test of whether the within and between effects of x are different from each other:

Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + group_mn_x + (1 | p)
   Data: data

REML criterion at convergence: 294

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.57439 -0.75814 -0.08622  0.68231  1.91484 

Random effects:
 Groups   Name        Variance Std.Dev.
 p        (Intercept)  5.715   2.391   
 Residual             21.468   4.633   
Number of obs: 50, groups:  p, 5

Fixed effects:
            Estimate Std. Error t value
(Intercept)    1.242      2.537   0.490
x             -1.683      0.688  -2.446
group_mn_x     3.065      1.030   2.975

Correlation of Fixed Effects:
           (Intr) x     
x           0.000       
group_mn_x -0.647 -0.668

The group_mn_x coefficeint of 3.065 has a standard error of 1.3, suggesting that it is significantly different than 0 at a p-value of .05 or smaller. Just to test it against the model without the group_mn_x:

> anova(mixed.model, mixed.model.grpmn)
refitting model(s) with ML (instead of REML)
Data: data
Models:
mixed.model: y ~ x + (1 | p)
mixed.model.grpmn: y ~ x + group_mn_x + (1 | p)
                  Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
mixed.model        4 314.33 321.97 -153.16   306.33                            
mixed.model.grpmn  5 308.15 317.71 -149.07   298.15 8.1787      1   0.004238 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

This model is favored both by AIC, BIC, and the Chisq is significant at p = .0042. Therefore, you have significant differences in the within- and between-effects of x on y. Economists would argue that this is a reason to abandon a mixed model and go with a fixed effects model because of their concern about endogeneity, but that is really up to the user. Personally, I like how the mixed modeling framework allows you to investigate these types of differences.

Sorry this was so long, but hopefully it shows why you get different coefficients when you estimate the regression of y on x using different methods.