Random slope random intercept model for two grouping levels

I'm trying to see how boldness of individuals change with time (plasticity in personality). My data consists of boldness scores (response variable), Trials (across time), individuals nested within colonies. We know that colonies differ from each other in their average boldness i,e some colonies are bolder than others because they consist of greater proportion of bold individuals. Thus, an individual's colony of origin has to be factored in, as individuals are not independent. My data consists of around 240 individuals from 14 colonies.

Since I want to see how boldness of individuals change with time, I would like to make random intercept random slope plots. However, model with (Trial|colony) + (Trial|individualID) is a singular fit, so I went with the following model:

M1 <- lmer(Boldness~Trial+(1|colony)+(Trial|individualID),data=mydata).


Please note that (1|colony)+(1|individualID) is the same as (1|colony/individualID) in this case (because of this excellent answer here). I got the following plot using the sjPlot package:

Fig 1:

Next, I created a model without colony as random intercepts:

M2<-lmer(Boldness~Trial+(Trial|individualID),data=mydata). I got the following plot:


Figure 2:

As one can see, Fig 2 shows greater variation in intercepts between individuals while this variation is reduced after factoring in colony as random intercept (Fig 1). Also Fig 1 shows fanning out, i.e individuals with greater intercepts have larger slopes.

My question is:

Is M1 correct if one is interested to see how slopes of individuals vary with time? As evident from Fig 2, there exists individuals with low boldness scores on the intercept. These individuals appear to get a higher boldness scores on the intercept when colonies are included as random intercepts (Fig 1). Thus, individuals with low boldness scores belonging to colonies which have high average boldness scores (and vice-versa) will have greater negative slopes in Fig 1, while in reality, that individual's low boldness scores may show very little variability across time/trials. Because the individual's intercept will be "pulled up" to match the higher average boldness at the level of the corresponding colony, and it's slope will go down over trials (which may explain fanning out).

I would like to know if my assumption is correct. If it is correct, then M1 is not the best way to model within-individual variation over time.

Many thanks!

• Here are two plots of the lines getting from the fitted models. How about raw data? You can draw a plot by observed data. Then check which plot from the model is similar to observed data. – user158565 Dec 17 '18 at 22:57
• I plotted the raw data through time and connected the dots for each individual. As I have many individuals, the plots look messy and uninformative. – BP86 Dec 18 '18 at 0:41
• From raw data, Did you find that variance of boldness increase with trials, or no big change along the trials? – user158565 Dec 18 '18 at 0:56
• Hard to say anything from the raw data plots. I have around 20 individuals in each colony. Lines go zig-zag. – BP86 Dec 18 '18 at 1:42
• How about not connect the points, just scatterplot? it will be easy to see the variation along the trials. – user158565 Dec 18 '18 at 2:17