Solution (source code)

= Solution

Add a location-specific <random intercept> to the <random-slope linear mixed model>:
$$
\boxed{Y_i=\beta_0+\beta_1t_i+b_{j(i)}t_i+u_{\ell(i)}+\varepsilon_i,\qquad u_\ell\overset{\rm iid}{\sim}N(0,\omega^2).}
$$
Take the location effects, incubator slopes and individual errors to be mutually independent, with the original <variances> $\tau^2$ and $\sigma^2$. The <random intercept> models variation already present at hatching because of the origin of the eggs. Because larvae are subsequently assigned to incubators, location and incubator are <crossed random effects>, not automatically nested groups. The corresponding code is
``
fly.model.loc <- lmer(size ~ hours + (0+hours|incubator) + (1|location))
``
The resulting <covariance> between two observations contains $\omega^2$ if they share a location and $\tau^2t_it_k$ if they share an incubator. If the three named locations themselves are the only populations of interest, use a fixed location factor instead:
``
fly.model.loc.fixed <- lmer(size ~ hours + location + (0+hours|incubator))
``
Choosing between these interpretations depends on the inferential target; with only three locations, a location <variance> component will be weakly estimated.