# How are random effects being included in the linear mixed model along with the concept behind their calculation

I'm having difficulty in understanding the "process" that is going on behind how we are calculating all of our parameter estimates and how the random effects are used in our models.

To begin we can express the linear mixed model as:

$$y = X\beta + Zb + \epsilon \\b \sim N(0, \psi_\theta),\ \epsilon \sim N(0,\Lambda_\theta)$$

Where $$X\beta$$ would represent the fixed effects and $$Zb$$ would be representative of the random effects. How are both sets of effects being estimated? What I mean by this is I get that we will be using maximum likelihood methods to estimate the parameters formally. But what is the process? Are we estimating the fixed effects $$X\beta$$ and the random effects $$Zb$$ separately and then bring them together in our model?

I was playing around in R to try and understand more, but I'm still stick in making the leap. Here is an output I got from fitting a LMM to the iris data set:

Linear mixed model fit by REML ['lmerMod']
Formula: Petal.Width ~ Sepal.Width + Sepal.Length + (1 + Sepal.Length |      Species)
Data: iris

REML criterion at convergence: -66.7

Scaled residuals:
Min       1Q   Median       3Q      Max
-2.76049 -0.54295 -0.08282  0.55066  2.74867

Random effects:
Groups   Name         Variance Std.Dev. Corr
Species  (Intercept)  0.24352  0.4935
Sepal.Length 0.01377  0.1173   0.40
Residual              0.03091  0.1758
Number of obs: 150, groups:  Species, 3

Fixed effects:
Estimate Std. Error t value
(Intercept)   0.20829    0.33520   0.621
Sepal.Width   0.27272    0.05248   5.196
Sepal.Length  0.01796    0.07726   0.232

Correlation of Fixed Effects:
(Intr) Spl.Wd
Sepal.Width -0.116
Sepal.Lngth  0.134 -0.279
> coef(ir_lme2)
$Species (Intercept) Sepal.Width Sepal.Length setosa -0.1400153 0.2727164 -0.10956019 versicolor 0.0573940 0.2727164 0.08658932 virginica 0.7074925 0.2727164 0.07684240 attr(,"class") [1] "coef.mer  So I have fixed effects and random effects as can be seen at the bottom of the outputs. With regards to the random effects I get that we "group" our observations and then use those groupings to get group estimates of our parameters which are the random effects below. Were those random effects estimated in isolation and if so how? Same question with respect to the fixed effects. Another question is how are those random effects playing into the fixed effect estimates? Are the random effects contributing to the values we see in the fixed effects read outs? I read previous articles on the site about the ideas: What is the difference between fixed effect, random effect and mixed effect models? What is a difference between random effects-, fixed effects- and marginal model? I had also asked a previous question about this, but I might delete it because it is muddied in my confused understanding of the concept. As you can see my head is all over and I'm very confused about how things are being put together in this model. Any help to clarify things would be appreciated. Even in chat because I feel the things I'm not getting should be easy to clear up. EDIT: So I attempted to get some clarification from a TA. In my example let's say we end up with the expression: $$y = \beta_0 +\beta_{SW} \cdot SW + \beta_{SL} \cdot SL + (\alpha_1 + \alpha_2 \beta_{SL}SL)$$ where $$\beta_i$$ correspond to the fixed effects and $$\alpha$$ corresponds to the random effects, $$SW =$$ sepal.width, $$SL =$$ sepal.length, and $$\alpha =$$ random effect from Species. So if I understand this correctly for the random effect of random slope, we would group by species, take all of the Sepal.Length values by species (let's use setosa as a concrete example), compute an estimate for the variance, use this estimate for the variance in a normal distribution $$\alpha \sim N(0, \sigma_{setosa}^2)$$ from which we would draw a random value for Sepal.Length, and then this would serve as the random factor $$\alpha_2$$ which we would multiply by $$\beta_{SL}$$ to get our value for the random slope? Not looking for the precise mathematics yet, just an understanding. • Is your question about the process of estimation or about how to intuitively view the process of data generation? Mar 15, 2022 at 13:38 • would something like the example at the end of this question/answer help? I could remake it based on the formulas that you have with the iris flowers. Mar 15, 2022 at 18:38 • That example has only random intercepts. The slope is the same. So those "lines" are made by a (random) shift in vertical direction for points that belong to the same group. Mar 15, 2022 at 19:14 • An example of random slopes is shown here note the plot with the 35 different plots. Mar 15, 2022 at 19:16 • Do you get the expression$Zb$if it would be a fixed effect? One source of confusion might be that$Z$is a categorical variable, and so possibly it is difficult to imagine that. Mar 15, 2022 at 19:19 ## 2 Answers You can think of mixed models as a two stage modeling approach. Firstly, you fit a model irrespective of the random effects; secondly you model the effect for each level of the grouping factors (random effects) via an approach known as partial pooling, see here and here for some more explanations and details. Finally, you adjust the fixed effects model based on the random effects. All of this happens together when running a mixed-effects model. Here is an example using the sleepstudy data in R: > m <- lmer(Reaction ~ Days + (1|Subject), data = sleepstudy) > fixef(m) (Intercept) Days 251.40510 10.46729  Which is the same as this: > m2 <- lm(Reaction ~ Days, data = sleepstudy) > coef(m2) (Intercept) Days 251.40510 10.46729  Going back to the lmer model, the random slopes estimates are: > ranef(m)$Subject
(Intercept)
308   40.783710
309  -77.849554
310  -63.108567
330    4.406442
331   10.216189
332    8.221238
333   16.500494
334   -2.996981
335  -45.282127
337   72.182686
349  -21.196249
350   14.111363
351   -7.862221
352   36.378425
369    7.036381
370   -6.362703
371   -3.294273
372   18.115747

with conditional variances for “Subject”


Now you adjust the average (fixed effect) intercept for each subject based on the estimated random effect, let's look at Subject 308 as an example:

$$251.40510 + 40.783710 = 292.1888$$

The result can also be checked by looking at coef(m) as well:

> coef(m)
$Subject (Intercept) Days 308 292.1888 10.46729 309 173.5556 10.46729 310 188.2965 10.46729 330 255.8115 10.46729 331 261.6213 10.46729 332 259.6263 10.46729 333 267.9056 10.46729 334 248.4081 10.46729 335 206.1230 10.46729 337 323.5878 10.46729 349 230.2089 10.46729 350 265.5165 10.46729 351 243.5429 10.46729 352 287.7835 10.46729 369 258.4415 10.46729 370 245.0424 10.46729 371 248.1108 10.46729 372 269.5209 10.46729  And now here's your example that also includes random slopes: > m3 <- lmer(Petal.Width ~ Sepal.Width + Sepal.Length + (1 + Sepal.Length |Species), data = iris) > fixef(m3) (Intercept) Sepal.Width Sepal.Length 0.20829042 0.27271644 0.01795717 > ranef(m3)$Species
(Intercept) Sepal.Length
setosa      -0.3483057  -0.12751737
versicolor  -0.1508964   0.06863214
virginica    0.4992021   0.05888522

with conditional variances for “Species”


Let's combine the fixed and random intercept and slopes together for setosa as an example:

Intercept:

$$0.20829042 - 0.3483057 = -0.1400153$$

Sepal.Width stays the same (no adjustment, i.e. not included in random effect terms):

$$0.27271644$$

Sepal.Length:

$$0.01795717 - 0.12751737 = -0.10956019$$

Let's check with:

> coef(m3)
\$Species
(Intercept) Sepal.Width Sepal.Length
setosa      -0.1400153   0.2727164  -0.10956019
versicolor   0.0573940   0.2727164   0.08658932
virginica    0.7074925   0.2727164   0.07684240


Here is a good explanation and an example that is not too complicated to follow and easy to understand: https://m-clark.github.io/mixed-models-with-R/random_intercepts.html#the-mixed-model

Another helpful link understanding complete pooling, no pooling and partial pooling: https://www.r-bloggers.com/2017/06/plotting-partial-pooling-in-mixed-effects-models/

• Aw man...I wish I gave you the bounty because your answer brings so much more clarity to the situation. One more question with regards to the random estimates. When we estimate the random effect (be it intercept, slopes), are we calculating it by "restricting" our data to the specific grouping and then calculating the parameter? For example in my model I restricted by species, so we restrict our data to the setosa species and perform the necessary calculations to get a slope and intercept and then those are the random effect slope and intercept we added to the fixed. Mar 23, 2022 at 0:23
• @dc3rd Have a look at the links above on complete pooling, no pooling and partial pooling. Suppose we have 10 subjects and 4 to 8 measurements for each subject. We can now calculate the mean of the whole dataset (complete pooling) which would give us a reliable estimate for the outcome but ignoring the contribution of each subject. Next, we could also calculate the mean for each subject and end up with 10 means (no pooling) which are likely not as reliable as well as susceptible to outliers since only 4 to 8 measurements went into the calculation [...] Mar 23, 2022 at 4:08
• [...] partial pooling is a "comprise" between complete pooling and no pooling in that it shrinks group observations with fewer observations or larger variances towards the complete pooling average. How exactly the group level estimates/adjustments for the random term components are calculated, I am not sure. For me this is all I know about this process conceptually. Here's is another link that visualizes this tjmahr.com/plotting-partial-pooling-in-mixed-effects-models Mar 23, 2022 at 4:08
• I read your mixed models link from Clark...helped a lot. Will check these ones too. Thanks for the help my fellow Canuck. Mar 23, 2022 at 6:45
• No problem @dc3rd :D I am glad you found this information useful. Mar 23, 2022 at 13:36

The long and short of the answer is: yes, the random and fixed effects are calculated/estimated at the same time. However, my sense is that you are not particularly interested in the mathematics behind that process. Instead, it seems that you are looking for a conceptual understanding of the estimation process. For this, it is OK to think of the processes happening separately.

And for this explanation, you would want to think of the data in groups. But first start with everything in one big group; here we would obtain the fixed effects. Now, inside of each group, we would obtain a solution for the remaining effects that weren't estimate in the big group (the random effects)...but, as we are just focusing on one group, these could also be seen as the "rest of the effects" (as opposed to the random effects). In other words, within each group estimation, there is nothing "random" happening with regards to the solution.

Now, once you've done this for all groups, each group will have the same fixed effects estimates, and each group will have distinct effects estimates for all the other effects. These are collectively the random effects. If you wanted to, you could then take any one of your random effects, compile the by-group parameters into one data set, and the mean and variance of this is (approximately) the mean and variance reported for your random effect in the output above.

Again, to clarify, this is not how the actual estimation process obtains the parameters for the random effects (namely the mean and variance of the parameter), but it conceptually is what is "happening under the hood" of the engine generating these results. (The estimation process uses a more sophisticate approach so as to obtain even better estimates for the variance and associated standard errors...and I will leave the explanation of that to someone else.)

Hope you find this useful.

• I did. Thank you very much. I know there is more to it, but why dorwned myself in the mathematics if I have no idea of what is somewhat happening conceptually. Now once I read this explanation over again and the suggested other posts, then I can drown myself in the "fun" stuff. Cheers. Mar 22, 2022 at 1:44