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 (using
bubble() 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
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.