# How should I deal with spatial autocorrelation in beta GLMM (glmmTMB)? (Bird diversity)

I am trying to compare certain a taxonomic diversity index (for bird communities) calculated for squares in a map grid system (a total of 34 squares) for two different years: 1998 and 2018. My responses (the diversity index) had values between 0 and 1, so I used beta distributed errors. Each square had to be compared to its future/past-self (obligatory - because of the way the index was calculated), so I defined square ID as a random effect (has 34 levels) and year as a fixed effect (2 levels).

My model structure was:

m1 <- glmmTMB(diversity_index ~ year + (1 | square_ID), family =
beta_family(), data = data)


Then, I realized that I have a spatial structure in my residuals, as shown by the distribution of residuals over space (usingbubble() function from sp package), variogram (using variogram() function from gstat package), and correlogram (correlog() function from ncf package) below. I produced different figures for the year 1998 and 2018, as I have two pieces of data for each point (one for 1998 and another for 2018 - I don't know if I was right to do that).

I have the coordinate (as longitude and latitude) data in data$longitude and data$latitude.

My question is what my model structure should look like if I were to model this data? Please keep in mind that my response is between (0,1) and each square needs to be compared to its own past/future-self (main aim: to tell whether site-based diversity increased or decreased in my study area). I'm open to any R package and any approach.