# The correct random slope model for nested data

I'm trying to see how personalities of individuals change with time. The variables in my data are: 1. latency to emerge (response variable in continuous scale) measured for 204 individuals from 14 colonies. Each individual was measured for 6 times (6 trials). The 204 individuals came from 14 different colonies, so individuals (ID) are nested within colonies. Individuals have distinct IDs in the data set such that each ID can belong to only one of the 14 colonies (for example: IDs 1-20 in colony 1, 21-40 in colony 2 and so on). In this scenario, there is an excellent link here which shows that the model output will be the same for (1|colony)+(1|ID) and (1|colony/ID).

I would like to see how the latency to emerge of individuals (IDs) change with time/trials. Remember, individuals are not independent, but nested within colonies. Therefore, I made two models:

m1 <- lmer(Latency~Trial+(1|colony)+(Trial|ID),data=mydata) and

m2 <- lmer(Latency~Trial+(Trial|ID),data=mydata).


I know that model m2 assumes that individuals are independent of colonies (while in reality, individuals are non-independent of their respective colonies). Does model m1 account for this lack of independence? How is model 1 different from model 2 in terms of the random slopes plotted against each individual?

The model you want is probably

Latency~Trial+(Trial|colony)+(Trial|ID)


(or equivalently Latency ~ Trial + (Trial|colony/ID)): the corresponding statistical model for observation $$k$$ of individual $$j$$ in colony $$i$$ (at time $$t$$) is

$$\begin{split} y_ijk(t) = & (\beta_0 + \epsilon_{c0,i} + \epsilon_{d0,j}) + \\ & (\beta_1 t + \epsilon_{c1,i} + \epsilon_{d1,j}) t + \epsilon_{r,k} \end{split}$$

where $$\{\epsilon_{c0},\epsilon_{c1}\}$$ (colony-level random effect on intercept and slope) and $$\{\epsilon_{d0},\epsilon_{d1}\}$$ (individual-level random effect on intercept and slope) are bivariate Normal deviates, and $$\epsilon_r$$ (the residual error) is univariate Normal.

Your first model doesn't allow for among-colony variation in slopes; as you state in your question, your second model allows for neither intercept nor slope variation among colonies.

Depending on what your data looks like and what's typical in your field, you might want to consider a log-Normal model (use log(Latency) as your response variable) or a Gamma GLM with a log link ...

• Thanks much for the answer @Ben Bolker. Trial|colony+Trial|ID would generate random slopes for both colonies and IDs. I would like random slopes for only ID after controlling for random intercepts of colonies. Also, the models becomes a singular fit for 1|colony+Trial|ID. Is there anyway to deal with this. I know MCMCglmm can probably deal with this, but since I would like to have a random slope random intercept plot, I guess I have to stick to lmer? Any leads into this will be much appreciated. Thanks! – BP86 Dec 14 '18 at 17:21
• I'm curious about your reason for not wanting random slopes by colonies - it seems to be part of the maximal model, so should in principle be left in if there's enough data to fit it. Dealing with singular fits: you could use blme::blmer to regularize ... also, MCMCglmm (or brms or rstanarm) should I think be able to handle any of these models (including one with random intercepts+slopes at one level and random intercepts only at another level) ... – Ben Bolker Dec 14 '18 at 21:20
• Thanks much Ben Bolker. The question you raised is interesting. We know that colonies differ consistently from one another in their average behavioural response over time (which is depicted by the y-intercepts; also know as personality). Hence I would like (1|colony) in the model. I'm mainly interested how individuals vary in personality over time (which is depicted by slopes; also known as behavioural plasticity). – BP86 Dec 15 '18 at 12:54
• So a random intercept random slope plot with (1|colony)+(Trial|ID) will show the plasticity in behaviour of individuals over time, while normalizing for the personality (y-intercept) attributed at the level of colonies. I would also show a random intercept random slope plot at the level of colonies (Trial|ID). Ideally, it would be nice to subset the data into colonies and show (Trial|ID) for each colony. However, I encounter this problem with singular fits. – BP86 Dec 15 '18 at 12:54
• I tried MCMcglmm followed by plot_model (from sjPlot package) to generate the random intercept random slope plot, but instead of lines for each ID or colony, I get a figure that looks something like the one depicting a distribution. – BP86 Dec 15 '18 at 12:54

Let analyze the random effects in the model.

From the description, the response variable from different colonies are independent. Let check the 2 individuals from the same colony and measurements at time 1 and 6. $$Y_{11}, Y_{16}, Y_{21}, Y_{26}$$.

The model m2 is: $$Y_{ij} = X\beta + \gamma_{1i} + j\gamma_{2i} + \epsilon_{ij}$$ where $$Var(\gamma_{1i}) = \sigma_1^2$$, $$Var(\gamma_{2i}) = \sigma_2^2$$, $$Var(\epsilon_{ij}) = \sigma_0^2$$, Then the variance-covariance matrix of four $$Y$$s is $$\Sigma= \left(\begin{matrix}\sigma_1^2+\sigma_2^2+\sigma_0^2&\sigma_1^2+6\sigma_2^2&0&0\\ \sigma_1^2+6\sigma_2^2&\sigma_1^2+36\sigma_2^2+\sigma_0^2&0&0 \\ 0& 0 & \sigma_1^2+\sigma_2^2+\sigma_0^2&\sigma_1^2+6\sigma_2^2\\ 0&0&\sigma_1^2+6\sigma_2^2&\sigma_1^2+36\sigma_2^2+\sigma_0^2\\ \end{matrix} \right)$$

In this model, there are two problems: (1) individual 1 and 2 come from the same colony but covariance is zero, (2) $$Y$$ at the time 1 and 6 have different variance.

The model m1 is:

$$Y_{ij} = X\beta + \gamma_{1i} + j\gamma_{2i} + \gamma_3 +\epsilon_{ij}$$ Let $$Var(\gamma_3) = \sigma_3^2$$.

Then the ariance-covariance matrix of four $$Y$$s under m1 is adding $$\sigma_3^2$$ to all of 16 elements of $$\Sigma$$. So the first problem in m2 is resolved. But the second problem is still there. If you believe that variance of response variable increases along the time, m1 is model that you wanted. If you think it is unreasonable that variance of response variable increases along the time, then the model m3 maybe more reasonable.

$$Y_{ij} = X\beta + \gamma_{1i} + \gamma_3 +\epsilon_{ij}$$ Under model m3, the ariance-covariance matrix of four $$Y$$s is: $$\Sigma_3= \left(\begin{matrix}\sigma_1^2+\sigma_3^2+\sigma_0^2&\sigma_1^2+\sigma_3^2&0&0\\ \sigma_1^2+\sigma_3^2&\sigma_1^2+\sigma_3^2+\sigma_0^2&0&0 \\ 0& 0 & \sigma_1^2+\sigma_3^2+\sigma_0^2&\sigma_1^2+\sigma_3^2\\ 0&0&\sigma_1^2+\sigma_3^2&\sigma_1^2+\sigma_3^2+\sigma_0^2\\ \end{matrix} \right)$$
So when you construct the random effect part in the model, need to check which one in more suitable for your situation.

In R, m3 canbe m3<-lmer(Latency~Trial+(1|colony)+(1|ID),data=mydata)

The discussion above is based on the assumption that value of the time (variable trial) are 1,2,3,4,5 and 6 and is continue variable in the model.