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?


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):

enter image description here

From this code:


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

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


data$p = factor(data$p)

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

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

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:
x -0.029

Fixed Model

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

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

  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

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 
  • $\begingroup$ +1 and I don't see why you say that this definition is "refuted" by the linked answer. $\endgroup$
    – amoeba
    Dec 19 '17 at 14:27
  • $\begingroup$ @amoeba -- thank you, this question came about by my trying to understand a number of sources that seemed to me to be trying to move away from that definition. In the linked answer, for example, it says "Unfortunately, users of mixed effect models often have false preconceptions about what random effects are and how they differ from fixed effects. People hear "random" and think ... fixed effects have to be used when something is "fixed" while random effects have to be used when something is "randomly sampled". But there's nothing particularly random about..." $\endgroup$
    – justme
    Dec 19 '17 at 15:03
  • $\begingroup$ OK, but in the next paragraph Paul (the author of that answer) says that he thinks Gelman's definition #4 is the most accurate, and that one reads "If an effect is assumed to be a realized value of a random variable, it is called a random effect." I think this is largely semantics; mathematically the difference is very clear. $\endgroup$
    – amoeba
    Dec 19 '17 at 16:03
  • $\begingroup$ Admittedly there is a gap between the high-level conceptual understanding I provide in that answer and the incredibly popular implementation provided in Doug Bates's lme4 package. I hope people are not relying on that answer to perform modeling work and then blindly quote standard errors, which are clearly being computed by lme4 with more assumptions about what the user wants than my answer makes. $\endgroup$
    – Paul
    Dec 19 '17 at 17:32
  • 1
    $\begingroup$ @amoeba all I meant is that my answer pushes the viewpoint that there's nothing really all that "random" about random effects, beyond that they follow a distribution. But the lme4 standard errors show that some important people think of other aspects of randomness as being implicit in mixed models, like the groups being randomly sampled from some larger population. From my perspective, using a random effect doesn't mean you should be forced to quote a standard error that includes the randomness of sampling groups from a larger population, but I don't know enough to prove that right now. $\endgroup$
    – Paul
    Dec 19 '17 at 22:37

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.

  • $\begingroup$ apologies, stupidly I missed when I re-ran the complete output for you that I had commented out my "contrasts" line... I've updated it. I'm not sure if it changes much in your answer apart from the statement "the (intercept) is 12" -- which is now 21 as it was meant to be... I'm still trying to wrap my head around your answer, but just wanted to update to the correct output. $\endgroup$
    – justme
    Dec 19 '17 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.