Troubleshooting linear mixed effects model plot and trendline slopes

Question

I'm working with multiple years of temperature data from different sites in a watershed. I'm using the 7DADM measurements i.e. one temperature reading for each site per year of monitoring. This is what my data looks like

ggplot(use_sevdadm, aes(x = year, y = maxroll, colour = sitename)) +
  geom_point(alpha = 0.5) +
  theme_classic() +
 theme(legend.position = "none",
        panel.spacing = unit(2, "lines"))

I've written a mixed effects model using R to see if there is a trend for each site for the years of data I have (15 years, 2008-2022):

glmmTMB(data = df, formula = temp ~ year + (1+ year|sitename))->mod1

As I understand it, I'm saying for each temperature what is the interaction between the year grouped by each site. And the model that it produces should have a random intercept and random slope for each site.

This results in:

summary(mod1)

Formula:          maxroll ~ year + (1 + yr_step | sitename)
Data: df

     AIC      BIC   logLik deviance df.resid 
   776.8    797.7   -382.4    764.8      232 

Random effects:

Conditional model:
 Groups   Name        Variance  Std.Dev.  Corr  
 sitename (Intercept) 4.142e+00 2.0351384       
          year        7.289e-07 0.0008537 -0.97 
 Residual             1.096e+00 1.0468734       
Number of obs: 238, groups:  sitename, 17

Dispersion estimate for gaussian family (sigma^2):  1.1 

Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 16.08723    0.51156  31.447  < 2e-16 ***
year      0.05854    0.01601   3.657 0.000255 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

using ranef(mod1) I see:

              (Intercept)         year
S1          -2.619128562  1.061847e-03
S2           0.016190616 -6.782850e-06
S3           1.571495483 -6.369153e-04
S4           6.298859587 -2.553502e-03
S5          -1.338048169  5.424395e-04
S6           0.483008860 -1.969374e-04
S7          -2.688549545  1.092067e-03
S8          -0.881644350  3.573641e-04
S9           0.512919081 -2.081534e-04
S10         -0.973586630  3.951104e-04
S11          0.008698755 -3.135866e-06
S12          1.446483921 -5.874260e-04
S13          0.818152768 -3.314534e-04
s14          0.371879453 -1.516007e-04
S15         -2.392803620  9.706963e-04
S16         -0.235775404  9.517781e-05
S17         -0.398206974  1.612278e-04

To me this means that the slopes (under year) are different for each site which should be visible when I plot the model with the data however when I do plot the data using ggplot() I get this:

ggplot(use_sevdadm, aes(x = year, y = maxroll, colour = sitename)) +
  geom_point(alpha = 0.5) +
  theme_classic() +
  geom_line(data = cbind(df, pred = predict(mod1)), aes(y = pred), size = 1) +  # adding predicted line from mixed model 
  theme(legend.position = "none",
        panel.spacing = unit(2, "lines"))

See how the slopes all look the same? My question is why are the lines not plotting the random slopes of each site that the model is giving me? Is there a way to figure out the problem or force it to plot the lines with the ranef() I have? Some thoughts I've had that could be why this is happening: the slopes are so small there's an imperceptible change amongst the lines; I wrote the formula wrong in the model;I should use the coef() function instead; I need more data points? I'm using 7 day maximums for each site/ year that I've extrapolated from continuous summer daily maximums so there is more data I can use but it just gets unwieldy quickly.

I'm envisioning something that would look like this figure on the right (B) but with my data:

You can find the article this figure came from here

Thank you for any and all help!

Answer 1

It looks like there just isn't much variability in the slopes between different sites. You can see this in two ways:

• The random slopes that you see under year when you run ranef(mod1) are all extremely small. These random effects represent how much (or little, in your case) each site's slope deviates from the global (fixed) slope (0.05854)
• Take the variance of the random intercepts (4.142e+00 = 4.142), the variance of the random slopes (7.289e-07), and the residual variance (1.096e+00 = 1.096). The variance of the random slopes only accounts for $\frac{7.289e-07}{4.142 + 7.289e-07 + 1.096}=1.391561e-07\approx0\%$ of the total variance. Meanwhile, the variance of the random intercepts accounts for $\approx79\%$ of the total variability

Therefore, you should not expect to see different slopes for the different sites. You could opt for a more parsimonious model and drop the random slopes i.e.

glmmTMB(data = df, formula = temp ~ year + (1|sitename))->mod1