I am running a linear mixed effects model for time series data using R-INLA. My response variable is normally distributed. The model has a random intercept, and temporal autocorrelation between data points is modelled with a first order autoregressive structure. The parameter estimates from the model match our hypotheses. However, when I check diagnostic plots for the model, the residuals vs fitted values plot has a very clear pattern. Other diagnostic plots (like the QQ plot, etc.) look fine. I’m assuming that this residual vs fitted plot is telling me that there is a non-linear relationship between my response variable and my fixed effects variables. But I am not sure how to fix this. I was wondering if you could please give me some suggestions.

Here is the plot: enter image description here

Here is the code for the model:

m_inla1 <- inla(y~ x1 +
                f(group, model = "ar1"),
                data = dat,
                control.compute = list(waic = TRUE, dic = TRUE, cpo = TRUE, return.marginals.predictor=TRUE))

Where y is a continuous variable. Four of the covariates (x1, x2, x3, and x4) are also continuous variables. These have been standardised using the r function ‘scale’, which standardises each value by subtracting the mean and dividing by the standard deviation. The final covariate (x5) is a binary categorical variable that takes the value 0 or 1. The variable ‘group’ is the random intercept and has 48 levels.

Thanks so much!

  • $\begingroup$ Welcome to CV! It would really help if you provided a lot more information about how you coded your model as well as some descriptions of your variables (which may give clues on what should be done with this mixed model). $\endgroup$ Commented Apr 16, 2023 at 14:05
  • $\begingroup$ Thanks Shawn, I've added some more info about the model and the data :) $\endgroup$ Commented Apr 17, 2023 at 6:17

1 Answer 1


This pattern of data points that seem to be aligned on a straight line resembles the pattern that might occur with the Yule-Simpson effect.

The effect within the groups is different from the main effect.

example of Simpson Yule effect

### create data
m = 10
n = 10

z = rnorm(m)
z = rep(z,each=n)
x = rnorm(m*n,z)
y = 3*z-x

### plot data
plot(x,y, main = "x versus y")

### add linear fit
mod = lm(y~x)

### plot residuals 
plot(predict(mod),y-predict(mod), main = "prediction versus residual")
  • $\begingroup$ By Simpson-Yule, you are referring to this effect, right? If so, it may be helpful to include that in the answer. $\endgroup$ Commented Apr 17, 2023 at 7:55
  • 1
    $\begingroup$ @ShawnHemelstrand yes I was referring to that effect. I didn't had time to search for an appropriate link because my phone was running out of battery. $\endgroup$ Commented Apr 17, 2023 at 8:52
  • $\begingroup$ Thanks both :) I've included 'group' as a random intercept in the model. Would including it as a random slope alleviate the problem? Or is there another way you could suggest to deal with the Simpson-Yule effect? $\endgroup$ Commented Apr 17, 2023 at 10:42
  • $\begingroup$ +1 thanks for clarifying $\endgroup$ Commented Apr 17, 2023 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.