# Mixed effects models - clarification on random effects

I'm hoping to get some clarification on a "mismatch" between some simulated data I'm creating and the resulting model fit by R's nlme library. Specifically the random effects parameter estimates. I'm desperately trying to balance the length of the query with enough details to be complete and helpful to others who might stumble upon a similar situation...!

The experimental setup is as follows:

• Two groups of animals (fasted or fed).
• Multiple animals in each group.
• Multiple images (microscopy) taken within each animal. (Nesting!)
• Each of those images measures the intensity of a biological marker stain as a function of normalized distance ($$x$$) from a feature of interest. For ease, say 10 measurements are taken for each image and assume the relationship between the stain intensity and distance is a straight line.

So the primary objective is to assess the "treatment" (fed/fasted) effect - does it change that "stain profile". I'm happy to say that the simulation and fitted model are in agreement, as expected. However, I'm not sure if the lack of agreement when it comes to random effects is expected, or due to something I'm missing.

To simulate, I start with the fixed treatment effects. Namely, the "fixed portion" of the response variable (stain intensity) $$y_{ijk}$$ (for animal $$i$$, image $$j$$, and the k=1,...,10 sampled distances) is given by

$$y^{fixed}_{ijk} = \beta_0 + \beta_1I(diet=fed) + \beta_2x_{ijk} + \beta_3I(diet=fed)x_{ijk}$$

where $$I(a=b)$$ is an indicator function for the diet (returns 1 if a=b, else 0) and $$x_{ijk}$$ are the normalized distances for measurement $$k$$ in image $$j$$ of animal $$i$$. Since I'm simulating the data, I set those $$\beta$$ values.

On top of those fixed effects, I add random intercepts and slopes for each animal, both randomly sampled from normal distributions. Concretely, for each animal $$i$$, I add intercept $$\alpha^{(0)}_i \sim \mathcal{N}(0,\eta_0^2)$$ and slope $$\alpha^{(1)}_i \sim \mathcal{N}(0,\eta_1^2)$$. Given that these random draws are performed independently, I take this to mean I have no correlation between the intercept and slope, i.e. $$\mathrm{cov}(\alpha^{(0)}_i, \alpha^{(1)}_{i}) = 0$$. The variances $$\eta_0^2$$, $$\eta_1^2$$ are known.

After generating the random effects at the animal level, I turned to the "image" level; recall images are nested within animals by design. Similarly, I generate image-specific random intercepts ($$\gamma^{(0)}_{ij} \sim \mathcal{N}(0,\omega_0^2)$$) and slopes ($$\gamma^{(1)}_{ij} \sim \mathcal{N}(0,\omega_1^2)$$) for image $$j$$ in mouse $$i$$. As above, these draws are performed independently and we know $$\omega_0^2$$, $$\omega_1^2$$. Additionally, there is no correlation between the animal and image levels, i.e. $$\mathrm{cov}(\alpha^{(0)}_i, \gamma^{(0)}_{ij}) = 0$$

Finally, I add these all up:

$$y_{ijk} = y^{fixed}_{ijk} + \alpha_{i}^{(0)} + \gamma^{(0)}_{ij}+ ( \alpha_{i}^{(1)} + \gamma^{(1)}_{ij})x_{ijk} + \epsilon_{ijk}$$ where $$\epsilon_{ijk}$$ is iid sampled from $$\mathcal{N}(0,\sigma^2)$$. Since the $$\alpha$$ and $$\gamma$$ were normal iid, we can combine those terms, e.g. $$\delta_{ij}^{(q)}=\alpha_{i}^{(q)} + \gamma^{(q)}_{ij} \sim \mathcal{N}(0, \eta^2_q + \omega^2_q)$$ for $$q=0,1$$. Based on the known variance parameters I'm using for the simulation, I'm expecting the estimated random intercept/slope parameters to be relatively close to these $$\delta$$'s

Currently, I've tried a couple models, all of which recover the expected fixed effect $$\beta$$ parameters (great!). However, none of the random effect parameter estimates come close to the expected structure I have "baked into" the simulation. This is where I'm a bit stuck. Or is it too much to hope that repeated use of rnorm will produce truly uncorrelated random intercepts/slopes such that the fitted model will be "close" to my known variance parameters?

The first model (which I know assumes a general covariance structure) is:

> model.1 <- lme(y ~ diet*x, random= ~ x | animal/image, data=mock_data)


To put concrete numbers, I used:

N_animals <- 50
N_images_per_animal <- 10
sampled_points <- 10

animal_intercept_variance <- 0.015^2 # \eta0
animal_slope_variance <- 0.02^2 # \eta1

image_intercept_variance <- 0.06^2 # \omega0
image_slope_variance <- 0.08^2 # \omega1

sigma = 0.085


and the model reported (in part)

> summary(model.1)
Linear mixed-effects model fit by REML
...
Random effects:
Formula: ~x | animal
Structure: General positive-definite, Log-Cholesky parametrization
StdDev      Corr
(Intercept) 0.002882329 (Intr)
x           0.004880377 -0.166

Formula: ~x | image %in% animal
Structure: General positive-definite, Log-Cholesky parametrization
StdDev     Corr
(Intercept) 0.05612260 (Intr)
x           0.07762704 -0.379
Residual    0.10186092
...


(the estimates are quite variable which probably is a clue that my simulation is not matching the theoretical model with perfectly uncorrelated random variables)

Alternatively, I tried a model where I dictate the lack of correlation per my simulation's construction. Since random = ~ x || animal/image (note the double-bar) is not accepted by nlme, I encoded the image factor so that the names were unique (e.g. "animal_1_image_3") per Crossed vs nested random effects: how do they differ and how are they specified correctly in lme4?. I then created the following model:

> model.2 <- lme(y ~ diet*x, random= list(animal = pdDiag(~x), image = pdDiag(~x)), data=mock_data)
> summary(model.2)
Linear mixed-effects model fit by REML
...
Random effects:
Formula: ~x | animal
Structure: Diagonal
(Intercept)           x
StdDev:  0.01666036 0.008468225

Formula: ~x | image %in% animal
Structure: Diagonal
(Intercept)          x   Residual
StdDev:  0.04158953 0.05053244 0.09226515
...


this also does not produce estimates close to my known params. However, I'm also uncertain what the full implication of that model is. As far as I understand, this will assume no correlation between the intercept and slopes at the animal and image levels ($$\mathrm{cov}(\alpha^{(0)}_i, \alpha^{(1)}_{ij}) = 0$$, $$\mathrm{cov}(\gamma^{(0)}_i, \gamma^{(1)}_{ij}) = 0$$), but will it also assume no correlation between the two? e.g. that $$\mathrm{cov}(\alpha^{(0)}_i, \gamma^{(0)}_{ij}) = 0$$?

So the bottom-line is either a) is there a better way to model the simulation I'm describing or b) is my expectation of "closely" capturing the expected random effects params unrealistic in practice. Or both!? Thanks!

• At first sight, it looks like the stddev.'s for images are okay, but not the ones for animals. Try more animals, 500 e.g.
– BenP
Commented Mar 19 at 10:44
• that does seem to inch it closer (and naturally decrease the variability of the estimates). Given how close the fixed params are to their true values, I was maybe expecting closer estimates for the random params. Unfortunately in practice, these experiments only ever have <10 animals... Commented Mar 19 at 12:43
• For the small number of animals 50, you could also run the simulation say 1000 times and calculate the mean of the 1000 estimates of the random effects. These should be close to the true values. I do not see anything wrong in your model description.
– BenP
Commented Mar 19 at 13:14

From the discussion in comments, I think that this is due to the high variance of sampling variance, even from a normal distribution. Remember that a "random effect" value is an estimate of the underlying variance around the fixed-effect value, based on the assumption of a normal distribution.

This is quite clear in sampling from a single normal distribution. Given a fixed number of observations, the sampling variance of a mean-value estimate (analogous to your fixed-effect estimates) is proportional to the variance in the population, but the sampling variance of the variance is proportional to the square of the variance in the population. See this page. With only 10 animals, you can't get very reliable estimates of underlying variances in practice.

• This is interesting and convincing (+1). Although, in the random effects model there are no directly observable random effects and hence the estimation of there variance is different from the situation in the paper, if I'm correct. But the idea that it is difficult to obtain a precise estimate of the variance with few "random animals" certainly makes sense. That is why it surprises me that often random effects (instead of fixed) are assumed and their variances estimated in situations with less than say 10 "animals". Isn't that strange???
– BenP
Commented Mar 19 at 19:40
• @BenP with small numbers of animals, a mixed model is presumably used to account for correlations within animals or within correlated groups like litters. I don't think much is (or at least should be) made of the estimates of random-effect variances; the important thing is trying to correct the covariance estimates for the fixed-effect coefficients.
– EdM
Commented Mar 19 at 20:31
• That's true, yes. In longitudinal data with say 3 fixed occasions it matters to model the covariances "correct" as far as possible. If one models compound symm. but the true covariances are AR1, that could lead to wrong p values for the fixed effects. Hence my doubts about a 'bad" estimate of the variance based on a small number of animals.
– BenP
Commented Mar 20 at 8:23