Solved – How to deal with autocorrelation in mixed models

autocorrelationmixed modelrrandom-effects-modelresiduals

I am trying to model a varible (maximum depth) as a function of type of dive and diel changes (day,night) with the individuals (whales in this case) as a random factor in R.

I tried first to apply a linear mixed model (lme) and I had a problem of autocorrelation and non-normality of residuals.

Next, I tried to apply a GLM with Poisson and negative binomial distribution. Both had the same problem of autocorrelation and/or non-normality of residuals.

What can I do after to model these variables or to correct these issues?

Thank you

Update

I tried using corAR1 (correlation = corAR1(form = ~ 1 | whale)) but when I plot ACF there is still a strong autocorrelation.

What can this mean and how could I proceed to solve it?

I could not include time in the corAR1 function because of repeated diel values (Error in Initialize.corAR1(X[[i]], …) : covariate must have unique values within groups for "corAR1" objects)

Thanks!

Best Answer

Mixed models are often a good choice when you have repeated measures, such as here, within whales. lme from the nlme package can fit mixed models and also handle autocorrelation based on a AR(1) process, where values of $X$ at $t-1$ determine the values of $X$ at $t$.

$$ X_{t}=c+\varphi X_{{t-1}}+\varepsilon _{t} $$

The function corAR1 handles discrete time and corCAR1 handles continuous time. For example:

mod <- lme(..., correlation = corAR1(form = ~ 1 | id))`

with equally spaced time intervals or:

lme(..., correlation = corAR1(form = ~ time | id))

where the time intervals are supplied by the time variable.

I tried to apply a GLM with Poisson and negative binomial distribution. Both had the same problem of autocorrelation and/or non-normality of residuals.

From your description, the outcome is continuous so a count model such as Poisson or negative binomial would not make sense

Related Solutions

Solved – Data normalization for lme mixed model_R

A few points of general advice:

1) The plot you have of residuals vs. index isn't particularly useful for assessing the normality of the residuals. It is more useful to use a q-q plot or a histogram.

2) Don't rely on statistical tests (Shapiro-Wilk, Anderson–Darling, Kolmogorov–Smirnov, et al.) to determine if data or residuals are normally-distributed. They are sensitive to sample size. If you have a lot of data, they are likely to find a significant deviation from normal even if that deviation is small.

Here is a small example of a q-q plot and histogram in R. The chosen transformation doesn't work all that well for this example.

library(nlme)

data(mtcars)

mtcars$mpg_trans = mtcars$mpg ^ 1.3 

model = lme(mpg_trans ~ disp + cyl, random=~1|gear,
            data=mtcars,
            method="REML")

x = residuals(model)

hist(x, prob=TRUE, col="darkgray")
Range = seq(min(x), max(x), length = length(x))
Norm = dnorm(Range, mean = mean(x), sd = sd(x))
lines(Range, Norm, col = "blue", lwd = 2)

qqnorm(x)
qqline(x, col="red")

plot(predict(model), residuals(model))

Solved – Non-normality in linear mixed models/GLMM

There are a few points here:

The research question concerns the association of maximum dive depth by whales with whether the dives are of foraging or non-foraging type.

After fixing problems with autocorrelation, the main problem concerns the normality of residuals. Even with non normal residuals a model may still be useful but in order to make valid inferences, we would like them to be independent and close to normally distributed.

The QQ plot and histogram shows clear departures from normality. A GLMM was considered but given that the outcome variable is dive depth there is no obvious distribution candidate for a GLMM. On thing to consider is whether the outcome is bounded. Obviously it is bounded by zero, but presumably it is also bounded by the ocean floor ! With bounded data, a gamma GLMM could be considered, but based on the plots, it would seem appropriate to try to improve the fit of the LMM first. Also, it seems obvious that the upper bound on the sea depth will be variable and possibly unknown, so that's another reason to keep it simple.

Log-transforming the response has improved the QQ plot and histogram of residuals considerably and these can be considered as close to normal.

The plot of residuals vs fitted values appears a little strange, mostly because there is a small gap in the middle. There are very few fitted values between approximately 4.85 and 5.1 (on the log scale presumably), so the model is predicting very few dive depths in this range. Could there be any physical/biological reason for this ? Another (possibly related) explanation is an omitted binary/categorical variable (something related to dive depth like shallow/deep but probably not that simple/obvious). Anyway, I don't consider this to be much of a problem.

Update

Best Answer

Related Solutions

Solved – Data normalization for lme mixed model_R

Solved – Non-normality in linear mixed models/GLMM

Related Question