# How to account for temporal autocorrelation in logistic regression with longitudinal data?

I am attempting to perform a logistic regression on longitudinal data (game camera footage of nesting birds, with a photo taken every 5 minutes for a period of 10-28 days, depending on whether the nest was abandoned part way through the nesting period or not). In total I have ~25,000 observations from 10 separate nests.

My overall goal is to determine the temperature at which these birds begin experiencing heat stress. I really don't care about the effect of time, I just want to control for the autocorrelation caused by the fact that this data is from a time-series. I've classified each photo in a binary manner according to whether the bird on the nest is displaying thermoregulatory behavior (0=no thermoregulatory behavior displayed, 1=thermoregulatory behavior displayed). For each photo, I have simultaneously collected weather station data, so that I can say "When it was X degrees out, the bird on the nest was displaying/not displaying thermoregulatory behavior." I have other climate data as well (humidity, wind speed, precipitation rate, etc.) that goes with each photo. I've used a hierarchical model selection approach (using mixed effects models to control for the random effect of nest site, since I have multiple photos of the individual at each site). I've found that my best model is:

Model8 <- glmer(ThermalResponse ~ Temperature + Humidity + (1 | Location),
data = thermoreg, family = binomial())


I've run into a few problems: Ideally, I would like to control for nesting territory (since I have multiple observations from each individual). However, there is definitely autocorrelation in my residuals (see plot below).

I've tried multiple approaches to correct for this, but from what I've read, it is not possible to incorporate both autocorrelation and a random effect into the same model. I tried doing so using the lme function from the nlme package (THe model I used was: Model1<-lme(ThermalResponse~Temperature+Humidity, data=thermoreg, random=~1|Location, correlation=corAR1(form~1|TimeSeries)) , but when I specified the model, it ran and ran (for over 24 hours!) and I eventually canceled it.

I've also tried the bild function (from the bild package), but I've received several errors that I cannot find documentation for anywhere.

In addition to those approaches, I've tried changing the structure of my data-instead of using the individual photos, I've averaged the response and predictor variables over different time periods (1 hour, 2 hours, 4 hours) and performed a proportional logistic regression in order to try to remove the autocorrelation, however it did not help with the autocorrelation (see the plots below, which show the acf plot for the 2-hour and 4-hour proportional logistic regression models).

I'm at a loss for what to do/try next. Any advice is greatly appreciated!

• How did you make the auto correlation plots? I guess with Pearson residuals in which case these should not be normally distributed when your outcomes are binary and you use logistic regression. Thus, I am not sure how informative the ACF plots are. – Benjamin Christoffersen Mar 17 at 9:24
• I used the acf function- so my code was acf(residuals(Model1)). I believe that uses the Pearson residuals. Is there another method of evaluating autocorrelation that might be more appropriate in this case? – C.H. Mar 18 at 13:42
• acf(residuals(Model1)) likely makes no sense. See my updated answer. – Benjamin Christoffersen Mar 18 at 16:54

If you would include the time variable in the specification of your random effects you would account for the correlations in the repeated measurements of your outcome variable ThermalResponse, e.g., something like

fm1 <- glmer(ThermalResponse ~ Temperature + Humidity + (TimeSeries | Location),
data = thermoreg, family = binomial())


You could study further evaluate if you the correlations in the data are more complex by including nonlinear terms of time in the specification of random-effects part, for example, using a second-degree polynomial , e.g.,

fm2 <- glmer(ThermalResponse ~ Temperature + Humidity + (poly(TimeSeries, 2) | Location),
data = thermoreg, family = binomial())


and compare with the previous model using a likelihood ratio test, i.e.,

anova(fm1, fm2)


Note that both glmer() and the glmmTMB() use the Laplace approximation in the calculation of the likelihood function of the model, which is known not to be that optimal for binary data. A better approach is to use the adaptive Gaussian quadrature approximation provided, for example, by the GLMMadaptive package.

• I believe he also concerned about temporal autocorrelation... GLMM ADAPTIVE can't apply different covariance structures. Unless I am wrong, obviously you would be the most likely to know. – OliverFishCode Mar 5 at 19:39
• Indeed GLMMadaptive does not yet have an auto-regressive structure, but my point was that by including the time variable as a random effect you similarly account for temporal autocorrelation. – Dimitris Rizopoulos Mar 5 at 19:41

If I understand you correctly, You want to argue that there is an effect of the temperature controlling for time and the location. Thus, you could use a fixed effect interaction between time and the location like in this example

n_times <- 100L # number of time periods
n_sites <- 10L  # number of sites
n_sub   <- 2L   # number of subjects / site

#####
# simulate data
set.seed(58550574)
df <- lapply(1:n_sites, function(site) {
# temporal site effect
ti <- Reduce(function(x, y) .5 * x + y, rnorm(n_times), accumulate = TRUE)
z <- runif(1, -1, 1) # location effect
do.call(rbind, lapply(1:n_sub, function(id) {
x <- runif(n_times, -1, 1) # variable of interest
y <- (1 + exp(-(z + .5 * x + ti)))^(-1) > runif(n_times)
data.frame(id = id, x = x, y = y, site = site, time = 1:n_times)
}))
})
df <- within(do.call(rbind, df), {
site <- as.factor(site)
time <- as.factor(time)
})

#####
# fit model and show variable of interest with Z-score
fit <- glm(y ~ x + site * time, binomial(), df)
summary(fit)$coefficients["x", , drop = FALSE] #R Estimate Std. Error z value Pr(>|z|) #R x 0.833 0.183 4.55 5.43e-06 # not sure this is a good idea but we do the plots anyway to show # that there is not sign of auto correlation as expected par(mar = c(5, 4, .5, .5), mfcol = c(3, 3)) tapply( residuals(fit, type = "pearson"), list(df$$site, df$$id), acf, main = "") # fit model faster library(speedglm) fit <- speedglm( y ~ x + site * time, family = binomial(), data = df, sparse = TRUE) summary(fit)$coefficients["x", , drop = FALSE]
#R   Estimate Std. Error z value Pr(>|z|)
#R x    0.833      0.183    4.55 5.43e-06


My overall goal is to determine the temperature at which these birds begin experiencing heat stress. I really don't care about the effect of time, I just want to control for the autocorrelation caused by the fact that this data is from a time-series.

The question is how much variation there is in the temperature which is unrelated to location and time. This seems like something that might be an issue regardless of the model you use.

## Update

It seems like you just did acf(residuals(<symbol for glmer output>)). What you get out of this will likely not make any sense and depends on the order of the original data. The code below shows an example hereof (this was too short to put into a comment)

library(lme4)
# ordered by events
cbpp <- cbpp[order(cbpp$$incidence/cbpp$$size), ]
gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), cbpp,
binomial, nAGQ = 0)
acf(residuals(gm1a))


# random order
set.seed(6490517)
cbpp <- cbpp[sample.int(nrow(cbpp)), ]
gm1a <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd), cbpp,
binomial, nAGQ = 0)
acf(residuals(gm1a))


What I think you want/intended to is to plot the within location auto-correlation of the Pearson residuals. Some caution may be needed in the interpretation of the plot as you observe binary events which does not justify a normal approximation of the Pearson residuals.

Switch to GLMMTMB package it allows for different covariance structures like AR1 with random effects and the use of multiple different distribution including the binomial. I use it all the time for logistic mixed models with temporal autocorrelation. Also, the reason LME freaked out is that your data are not normal (which it requires); thus, why you are using the binomial distribution.

https://rdrr.io/cran/glmmTMB/man/glmmTMB.html

Here is an example of someone using an AR1 in GLMMTMB on a random effect:

https://github.com/glmmTMB/glmmTMB/issues/329

Goodluck!

• Thank you for this response! I have a followup question, if you don't mind. In GLMMTMB, I fit the model as: Model<-glmmTMB(ThermalResponse~Temperature+Humidity+(1| Location)+ar1(TimeSeries-1|Location), family=binomial, data=thermoreg) However, I am getting the following error: Warning message: Model convergence problem; non-positive-definite Hessian matrix. This warming shows up even when I remove the fixed effect humidity and the random effect location. I think my problem might be the format my time series data is in. Currently, it ranges from 1-9611. Is this range too large? – C.H. Mar 4 at 8:35
• Are you referring to the response, grouping or time as 1-9611 – OliverFishCode Mar 4 at 13:42
• I assume it's time step...i.e. you have 9611 pictures – OliverFishCode Mar 4 at 14:17
• Yes, that's correct. I have 9611 pictures. – C.H. Mar 4 at 17:33
• Or rather, in my site with the most pictures, there are 9611 pictures. I originally assigned a time series based on the number of minutes that had elapsed since the first photo at any site was taken (i.e. if site 1 started at 12 PM on 6/1/2018 and site 2 started at 12 PM on 6/2/2018, the time series value for the first photo taken at site 1 was "1" and the time series value for the first photo from site 2 was 1441. As you can imagine, the range of values when I did it that way was VERY large, so I then assigned values within groups rather than overall and did it as a time step. – C.H. Mar 4 at 17:41