###### Richness vs Xt model #####
#Author: 
#Year: 2024

#The code runs an INLA model to analyze the association between species richness
#and _Lanice conchilega_ log abundance.
#This model includes an SPDE to explicitly deal with spatial autocorrelation. 
#The SPDE is allowed to change between years with an "AR 1" temporal structure.

# R version used: 4.2.2
# Rstudio version used: RStudio 2022.12.0+353 "Elsbeth Geranium" Release (7d165dcfc1b6d300eb247738db2c7076234f6ef0, 2022-12-03) for CentOS 7 Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/123.0.0.0 Safari/537.36

library(dplyr) # version 1.1.4
library(tidyr) # version 1.3.1
library(INLA)  # version 23.04.24


##### Data ####
spc_df <- read.csv("code/richness_data.csv") %>%  #poner el doi?
  #ensure year_chr is categorical
  mutate(year_chr = as.character(year_chr))

n_temporal <- length(unique(spc_df$year_chr))

###### Mesh specs #####
range_guess = 10 
max_edge = range_guess/5 

#transform m to km for easier handling
loc_complete <- cbind (spc_df$utm31n_x/1000, spc_df$utm31n_y/1000)

non_conv_hull <- inla.nonconvex.hull(loc_complete, convex =-0.04)

#prepare mesh
mesh1 <-  inla.mesh.2d ( boundary = non_conv_hull,
                         loc = loc_complete,
                         max.edge = c(1, 5)*max_edge,
                         offset =  c(0.5, 10),
                         cutoff = max_edge/5
)

A1 <- inla.spde.make.A(mesh1,
                       group = spc_df$ord_year,
                       loc = loc_complete)

####### SPDE & spatial field ######

spde <- inla.spde2.matern(mesh1, alpha = 2)

w_index <-  inla.spde.make.index(name = "w",
                                 n.spde = mesh1$n,
                                 n.group = n_temporal)

#### Covariable's data frame ####

#for fitting
covar_df <- model.matrix(~ - 1 + xt ,  
                         data = spc_df) %>% 
  data.frame()

#for predictions
pred_mat <- expand.grid(xt = seq(min(spc_df$xt), max(spc_df$xt), length.out = 20),
                        year_chr = "2019") %>% 
  data.frame() 


#### INLA Stack ####
#puts elements together

#stack for fitting
stack_1 <- inla.stack(
  tag = "fit",
  data = list(y = spc_df$sp_count),
  A = list (1,1,1, A1),
  effects= list(
    Intercept = rep(1, nrow(covar_df)),
    X = covar_df,
    year_chr = spc_df$year_chr,
    w = w_index)
)


#stack for prediction
stack_pred1 <- inla.stack(
  tag = "pred",
  data = list(y = NA),
  A = list (1,1),
  effects= list(
    Intercept = rep(1, nrow(pred_mat)),
    X = pred_mat)
)


all_stack <- inla.stack(stack_1, stack_pred1)


##### Formula ####
spatial_formula1 <- y ~ - 1 + Intercept + 
  xt + 
  #spatial
  f(w, model = spde,
    group = w.group,
    control.group = list (model = "ar1"))+
  #year as random factor
  f(year_chr, model = 'iid')

##### Run INLA model ####

richness_ar1 <-  inla(spatial_formula1,
                      family = "gaussian",
                      data = inla.stack.data(all_stack),
                      control.compute = list(
                        dic = TRUE,
                        waic = TRUE),
                      control.predictor = list(
                        A = inla.stack.A(all_stack)))


##### Results summary ####
summary(richness_ar1)