###### Rt vs average temperature model #####
#Author:  
#Year: 2024

#The code runs an INLA model to analyze the association between _Lanice conchilega_ log reproductive rate 
#and average air temperature during winter.
#This model includes an SPDE to explicitly deal with spatial autocorrelation. 
#The SPDE is allowed to change between years with a "replicated" temporal structure.

# R version used: 4.2.0
# Rstudio version used: RStudio 2023.06.2+561 "Mountain Hydrangea" Release (de44a3118f7963972e24a78b7a1ad48b4be8a217, 2023-08-23) for windows
# Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) RStudio/2023.06.2+561 Chrome/110.0.5481.208 Electron/23.3.0 Safari/537.36

library(dplyr) # version 1.1.2
library(tidyr) # version 1.3.0
library(INLA)  # version 22.07.23


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

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

##### Mesh and SPDE ####

range_guess = 10 
max_edge = range_guess/5 
#transform m to km for easier handling
loc_complete <- cbind (rt_df$utm31n_x/1000, rt_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,
                       repl = rt_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.repl = n_temporal)

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

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

#for predictions
pred_mat <- expand.grid(
  xt = c(.5, 1, 1.5, mean(rt_df$xt)),
  flood_time_y = c(.27, .5, 1, mean(rt_df$flood_time_y)),
  t1_avg_winter_t = c(0:7, 4.4)) %>% 
  data.frame() %>% 
  mutate(flood_time_y.t1_avg_winter_t = flood_time_y * t1_avg_winter_t)


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

#stack for fitting
stack_1 <- inla.stack(
  tag = "fit",
  data = list(y = rt_df$rt),
  A = list (1, 1, A1),
  effects= list(
    Intercept = rep(1, nrow(covar_df)),
    X = covar_df,
    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 + 
  flood_time_y+
  t1_avg_winter_t + 
  flood_time_y.t1_avg_winter_t + #interaction term
  #spatial
  f(w, model = spde,
    replicate = w.repl) 

##### Run INLA model ####
rt_temperature_repl <-  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(rt_temperature_repl)