# spatial autocorrelation in dynamic occupancy models via colext

In biostatistics, spatial autocorrelation can almost always be an issue. In single-season models, one can apply mixed models to take into account spatial cluster of study plots. I wonder whether there is a similar way to deal with spatial autocorrelation in multi-season models.

Consider the following fake data featuring presence/absence data of a bird in 40 study plots. The plots are spatially clustered in three regions A, B and C. This is an example in R:

library (unmarked)
library (tidyverse)

M <- 40 # number of Sites
J <- 1 # num secondary sample periods
Times <- 2 # num primary sample periods

year1 <- c (1,1,1,0,0,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0)
year2 <- c (0,1,0,1,1,1,0,0,0,0,1,1,0,0,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0)

bird_territories <- array (data = c (year1, year2), dim = c (M, Times)) %>% as.matrix()

region <- c (rep ("A", 15), rep ("B", 12), rep ("C", 13))
var_a <- rnorm(40)
var_b <- (year1 - year2 + 0.1) * (rnorm (40)^2)
var_c <- data.frame ("recent" = (year2+0.1)*rnorm (40)^2,
"historic" = (year1+0.1)*rnorm (40)^2)

covariatesSite <- data.frame (region = region, var_a = var_a , var_b = var_b)
covariatesSiteYear <- list (var_c = var_c)

#create unmarked data frame
bird_data <- unmarkedMultFrame (y = bird_territories,
numPrimary = Times,
siteCovs = covariatesSite,
yearlySiteCovs = covariatesSiteYear
)

summary (bird_data)


fitting a null-model:

m0 <- colext (psiformula= ~ 1  , gammaformula = ~ 1, epsilonformula = ~ 1,
pformula = ~ 1, data = bird_data, method="BFGS")
summary(m0)


I'd be satisfied if there were no spatial autocorrelation. However, given that the study plots are arranged in three cluster, it might be sensible to include spatial autocorrelation. In a mixed model, this would look something like:

m0 <- colext (psiformula= ~ 1 + (1|region)  , gammaformula = ~ 1, epsilonformula = ~ 1,
pformula = ~ 1, data = bird_data, method="BFGS")
summary(m0)


However, I could not find such an option for the colext function. So my question is: does anyone know of a sensible approach to spatial autocorrelation in dynamic occupancy models?