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!


 A: 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. 
A: 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.
Link with explanation of package:
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!
A: 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. 
