Random vs fixed effects: why is the standard error the same for the slope but wildly different for the intercept?

Question

I've noticed that when I compare my results from fixed and mixed effect models, I get very similar standard errors for the slopes of (fixed-effect) covariates, whereas on the intercepts, the fixed effects model can have far far smaller standard errors (while the parameter estimates themselves are almost identical). Why is this?

---

EXAMPLE

For example I get these results (abridged here; full results at the bottom):

Mixed Model
              Estimate Std. Error t value
 (Intercept) 21.040772   1.921112   10.95
 x           -1.007014   0.009933 -101.38

Fixed Effects Model
              Estimate Std. Error t value Pr(>|t|)    
 (Intercept) 21.043157   0.055888  376.52   <2e-16 ***
 x           -1.007443   0.009934 -101.41   <2e-16 ***

For data that looks like this (different symbols for different levels of p):

From this code:

set.seed(1)    

generatepoint = function(ppt,n){
  x = rnorm(n)+ppt
  y = -x + 10 + 2*ppt + rnorm(n,sd=.1)
  cbind(x=x,y=y,p=ppt)
}

data = as.data.frame(do.call(rbind,lapply(1:10,generatepoint,10)))

plot(y~x,data=data,pch=data$p)

data$p = factor(data$p)

cat("\n\nMixed Model\n\n")
mod = lmer(y~x+(1|p),data=data)
print(summary(mod))

cat("\n\nFixed Model\n\n")
contrasts(data$p) <- contr.sum(10)
fmod = lm(y~x+p,data=data)
print(summary(fmod))

Here's the full output:

Mixed Model

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

REML criterion at convergence: -88.6

Scaled residuals: 
  Min      1Q  Median      3Q     Max 
-2.6303 -0.6669 -0.1463  0.6713  2.2773 

Random effects:
  Groups   Name        Variance Std.Dev.
p        (Intercept) 36.87548 6.07252 
Residual              0.00807 0.08983 
Number of obs: 100, groups:  p, 10

Fixed effects:
  Estimate Std. Error t value
(Intercept) 21.040772   1.921112   10.95
x           -1.007014   0.009933 -101.38

Correlation of Fixed Effects:
  (Intr)
x -0.029


Fixed Model


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

Residuals:
  Min       1Q   Median       3Q      Max 
-0.23547 -0.06006 -0.01333  0.06017  0.20400 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
(Intercept) 21.043157   0.055888  376.52   <2e-16 ***
  x           -1.007443   0.009934 -101.41   <2e-16 ***
  p1          -9.009845   0.051524 -174.87   <2e-16 ***
  p2          -7.017192   0.045469 -154.33   <2e-16 ***
  p3          -5.005483   0.036105 -138.64   <2e-16 ***
  p4          -3.034799   0.029087 -104.34   <2e-16 ***
  p5          -0.993759   0.027279  -36.43   <2e-16 ***
  p6           0.992106   0.028082   35.33   <2e-16 ***
  p7           2.966062   0.030515   97.20   <2e-16 ***
  p8           4.987105   0.032854  151.79   <2e-16 ***
  p9           7.100627   0.045236  156.97   <2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08983 on 89 degrees of freedom
Multiple R-squared:  0.9993,    Adjusted R-squared:  0.9992 
F-statistic: 1.194e+04 on 10 and 89 DF,  p-value: < 2.2e-16

Answer 1

OK, I think I understand this now. I figured one way to understand it might be through a parametric bootstrap of mod using bootMer as an alternative way of generating the standard errors. This function works by randomly sampling new datasets according to the parameters fitted in a given mixed model from lmer, then fitting a new model to that simulated data, to see how the parameter estimates vary for different samples from the same population. It offers one argument in particular, use.u, which allows the user to choose whether the random effects are resampled for each sample, or not. If we choose not to resample them, then we are simulating a world where the random effects come out at the same values each time -- sounds more like fixed effects? -- well, choosing that setting gives us the same standard errors as from the fixed effects model! On the other hand, choosing to simulate u, we get the same standard errors as from the mixed-effects model. (See below for code and output!)

Of course, this fits with the idea of a random effect as one that will come out with different levels each time, as opposed to a fixed effect that does not change (albeit my only confusion in answering this question is that that seems to be a definition that is refuted in some highly rated StackExchange answers such as this one). But that then makes sense of the standard errors presented in the question: if fixed effects are indeed, "fixed", then we are making inference in the context of those specific levels of the fixed effects: in the example in the question, in a repeat of the experiment, we would expect the dots to jiggle around within the bounds of the lines they are currently in, but those lines themselves would not move much. On the other hand, if the intercepts of those lines are random effects, then the lines should be expected to appear in different places in each new sample, and our idea of where the "middle" is, is far less clear: for all we know, this sample might have come out with all ten levels of the random effect being below the mean, for example, in which the case the true mean could even be above all the lines in the plot. On the other hand, the lines themselves are consistently pointing at a gradient of -1 and whether we view the heights of those lines as fixed or random, we see the same behaviour within a given line, so we should expect to draw the same inference around the gradient in both models.

---

Code and output:

> booted.simulate.u = bootMer(mod,fixef,nsim=10000,use.u=FALSE)
> booted.fix.u      = bootMer(mod,fixef,nsim=10000,use.u=TRUE )

> print(apply(booted.simulate.u$t,2,sd))
(Intercept)           x 
1.915025294 0.009903413 

> print(apply(booted.fix.u$t,2,sd))
(Intercept)           x 
0.055405505 0.009868006 

Answer 2

The rationale for fitting a fixed or random intercept is much the same. The choice ultimately boils down to power: considerably more power is required to fit the fixed effect with one level for each cluster. You have 10 clusters with 90 observations.

The term called "intercept" in the fixed and mixed effects model has a different interpretation. For the fixed model the intercept is generated by contrast to estimate the predicted mean in the $p_1$ cluster with $x=0$. In the mixed model, the intercept term is the expected value for a person with $X=0$ drawn from no cluster in particular (or, similarly, averaging up over all the possible clusters to which that person may belong).

You'll note in the fixed effects output that $p_1$ is not included as a term. The (intercept) is 12. The subsequent $p_2, p_3, \ldots$ terms provide the contrast as a mean difference from $p_1$. So they tend to be 2 units higher per your simulation.

Your use of contrasts then marginalizes the fixed intercepts into a grand mean which is centered at 21 which you'll note is about the average of the sequence of generated means for $p$ of [12, 14, 16, 18, 20, 22, 24, 26, 28, 30]. The standard error for this marginalized coefficient is going to be more precise than the value you would get from ignoring clustering and merely looking at the intercept term in a fixed effect model which only adjusts for $x$.

So this value has the same interpretation as the mixed model. Why is the inference different? The Wald based inference (obtained from summary.merMod) is inappropriate in this case. The documentation in R points to using bootMer or confint with profile likelihood. The joint distribution of fixed and random effects is hierarchical, so as with the fixed effects model, you need to perform some kind of marginalization. The bootstrap and profile likelihood methods are akin to MCMC and numerical integration. That's why these results finally agree.