Inspecting and modeling residual autocorrelation with gaps in linear mixed effects models

Question

Here I generate a dataset where measurements of response variable y and covariates x1 and x2 are collected on 30 individuals through time. Each individual is denoted by a unique ID. The observations are collected in hourly increments, but are only available for given individuals (IDs) when they are present during the respective hour (thus creating irregularities in each time series).

library(tidyverse)
library(lubridate)
library(data.table)
set.seed(123)
TimeSeries <- data.table(tm = 
               rep(
               seq(as.POSIXct("2021-01-23 01:00"), as.POSIXct("2021-10-27 17:00"), by="hour"),30),
             ID = factor(rep( paste("ID_",c(1:30), sep = ""), each = 6664)),
             Obs = sample(c(NA, 1), 199920,prob = c(0.7,0.3), replace = TRUE))
TimeSeries<- TimeSeries[Obs == 1]
#explicitly making large gaps at the beginning of some time series to illustrate that some individuals are present (for  the first time) until later in the time series  
TimeSeries <- 
TimeSeries%>%
  dplyr::filter(!(ID== "ID_2" & tm < as.POSIXct("2021-5-27 10:00")), 
         !(ID== "ID_5" & tm < as.POSIXct("2021-3-10 15:00")),
         !(ID== "ID_6" & tm > as.POSIXct("2021-3-10 15:00")),
         !(ID== "ID_5" & tm > as.POSIXct("2021-6-10 23:00")))
#response variable
TimeSeries[,y:= rnorm(nrow(TimeSeries))]
#predictors x1 and x2
TimeSeries[,x1:= rnorm(nrow(TimeSeries))]
TimeSeries[,x2:= rnorm(nrow(TimeSeries))]
#now irrelevant so remove:
TimeSeries[,Obs:= NULL]

we wish to fit a linear mixed effects model to determine if changes in x1 and x2 have an effect on the response y while allowing for variation across indivduals with a random intercept. I demonstrate this with nlme:

mod <- lme(y ~ x1+x2, random = ~1|ID, data = TimeSeries, method = "ML")

However, we suspect that when IDs are present in consecutive (or close) hours, the residuals will be heavily autocorrelated (within IDs).

Thus we would like to check this assumption with ACF/PACF plots, and explore different correlation structures if it is a problem. (note, obviously this will not be true for the simulated data above, it was just to illustrate the structure of my data)

I am unsure of the appropriate method to look at autocorrelation in this case, and calculate confidence bands. My understanding is that the nlme::ACF function will respect the grouping structure of random effects, but does not calcualte the correct autocorrelation function with irregular or missing values (the later of which is the apparent issue here). Is this true even with the inclusion of na.action = na.omit in the ACF call? Or is there a more appropriate method?

Moreover, assuming autocorrelation is an issue, does the structure of this data require the use of continuous correlation structures (e.g., corCAR1) or is it reasonable to use corAR1 and corARMA if we do something like specify the time variable as the number of hours since the first observation? for example:

TimeSeries[,TimeSinceFirst:= as.numeric(difftime(tm, min(TimeSeries$tm), units = "hours"))]
lme(y ~ x1+x2, random = ~1|ID, correlation = corAR1(form = ~TimeSinceFirst | ID),data = TimeSeries, method = "ML")

I have worked through several examples (including those available here, here , and [here] (https://bbolker.github.io/mixedmodels-misc/ecostats_chap.html) ), but I cant seem to find any that deal with gaps and irregularities within each time series like I have presented above.

Answer 1

There are two different questions at play here

1. How to model the autocorrelation function, $\rho(t)$, for irregularly spaced data
2. How to adjust for autocorrelation for irregularly spaced data

You've somewhat answered the second question with the use of corAR1, but note that this still assumes the autocorrelation function $\rho(t) = \alpha^t$ for any $t \in \mathbb{R}$ (instead of $t \in \mathbb{Z}$), and can still potentially face model misspecification.

Question 1

Assuming stationarity, I personally like to work with the semi-variogram function, $\gamma(t)$, instead: \begin{align*} \gamma(t) = \frac{1}{2}\mathbb{E}[(R_{t+s} - R_s)^2] \end{align*} where $R(t)$ is the residual at time $t$ within the same ID. Let $\sigma^2 = \text{Var}(R_t)$ be the common variance at all time points (implied by stationarity), we have \begin{align*} \rho(t) = 1 - \frac{\gamma(t)}{\sigma^2} \end{align*} To estimate $\gamma(t)$, we can fit a smoothing curve to the observed half-squared differences between residuals vs the time differences: \begin{align*} \widehat{\gamma}_{i,jj'} = \frac{1}{2}(R_{ij} - R_{ij'})^2 \quad \text{vs} \quad t_{i,jj'} = |\tau_{ij} - \tau_{ij'}| \end{align*} where $\tau_{ij}$ are the real times at which the $j$th measurement on the $i$th individual is observed.

library(nlme)
library(npreg)

#Half pairwise differences of residuals
fit_conditional = lme(y ~ x1+x2, random = ~1|ID, data = TimeSeries, method = "ML")
TimeSeries$resid = resid(fit_conditional)
semivariogram.data = TimeSeries %>% group_by(ID) %>% summarize(gamma = c(dist(resid)^2/2),
                                          t = c(dist(tm)))

#Normalize time scale to hours
semivariogram.data$t = semivariogram.data$t/3600

#Residual variance
s2 = var(TimeSeries$resid)

#Smoothing spline
mod.smooth = ss(semivariogram.data$t, semivariogram.data$gamma, nknots = 10)

#Predict smoothed semivariogram
myfit = predict(mod.smooth, x = 1:6663, se = F)
names(myfit) = c("t", "gamma")

#Compute autocorrelation function from this
myfit$rho = 1 - myfit$gamma/s2
ggplot(myfit, aes(x=t, y=rho)) +
  geom_line() +
  ylim(-1,1)

The autocorrelation here is essentially zero except near the very end, where variability is introduced due to the fewer pairs of half-squared differences.

Question 2

Besides the lme approach, you can also consider a GEE approach, where even if you get the correlation structure wrong, you're guaranteed correct standard errors as long as your mean model is correctly specified. So, we can still adjust for an AR(1) structure, but if it's misspecified, no worries.

library(geepack)
mod = geeglm(y ~ x1+x2, id=ID, data=TimeSeries, corstr="ar1", waves=tm)