AIM of script

Fit non-spatial, spatial and spatio-temporal INLA models Raja clavata in the Greater North Sea, run predictions for best model and prepare best model for simulation of predictions.

Setup

clear global environment

rm(list = ls()); gc()
getwd()

load libraries

library(tidyverse)
library(INLA)
library(inlabru)
library(sf)
library(ggplot2)
library(rnaturalearth)

Set date for files surveys

Timestamp when saving dowloaded exchange data.

dl_date <- "14_06_2024"

Setup file paths:

  • data_path: file path to data
  • output_path: file path for outputs from this script
  • model_output_path: file path for fitted models
  • results_path: file path for results

Supporting fuctions

Custom functions that help summarise model outputs.

## flextable autofit to width
# function to autofit flextable to desired word page width in cm
FitFlextableToPage <- function(ft, pgwidth = 17){ # desired pgwidth in cm
  ft_out <- ft %>% autofit()
  ft_out <- flextable::width(ft_out, width = dim(ft_out)$widths * pgwidth * 0.3937 /(flextable_dim(ft_out)$widths))
  return(ft_out)
}

Load data required for fitting INLA models:

  • model_data: processed data of Datras haul and catch survey data by size class for Raja clavata
  • NS_EEC_poly: polygons as simple_features for Greater North Sea study area
  • NS_EEC_poly_simpl_UTM_km: same as above NS_EEC_poly but with simplified coastlines
  • NS_EEC_simpl_EUNIS_2021: simple feature collection of EUNIS_2021 for study area
  • bathymetry_GNS: Bathymetry data for Greater North Sea study area for prediction grid
  • gear_efficiency: gear efficiency data for Raja clavata obtained from Walker et al. 2017: https://doi.org/10.1093/icesjms/fsw250
  • exploit_maturity_length_class_stats: data of exploitation_maturity mean weight, mean length, gear, efficiency
  • exploit_length_class_stats: data of exploitation mean weight, mean length, gear, efficiency
# To load the data again
load(file.path(data_path, "data_RJC_INLA.RData"))

get final year of data that is used for naming results files

last_dat_yr <- max(model_data$Year)

set colour pallettes

# 18 distinct colours from 99% 
cb18_99_pal <- c('#e6194B', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#42d4f4', '#f032e6', '#fabed4', '#469990', '#dcbeff', '#9A6324', '#fffac8', '#800000', '#aaffc3', '#000075', '#a9a9a9', '#ffffff', '#000000')

Load North Sea polygon and GNS_countries polygons

NS_EEC_poly <-  read_rds(file.path(data_path, "NS_EEC_poly.RDS")) %>% st_union() %>% st_set_crs(3857) %>% st_transform(crs = "+proj=utm +zone=31 +units=km") # transform to UTM 31N with km as units

GNS_countries_UTM_km <- ne_countries(continent = "Europe", scale = "medium") %>% st_as_sf() %>% st_transform(crs = "+proj=utm +zone=31 +units=km") # transform to UTM 31N with units as km

Create meshes for INLA models

create nonconvex mesh: ConvHull

Loc <- with(model_data[!duplicated(model_data$HaulCode),],cbind(XkmRounded, YkmRounded))

## JJ simple mesh 
ConvHull_JJ <- inla.nonconvex.hull(points = Loc, convex = -0.02 , resolution = 55)

nonconvex_hull_JJ_poly <- ConvHull_JJ$loc %>% 
  as_tibble() %>% 
  dplyr::rename(X = V1, Y = V2) %>% 
  st_as_sf(coords = c("X", "Y")) %>% 
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON") %>% 
  st_set_crs("+proj=utm +zone=31 +units=km")

mesh_ConvHull_JJ <- inla.mesh.2d(boundary = ConvHull_JJ, max.edge = 45, cutoff = 10)
mesh_ConvHull_JJ$n

## plot mesh with with land mass
model_data %>% distinct(HaulCode, .keep_all = TRUE) %>% 
ggplot() + 
geom_point(aes(x = XkmRounded, y = YkmRounded, colour = Survey)) +
geom_sf(data = NS_EEC_poly, colour = "green", fill = NA) +                      #plot NS_EEC_land boundary shape
geom_sf(data = GNS_countries_UTM_km) +                                                       # plot land
gg(mesh_ConvHull_JJ, edge.color = "royalblue1", ext.color = "royalblue1", alpha = 0.3) + #plot mesh
coord_sf(xlim = c(-420, 1400),
         ylim = c(5500, 7000)) +
theme_classic()

create mesh with nonconvex + land boundary and outer bufferzone

create mesh based on landboundary with outerbuffer based on previous non-convex mesh (this ensures we are not digging into Norwegian unsampled area but follow land boundaries where possible)

# Best guess of range: formed from distance histogram & variogram & biology of study animals
range_guess <- 130 #

## mesh with land as boundary + Outer bufferzone
bound_outer <- range_guess # should be at least as large as the range (distance at which the dependency diminishes)

# assign mesh parameters
max_edge_UTM_km <- range_guess/5; max_edge_UTM_km

# create combination of non-convex and land-boundary and save to file
nonconvex_landboundary_poly <- nonconvex_hull_JJ_poly %>% st_intersection(NS_EEC_poly) #NS_EEC_poly_simpl_UTM_km

mesh_nonconvex_landboundary_outerbuffer_UTM_km <- inla.mesh.2d(
  loc = Loc,
  boundary = as_Spatial(nonconvex_landboundary_poly), 
  max.edge = c(1, 3) * max_edge_UTM_km, 
  cutoff = 20,  # min distance at which points get unique vertices (if within the distance they share a vertice)
  offset = c(max_edge_UTM_km, bound_outer))

mesh_nonconvex_landboundary_outerbuffer_UTM_km$n
plot(mesh_nonconvex_landboundary_outerbuffer_UTM_km, asp = 1)

#write_rds(mesh_nonconvex_landboundary_outerbuffer_UTM_km, file = file.path(results_path, "mesh_nonconvex_landboundary_outerbuffer_UTM_km.rds"))

## plot mesh with with land mass
model_data %>% distinct(HaulCode, .keep_all = TRUE) %>% 
  ggplot() + 
  geom_point(aes(x = XkmRounded, y = YkmRounded), colour = "red") +
 
  geom_sf(data = NS_EEC_poly, colour = "green", fill = NA) +                      #plot NS_EEC_land boundary shape
  geom_sf(data = GNS_countries_UTM_km, alpha = 0.2) + # plot land
  gg(mesh_nonconvex_landboundary_outerbuffer_UTM_km, edge.color = "royalblue1", ext.color = "royalblue1", alpha = 0.3) + #plot mesh
  geom_sf(data = NS_EEC_poly, fill = "orange", alpha = 0.8) +
  coord_sf(xlim = c(-420, 1400),
           ylim = c(5500, 7000)) +
  theme_classic() 

INLA models

Unexploited_Immature = Size class 1 Exploited_Immature = Size class 2 Exploited_Mature = Size class 3

################################### Stack based syntax for non-spatial INLA model with 3 RWs for Year specific to Exploit_Maturity classes
NoPerHaul <- model_data$NoPerHaul
Unexploited_Immature <- model_data$Unexploited_Immature
Exploited_Immature <- model_data$Exploited_Immature
Exploited_Mature <- model_data$Exploited_Mature
Year_Unexploited_Immature <- model_data$Year_Unexploited_Immature
Year_Exploited_Immature <- model_data$Year_Exploited_Immature
Year_Exploited_Mature <- model_data$Year_Exploited_Mature
logSweptAreakm2 <- model_data$logSweptAreakm2
Depth_Unexploited_Immature <- model_data$Depth1m_Unexploited_Immature
Depth_Exploited_Immature <- model_data$Depth1m_Exploited_Immature
Depth_Exploited_Mature <- model_data$Depth1m_Exploited_Mature

non-spatial model

IM_NS1 Stack based syntax 3RW Year

model_name <- "IM_NS1_no_priors"
description_IM_NS1 <- "non-spatial 3RWs of Year Exploit_Maturity specific"

################################### Stack based syntax for non-spatial INLA model with 3 RWs for Year specific to Exploit_Maturity classes
# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity * Exploit_Maturity, data = model_data)
X <- as.data.frame(X[,-1]) # remove Intercept 
names(X)<- c(gsub("[// |//_ |//:]", "",names(X))) # remove special characters

# create Stack                    )
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature), 
                                       1, 
                                       1),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature     = Year_Exploited_Immature,
                                       Year_Exploited_Mature       = Year_Exploited_Mature, 
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2
                                       )
                        )

formula <- as.formula(paste('NoPerHaul ~ -1 + Intercept',
                             paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates,
                             'f(Year_Unexploited_Immature,  model="rw1", scale.model = TRUE) + f(Year_Exploited_Immature,  model="rw1", scale.model = TRUE) + f(Year_Exploited_Mature,  model="rw1", scale.model = TRUE)',
                             'offset(logSweptAreakm2)',
                             sep = " + " )
                       ); formula

# save INLA model setup objects
model_setup <- list("description" = description_IM_NS1,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = NA, 
                                       "mesh" = NA, 
                                       "y_mesh" = NA, 
                                       "A" = NA, 
                                       "spde_name" = NA, 
                                       "spde" = NA, 
                                       "w.st" = NA, 
                                       "X" = X, 
                                       "Stack" = mystack)

# fit IM
rm(IM)
IM = inla(formula,
           data              = inla.stack.data(mystack),
           family            = "nbinomial",
           control.predictor = list(A = inla.stack.A(mystack), compute = TRUE),
           control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #
           control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
           #control.mode      = list(fixed=T, theta = init_hyper), # set initial values for hyperparameters
           verbose           = TRUE,
           inla.mode         = "experimental")
 
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0

# save inla model & model setup
IM_NS1 <- IM

IM_filename <- paste(model_name, "3RWsYear_1988",last_dat_yr,".rds", sep="_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

summary(IM)

plots_RandomEffects <- list()
for(i in 1:length(IM$summary.random)){
  covariate <- names(IM$summary.random)[i]
  plots_RandomEffects[[covariate]] <- IM$summary.random[[covariate]] %>%
  ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(sides="b") +
  theme_bw() +
  ggtitle(covariate)
}

patchwork::wrap_plots(plots_RandomEffects) + patchwork::plot_layout(nrow = 3, ncol =1)

IM_NS2 Stack based syntax 3RW Year + 3RW Depth

model_name <- "IM_NS2_no_priors"
description_IM_NS2 <- "non-spatial 3RWs of Year Exploit_Maturity specific & 3 RWs of Depth Exploit_Maturity specific"

################################### Stack based syntax for non-spatial INLA model with 3 RWs for Year specific to Exploit_Maturity classes
# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity * Exploit_Maturity, data = model_data)
X <- as.data.frame(X[,-1]) # remove Intercept 
names(X)<- c(gsub("[// |//_ |//:]", "",names(X))) # remove special characters

# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature,
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2
                                       )
                        )

## scaled rw
formula = as.formula(str_c("NoPerHaul",
                    str_c( "-1 + Intercept +",
                      paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                      paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                      #paste0('f(w, model = spde, group = w.group, control.group = list(model = "ar1"))'), 
                      paste0('offset(logSweptAreakm2)'),
                      sep = " + "), sep = " ~ ")); formula

# save INLA model setup objects
model_setup <- list("description" = description_IM_NS2,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = NA, 
                                       "mesh" = NA, 
                                       "y_mesh" = NA, 
                                       "A" = NA, 
                                       "spde_name" = NA, 
                                       "spde" = NA, 
                                       "w.st" = NA, 
                                       "X" = X, 
                                       "Stack" = mystack)

# fit IM
rm(IM)
IM = inla(formula,
           data              = inla.stack.data(mystack),
           family            = "nbinomial",
           control.predictor = list(A = inla.stack.A(mystack), compute = TRUE),
           control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #
           control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
           verbose           = TRUE,
           inla.mode         = "experimental")
 
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0 # should be TRUE -> different to 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ # as long as this equals 0 rerun the model 
  (i = i+1)
  IM <- inla.rerun(IM)
}

# save inla model & model setup
IM_NS2 <- IM

IM_filename <- paste(model_name, "3RWsYear_1988",last_dat_yr,".rds", sep="_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

summary(IM)

plots_RandomEffects <- list()
for(i in 1:length(IM$summary.random)){
  covariate <- names(IM$summary.random)[i]
  plots_RandomEffects[[covariate]] <- IM$summary.random[[covariate]] %>%
  ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(sides="b") +
  theme_bw() +
  ggtitle(covariate)
}

patchwork::wrap_plots(plots_RandomEffects) + patchwork::plot_layout(nrow = 3, ncol =2)

IM_NS3 Stack based syntax 3RW Year + 3RW Depth + EUNIS2019C stack approach

model_name <- "IM_NS3_no_priors"
description_IM_NS3 <- "non-spatial 3RWs of Year Exploit_Maturity specific & 3 RWs of Depth Exploit_Maturity specific and EUNIS2019C"

################################### Stack based syntax for non-spatial INLA model with 3 RWs for Year specific to Exploit_Maturity classes
# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity * Exploit_Maturity + fEUNIS2019C, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//_ |//:]", "",names(X))) # remove special characters
head(X)

# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                       A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2
                                       )
                        )

formula = as.formula(str_c("NoPerHaul",
                    str_c( "-1 + Intercept +",
                      paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                      paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                      paste0('offset(logSweptAreakm2)'),
                      sep = " + "), sep = " ~ ")); formula

# save INLA model setup objects
model_setup <- list("description" = description_IM_NS3,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = NA, 
                                       "mesh" = NA, 
                                       "y_mesh" = NA, 
                                       "A" = NA, 
                                       "spde_name" = NA, 
                                       "spde" = NA, 
                                       "w.st" = NA, 
                                       "X" = X, 
                                       "Stack" = mystack)

# fit IM
rm(IM)
IM = inla(formula,
           data              = inla.stack.data(mystack),
           family            = "nbinomial",
           control.predictor = list(A = inla.stack.A(mystack), compute = TRUE),
           control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #
           control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
           verbose           = TRUE,
           inla.mode         = "experimental")
 
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0 # should be TRUE -> different to 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while(sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ # as long as this equals 0 rerun the model 
  (i = i+1)
  IM <- inla.rerun(IM)
}

# save inla model & model setup
IM_NS3 <- IM

IM_filename <- paste(model_name, "3RWsYear_1988",last_dat_yr,".rds", sep="_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

summary(IM)

plots_RandomEffects <- list()
for(i in 1:length(IM$summary.random)){
  covariate <- names(IM$summary.random)[i]
  plots_RandomEffects[[covariate]] <- IM$summary.random[[covariate]] %>%
  ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(sides="b") +
  theme_bw() +
  ggtitle(covariate)
}

patchwork::wrap_plots(plots_RandomEffects) + patchwork::plot_layout(nrow = 3, ncol =2)

IM_NS4 Stack based syntax 3RW Year + 3RW Depth + fSubstrate stack approach

model_name <- "IM_NS4_no_priors"
description_IM_NS4 <- "non-spatial 3RWs of Year Exploit_Maturity specific & 3 RWs of Depth Exploit_Maturity specific and Substrate as provided  by eunis data"

################################### Stack based syntax for non-spatial INLA model with 3 RWs for Year specific to Exploit_Maturity classes
# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity * Exploit_Maturity + fSubstrate, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//:]|[[]|[]]", "",names(X))) # due to " " (space) & [ & ] in Substrate factor levels (spaces in factor levels)
head(X)

# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                       A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2
                                       )
                        )

# define formula without any prior
formula = as.formula(str_c("NoPerHaul",
                    str_c( "-1 + Intercept +",
                      paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                      paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                      paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                      paste0('offset(logSweptAreakm2)'),
                      sep = " + "), sep = " ~ ")); formula

# save INLA model setup objects
model_setup <- list("description" = description_IM_NS4,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = NA, 
                                       "mesh" = NA, 
                                       "y_mesh" = NA, 
                                       "A" = NA, 
                                       "spde_name" = NA, 
                                       "spde" = NA, 
                                       "w.st" = NA, 
                                       "X" = X, 
                                       "Stack" = mystack)

# fit IM
rm(IM)
IM = inla(formula,
           data              = inla.stack.data(mystack),
           family            = "nbinomial",
           control.predictor = list(A = inla.stack.A(mystack), compute = TRUE),
           control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #
           control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
           verbose           = TRUE,
           inla.mode         = "experimental")
 
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0 # should be TRUE -> different to 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ # as long as this equals 0 rerun the model 
  (i = i+1)
  IM <- inla.rerun(IM)
}

# save inla model & model setup
IM_NS4 <- IM

IM_filename <- paste(model_name, "3RWsYear_1988",last_dat_yr,".rds", sep="_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

summary(IM)

plots_RandomEffects <- list()
for(i in 1:length(IM$summary.random)){
  covariate <- names(IM$summary.random)[i]
  plots_RandomEffects[[covariate]] <- IM$summary.random[[covariate]] %>%
  ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
  geom_line() +
  geom_ribbon(alpha = 0.2) +
  geom_rug(sides="b") +
  theme_bw() +
  ggtitle(covariate)
}

patchwork::wrap_plots(plots_RandomEffects) + patchwork::plot_layout(nrow = 3, ncol =2)

spatial models (mesh_nonconvex_landboundary_outerbuffer_UTM_km)

set up mesh and w for spatial models

# assign mesh
mesh_name <- deparse(substitute(mesh_nonconvex_landboundary_outerbuffer_UTM_km))
mesh <- mesh_nonconvex_landboundary_outerbuffer_UTM_km

# Compute projection matrix
obsloc <- with(model_data,cbind(XkmRounded, YkmRounded))
A <- inla.spde.make.A(mesh = mesh, loc = obsloc)
dim(A)

# Defining spatial field w
spde <- inla.spde2.matern(mesh)

# define index vector for spde
w.st <- inla.spde.make.index('w', n.spde = spde$n.spde)

IM_S5 spatial model 3 RWs Year stack approach

model_name <- "IM_S5_no_priors"
description_IM_S5 <- "non-spatial 3RWs of Year Exploit_Maturity specific"

# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity, data = model_data)
X <- as.data.frame(X[,-1]) # remove Intercept 
names(X)<- c(gsub("[// |//_ |//:]", "",names(X))) # remove special characters

# design matrix
# create design matrices for covariates
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                       effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                                      )
                       ) 

formula = as.formula(str_c("NoPerHaul",
                    str_c( "-1 + Intercept +",
                      paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                      paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                      paste0('f(w, model = spde)'), 
                      paste0('offset(logSweptAreakm2)'),
                      sep = " + "), sep = " ~ ")); formula

model_setup <- list("description" = description_IM_S5,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = NA, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X,
                                       "Stack" = mystack)

# fit inla model
IM = inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A = inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE,  config = TRUE),
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               #control.mode      = list(restart=T, theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode = "experimental")

# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0

i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

################
# save inla model & model setup
IM_S5 <- IM

IM_filename <- paste(model_name,"3RWsYear", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

############
summary(IM)

IM_S6 spatial model 3 RWs Year + 3RWs Depth stack approach

model_name <- "IM_S6_no_priors"
description_IM_S6 <- "non-spatial 3RWs of Year & 3RWs of Depth1m Exploit_Maturity specific"

# create design matrix for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity, data = model_data)
X <- as.data.frame(X[,-1]) # remove Intercept 
names(X)<- c(gsub("[// |//_ |//:]", "",names(X))) # remove special characters

# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                       effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                                      )
                       )   

formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept ",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde)'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

model_setup <- list("description" = description_IM_S6,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = NA, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh, 
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,                   
                                       "y_mesh" = NA, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X,
                                       "Stack" = mystack)


# fit inla model
IM = inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A = inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE),
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               #control.mode      = list(restart=T, theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode = "experimental")

# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

################
# save inla model & model setup
IM_S6 <- IM

IM_filename <- paste(model_name,"3RWsYear_3RWsDepth", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

############
IM <- IM_S6
summary(IM)

IM_S7 spatial model 3 RWs Year + 3RWs Depth + EUNIS2019C stack approach

starting values from successful fit IM_S6

model_name <- "IM_S7_no_priors"
description_IM_S7 <- "non-spatial 3RWs of Year & 3RWs of Depth1m Exploit_Maturity specific & EUNIS2019C"

# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity + fEUNIS2019C, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//: |//_]", "", names(X)))
head(X)
    
# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                       effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                                      )
                       )   

# define formula
# no priors
formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept ",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde)'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

init_hyper_succ_fit_IM_S6 <- c(log(IM_S6$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature", "Precision for Year_Exploited_Mature", "Precision for Depth_Unexploited_Immature", "Precision for Depth_Exploited_Immature", "Precision for Depth_Exploited_Mature"), "mean"]), IM_S6$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"])

init_hyper <- init_hyper_succ_fit_IM_S6

model_setup <- list("description" = description_IM_S7,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = NA, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X,
                                       "Stack" = mystack)

# fit inla model
rm(IM)
IM = inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A = inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE),
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               control.mode      = list(restart=T, theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode = "experimental")

#check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

################
# save inla model & model setup
IM_S7 <- IM

IM_filename <- paste(model_name, "3RWsYear_3RWsDepth_substrate", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

############
summary(IM)

IM_S8 spatial model 3 RWs Year + 3RWs Depth + Substrate stack approach

starting values from successful fit IM_S6

model_name <- "IM_S8_no_priors"
description_IM_S8 <- "non-spatial 3RWs of Year & 3RWs of Depth1m Exploit_Maturity specific & Substrate"

# create design matrices for covariates
X <- model.matrix( ~ fSurveyQuarter * Exploit_Maturity + fSubstrate, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//:]|[[]|[]]", "",names(X))) # due to " " (space) & [ & ] in Substrate factor levels (spaces in factor levels)
head(X)

# create Stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                       effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                                       )
                       )   


# define formula without any prior
formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept ",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde)'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

IM_S6$summary.hyperpar

init_hyper_succ_fit_IM_S6 <- c(log(IM_S6$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature", "Precision for Year_Exploited_Mature", "Precision for Depth_Unexploited_Immature", "Precision for Depth_Exploited_Immature", "Precision for Depth_Exploited_Mature"), "mean"]), IM_S6$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"])

init_hyper <- init_hyper_succ_fit_IM_S6

model_setup <- list("description" = description_IM_S8,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = NA, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X,
                                       "Stack" = mystack)

# fit inla model
rm(IM)
IM = inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A = inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE),
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               control.mode      = list(restart=T, theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode = "experimental")

#check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) != 0

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

################
# save inla model & model setup
IM_S8 <- IM

IM_filename <- paste(model_name,"3RWsYear_3RWsDepth_substrate", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

############
IM <- IM_S8
summary(IM)

spatio-temporal model (fix SurveyQuarter to Survey!!!)

# set discrete or continuous temporal components
startyr <- min(model_data$Year)
endyr <- max(model_data$Year)

# for continuous temporal component
y_mesh <- inla.mesh.1d(loc = seq(startyr, endyr, length.out = length(startyr:endyr)/4)); y_mesh$n # here we control the resolution of the temporal mesh; set seq(startyr,endyr, length.out = length(startyr:endyr)/1) to get nods for each year,  change 1 to 2 for each 2nd year, to 5 for every 5th year

# define mesh
mesh_name <- deparse(substitute(mesh_nonconvex_landboundary_outerbuffer_UTM_km))
mesh <- mesh_nonconvex_landboundary_outerbuffer_UTM_km

# define projection matrix
A <- inla.spde.make.A(mesh = mesh,
                      loc = obsloc,
                      group = model_data$Year_Quarter, # use year_quarter resolution
                      group.mesh = y_mesh) # include for continuous temporal component
dim(A)

# define spde
spde <- inla.spde2.matern(mesh)

# define index vector for spde
w.st <- inla.spde.make.index('w', n.spde = spde$n.spde, n.group = y_mesh$n)

IM_ST9 spatio-temporal model with 3 RWs for Year

model_name <- "IM_ST9_no_priors"
description_IM_ST9 <- "spatio-temporal INLA with 3 RWs for Year Exploit_Maturity specific"

# design matrix
X <- model.matrix( ~ fSurvey * Exploit_Maturity, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//:]","",names(X))) # remove special characters from levels of EUNIS2019D_substrate
head(X)

# create stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature), 
                                       1, 
                                       1,
                                       A),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                                      )
                       )   

# define formula
formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept ",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde, group = w.group, control.group = list(model = "ar1"))'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

init_hyper_succ_fit_IM_S5 <-c(log(IM_S5$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature", "Precision for Year_Exploited_Mature"), "mean"]),
                              IM_S5$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"],
                              2.38) # here we provide a best guess from previous analyses for temporary correlation starting value

init_hyper <- init_hyper_succ_fit_IM_S5; init_hyper

model_setup <- list("description" = description_IM_ST9,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = y_mesh, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X,
                                       "Stack" = mystack)

# fit INLA model once
IM = inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A=inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE,  config = TRUE), #openmp.strategy = "pardiso",
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               control.mode      = list(restart=T,  theta = init_hyper), # set initial values for hyperparameters
               #control.family    = list(hyper=list(theta=list(prior= "normal", param= c(log(0.236),70)))), 
               verbose           = TRUE,
               inla.mode         = "experimental")

# alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

##################
IM_ST9 <- IM

IM_filename <- paste(model_name, "3RWsYear_9tmpknots", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM_ST9, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))


#################
summary(IM)

IM$dic$dic

plots_RandomEffects <- list()
sumrandom <- IM$summary.random
for(i in 1:length(sumrandom)){
  covariate <- names(sumrandom )[i]
  plots_RandomEffects[[covariate]] <- sumrandom[[covariate]] %>%
    ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(sides="b") +
    theme_bw() +
    ggtitle(covariate)
}


patchwork::wrap_plots(plots_RandomEffects[1:4]) + patchwork::plot_layout(nrow = 3, ncol =3)

IM_ST10 spatio-temporal 3 rWs year-Exploit_Maturity + 3 depth RW for Exploit_Maturity specific

model_name <- "IM_ST10_no_priors"
description_IM_ST10 <- "spatio-temporal INLA with 3 RWs for Year Exploit_Maturity specific & 3 Depth RW Exploit_Maturity specific"

# design matrix
X <- model.matrix( ~ fSurvey * Exploit_Maturity, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//:]","",names(X))) # remove special characters from levels of

# create stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                       w = w.st
                        )
)

formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept ",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde, group = w.group, control.group = list(model = "ar1"))'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

# inits coming from IM_ST9 and for Depth from IM_S6
init_hyper <- c(log(IM_ST9$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature", "Precision for Year_Exploited_Mature"), "mean"]),
                log(IM_S6$summary.hyperpar[c("Precision for Depth_Unexploited_Immature", "Precision for Depth_Exploited_Immature", "Precision for Depth_Exploited_Mature"), "mean"]),
                IM_ST9$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"], exp(IM_ST9$summary.hyperpar["GroupRho for w", "mean"]))

# store model setup details in a list
model_setup <- list("description" = description_IM_ST10,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh, 
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,                   
                                       "y_mesh" = y_mesh, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X, 
                                       "Stack" = mystack)


# fit INLA model once
rm(IM)
IM <- inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A=inla.stack.A(mystack), compute=TRUE),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #openmp.strategy = "pardiso",
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               control.mode      = list(restart=T,  theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode         = "experimental")

i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}

##################
# save model results and model setup
IM_ST10 <- IM

IM_filename <- paste(model_name, "3RWsYear_9tmpknots", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path,  paste("model_setup", IM_filename, sep="_")))

summary(IM)

IM$dic$dic

##################
IM <- IM_ST10

summary(IM)

IM$dic$dic
 
# set y scale for Random walks
y_scale_year <- range(cbind(IM$summary.random$Year_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_depth <- range(cbind(IM$summary.random$Depth_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_w <- range(cbind(IM$summary.random$w[,c("0.025quant", "0.975quant")], IM$summary.random$w[,c("0.025quant", "0.975quant")]))

plots_RandomEffects <- list()
sumrandom <- IM$summary.random
for(i in 1:length(sumrandom)){
  covariate <- names(sumrandom )[i]
  if(str_detect(covariate, pattern = "Year")){y_scale <- y_scale_year}
  if(str_detect(covariate, pattern = "Depth")){y_scale <- y_scale_depth}
  if(str_detect(covariate, pattern = "w")){y_scale <- y_scale_w}


  plots_RandomEffects[[covariate]] <- sumrandom[[covariate]] %>%
    ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(sides="b") +
    ylim(y_scale) +
    theme_bw() +
    ggtitle(covariate)
    # fix scale across RW Mature and Immature for year and depth seperately


}

patchwork::wrap_plots(plots_RandomEffects[1:6]) + patchwork::plot_layout(nrow = 3, ncol =3)

IM_ST11 spatio-temporal 3 RWs year-Exploit_Maturity + 3 depth RW for Exploit_Maturity specific + EUNIS2019C

model_name <- "IM_ST11_no_priors"
description_IM_ST11 <- "spatio-temporal INLA with 3 RWs for Year Exploit_Maturity specific & 3 Depth RW Exploit_Maturity specific + EUNIS2019C"

# design matrix
#X <- model.matrix( ~ fSurvey + Substrate, data=model_data)
X <- model.matrix( ~ fSurvey * Exploit_Maturity + fEUNIS2019C, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//_ | //:]","",names(X)))
head(X)

# create stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                      w = w.st
                        )
)

# define formula
## no priors
formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde, group = w.group, control.group = list(model = "ar1"))'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

# set starting values form IM_ST10
init_hyper <- c(log(IM_ST10$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature","Precision for Year_Exploited_Mature", "Precision for Depth_Unexploited_Immature", "Precision for Depth_Exploited_Immature","Precision for Depth_Exploited_Mature"), "mean"]), IM_ST10$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"], exp(IM_ST10$summary.hyperpar["GroupRho for w", "mean"]))

# store model setup details in a list
model_setup <- list("description" = description_IM_ST11,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = y_mesh, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X, 
                                       "Stack" = mystack)

#model_setup_list[[model_name]] <- model_setup_IM9
rm(IM)
IM <- inla(formula,
             data              = inla.stack.data(mystack),
             family            = "nbinomial",
             control.predictor = list(A=inla.stack.A(mystack), compute=TRUE),
             control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE,  config = TRUE), #openmp.strategy = "pardiso",
             control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
             control.mode      = list(restart=T,  theta = init_hyper), # set initial values for hyperparameters
             #control.family    = list(hyper=list(theta = list(prior = "normal", param= c(log(0.236), 70)))),
             verbose           = TRUE,
             inla.mode         = "experimental")
 
 # alternative while loop of fitting inla models repetitively until a good fit is found
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){
  (i = i+1)
  try(
    IM <- inla.rerun(IM)
    ) # try to continue loop with new attempt after an error encountered
}

##################
# save model results and model setup
IM_ST11 <- IM

IM_filename <- paste(model_name, "3RWsYear_9tmpknots", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

##################
# model results
##################
summary(IM)

IM$dic$dic
 
# set y scale for Random walks
y_scale_year <- range(cbind(IM$summary.random$Year_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_depth <- range(cbind(IM$summary.random$Depth_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_w <- range(cbind(IM$summary.random$w[,c("0.025quant", "0.975quant")], IM$summary.random$w[,c("0.025quant", "0.975quant")]))

plots_RandomEffects <- list()
sumrandom <- IM$summary.random
for(i in 1:length(sumrandom)){
  covariate <- names(sumrandom )[i]
  if(str_detect(covariate, pattern = "Year")){y_scale <- y_scale_year}
  if(str_detect(covariate, pattern = "Depth")){y_scale <- y_scale_depth}
  if(str_detect(covariate, pattern = "w")){y_scale <- y_scale_w}


  plots_RandomEffects[[covariate]] <- sumrandom[[covariate]] %>%
    ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(sides="b") +
    ylim(y_scale) +
    theme_bw() +
    ggtitle(covariate)
    # fix scale across RW Mature and Immature for year and depth seperately


}

patchwork::wrap_plots(plots_RandomEffects[1:6]) + patchwork::plot_layout(nrow = 3, ncol =3)

IM_ST12 spatio-temporal 3 RWs year-Exploit_Maturity + 3 depth RW for Exploit_Maturity specific + Substrate

model_name <- "IM_ST12_no_priors"
description_IM_ST12 <- "spatio-temporal INLA with 3 RWs for Year Exploit_Maturity specific & 3 Depth RW Exploit_Maturity specific + Substrate"

# create design matrices for covariates
X <- model.matrix( ~ fSurvey * Exploit_Maturity + fSubstrate, data=model_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X)<- c(gsub("[// |//:]|[[]|[]]", "",names(X))) # due to " " (space) & [ & ] in Substrate factor levels (spaces in factor levels)
head(X)

# create stack
mystack <-  inla.stack(data    = list(NoPerHaul = NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature), 
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                        effects = list(Intercept       = rep(1, length(NoPerHaul)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature, 
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = logSweptAreakm2,
                                      w = w.st
                        )
)

# define formula
## no priors
formula = as.formula(str_c("NoPerHaul",
                      str_c( "-1 + Intercept",
                        paste(names(X)[!names(X) == "logSweptAreakm2"], collapse = " + "), #paste covariates
                        paste0('f(Year_Unexploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Immature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Year_Exploited_Mature,  model = "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Unexploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Immature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(Depth_Exploited_Mature, model= "rw1", scale.model = TRUE)'),
                        paste0('f(w, model = spde, group = w.group, control.group = list(model = "ar1"))'), 
                        paste0('offset(logSweptAreakm2)'),
                        sep = " + "), sep = " ~ ")); formula

# set starting values form IM_ST10
init_hyper <- c(log(IM_ST10$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", "Precision for Year_Unexploited_Immature", "Precision for Year_Exploited_Immature","Precision for Year_Exploited_Mature", "Precision for Depth_Unexploited_Immature", "Precision for Depth_Exploited_Immature","Precision for Depth_Exploited_Mature"), "mean"]), IM_ST10$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"], exp(IM_ST10$summary.hyperpar["GroupRho for w", "mean"]))

# store model setup details in a list
model_setup <- list("description" = description_IM_ST12,
                                       "data" = model_data, 
                                       "Formula" = formula, 
                                       "init_hyper" = init_hyper, 
                                       "mesh_name" = mesh_name, 
                                       "mesh" = mesh,
                                       "interior_mesh_poly" = nonconvex_landboundary_poly,
                                       "y_mesh" = y_mesh, 
                                       "A" = A, 
                                       "spde_name" = NA, 
                                       "spde" = spde, 
                                       "w.st" = w.st, 
                                       "X" = X, 
                                       "Stack" = mystack)

#model_setup_list[[model_name]] <- model_setup_IM9
rm(IM)
IM <- inla(formula,
             data              = inla.stack.data(mystack),
             family            = "nbinomial",
             control.predictor = list(A=inla.stack.A(mystack), compute=TRUE),
             control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE,  config = TRUE), #openmp.strategy = "pardiso",
             control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
             control.mode      = list(restart=T,  theta = init_hyper), # set initial values for hyperparameters
             #control.family    = list(hyper=list(theta = list(prior = "normal", param= c(log(0.236), 70)))),
             verbose           = TRUE,
             inla.mode         = "experimental")
 
# check if hessian was successful/happy by checking if the code below evaluates to TRUE (if the hessian would be unhappy then it would equal 0)
i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){
  (i = i+1)
  try(
    IM <- inla.rerun(IM)
    ) # try to continue loop with new attempt after an error encountered
}

##################
IM_ST12 <- IM

# save model results and model setup
IM_filename <- paste(model_name, "3RWsYear_9tmpknots", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path, paste("model_setup", IM_filename, sep="_")))

##################
# model results
##################
summary(IM)

IM$dic$dic


# set y scale for Random walks
y_scale_year <- range(cbind(IM$summary.random$Year_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_depth <- range(cbind(IM$summary.random$Depth_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_w <- range(cbind(IM$summary.random$w[,c("0.025quant", "0.975quant")], IM$summary.random$w[,c("0.025quant", "0.975quant")]))

plots_RandomEffects <- list()
sumrandom <- IM$summary.random
for(i in 1:length(sumrandom)){
  covariate <- names(sumrandom )[i]
  if(str_detect(covariate, pattern = "Year")){y_scale <- y_scale_year}
  if(str_detect(covariate, pattern = "Depth")){y_scale <- y_scale_depth}
  if(str_detect(covariate, pattern = "w")){y_scale <- y_scale_w}


  plots_RandomEffects[[covariate]] <- sumrandom[[covariate]] %>%
    ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(sides="b") +
    ylim(y_scale) +
    theme_bw() +
    ggtitle(covariate)
    # fix scale across RW Mature and Immature for year and depth seperately


}

patchwork::wrap_plots(plots_RandomEffects[1:6]) + patchwork::plot_layout(nrow = 3, ncol =3)

model comparison

Load all Exploit_Maturity INLA models

# list files in model path
list.files(file.path(model_output_path))

# all INLA model files & model setup files
IM_models <- tibble(IM_file = list.files(file.path(model_output_path))) %>% 
  mutate(model_setup_file = paste0("model_setup_", IM_file)) %>%
  filter(str_starts(IM_file, "IM\\w+\\d+_no_priors")) %>% 
  mutate(model_name = case_when(str_detect(IM_file, "pred") == TRUE ~ paste0(str_extract(IM_file, "IM_[A-Z]+[0-9]+"), "_pred"),
                                TRUE ~ str_extract(IM_file, "IM_[A-Z]+[0-9]+"))) #str_subset("IM\\w+\\d+_no_priors")
IM_models

model_list <- list()
model_setup_list <- list()
for (model_file in IM_models$IM_file){
  model_name <- IM_models %>% filter(IM_file == model_file) %>% pull(model_name)
  model_list[[model_name]] <- readRDS(file.path(model_output_path, model_file))
  model_setup_list[[model_name]] <- readRDS(file.path(model_output_path, filter(IM_models, IM_file == model_file) %>% pull(model_setup_file) ))
}

create model comparison table and determine the best model for model fitting with predictions based on DIC, WIAC and LCPO

# function to calculate LCPO
# calculate LCPO
slcpo <- function(m, na.rm = TRUE) {
  - sum(log(m$cpo$cpo), na.rm = na.rm)
}

models_eval_tbl <- tibble("ModelID" = NA, "DIC" = NA, "WAIC" = NA, "LCPO" = NA)
for(model_name in names(model_list)){
  n <- which(names(model_list) == model_name)
  IM <- model_list[[model_name]]
  dic <- IM$dic$dic
  waic <- IM$waic$waic
  lcpo <- slcpo(IM)
  tmp_tib <- tibble("ModelID" = model_name, "DIC" = dic, "WAIC" = waic, "LCPO" = lcpo)
  if(n == 1){
    models_eval_tbl <- tmp_tib }
  else{
    models_eval_tbl <- models_eval_tbl %>% add_row(tmp_tib)
  }
}

models_eval_tbl
models_eval_flex <- models_eval_tbl %>% flextable::flextable() %>% flextable::width( width = 3, unit = "cm");# models_eval_flex
models_eval_flex %>% flextable::save_as_docx(path = file.path(model_output_path, "model_evaluation_summary.docx"))

fitting best model IM_ST10 with predictions

Set some parms and read data IM_ST10

mesh_name <- "mesh_nonconvex_landboundary_outerbuffer_UTM_km"
mesh_n <- mesh_nonconvex_landboundary_outerbuffer_UTM_km$n
last_dat_yr <- 2023

IM_ST10 <- read_rds(file = file.path(model_output_path, paste("IM_ST10_no_priors_3RWsYear_9tmpknots", mesh_name, mesh_n,"1988",last_dat_yr,".rds", sep = "_")))

prepare prediction dataframe for IM_ST10

#make grid for predictions
rangex <- range(model_data$XkmRounded)
rangey <- range(model_data$YkmRounded)
gridcellsize <- 5 # 5km

## sf object unique grid points
grid <- expand.grid(XkmRounded = seq(rangex[1], rangex[2], gridcellsize),
                    YkmRounded = seq(rangey[1], rangey[2], gridcellsize))

# wld_sf 
wld <- maps::map("world", fill=TRUE, col="transparent",  plot=FALSE,  xlim=c(-5,15), ylim=c(47,62))
wld_sf <- wld %>% st_as_sf() %>% st_transform("+proj=utm +zone=31 +units=km")

# only retain grid cells within inner mesh poly determined by NS_EEC_poly
grid_sf            <- grid %>% st_as_sf(coords = c("XkmRounded", "YkmRounded")) %>% st_set_crs("+proj=utm +zone=31 +datum=WGS84 +units=km +no_defs ")
grid_sf            <- grid_sf[as.data.frame(grid_sf %>% st_intersects(nonconvex_landboundary_poly))$row.id,] # only retain grid cells in convex hull polygon used for mesh
wld_cropped_sf     <- wld_sf %>% st_crop(st_bbox(nonconvex_landboundary_poly))                               # crop wld map to boundary of convHullPoly to speed up in water detection
grid_in_water_indx <- st_disjoint(grid_sf, wld_cropped_sf)                                                   # index of grid cells in water
grid_sf            <- grid_sf[lengths(grid_in_water_indx) == max(lengths(grid_in_water_indx)), ]             # only retain grid cells in water

# create prediction dataframe 
grid_pred <- grid_sf %>% st_coordinates() %>% as_tibble() %>% transmute(XkmRounded =X, YkmRounded = Y) %>%
  group_by(XkmRounded, YkmRounded) %>%
  expand_grid(Year_Unexploited_Immature = sort(unique(model_data$Year)),
              Unexploited_Immature = c(1, 0, 0)
              )

grid_pred <- grid_pred %>%
  mutate(NoPerHaul       = NA,
         logSweptAreakm2 = log(gridcellsize^2),
         Year_Exploited_Immature     = Year_Unexploited_Immature,
         Exploited_Immature = rep(c(0, 1, 0), length.out = nrow(grid_pred)),
         Year_Exploited_Mature     = Year_Unexploited_Immature,
         Exploited_Mature = rep(c(0, 0, 1),  length.out = nrow(grid_pred)),
         Year_Quarter    = Year_Unexploited_Immature,
         fSurvey         = case_when(Unexploited_Immature == 1 ~ "BTS",
                                     Exploited_Immature == 1  ~ "FRCGFS",
                                     Exploited_Mature == 1  ~ "FRCGFS"), # predictions for surveys are based on size class
         prediction = TRUE
         )

# add Exploit_maturity for interaction effect
grid_pred <- grid_pred %>% 
  mutate(Exploit_Maturity = case_when(Unexploited_Immature == 1 ~ "unexploited_immature",
                                      Exploited_Immature == 1 ~ "exploited_immature",
                                      Exploited_Mature == 1 ~ "exploited_mature",
                                      ))

# add grid cell id
grid_pred <- grid_pred %>% group_by(XkmRounded, YkmRounded) %>% mutate(ID_grid_cell = cur_group_id()) %>% relocate(ID_grid_cell)

# add Depth to grid_pred based on bathymetry map used for model
## read bathymetry GNS file
load(file.path(data_path, "bathymetry_GNS_low_res_1_5th.Rdata"))
bathymetry_GNS <- bathymetry_GNS %>% select(Xkm, Ykm, Depth1m) %>% mutate(Depth1m = -1 * Depth1m)

# only retain gridcells with depth used in model fitting
model_data_depth_range <- sort(unique(model_data$Depth1m))
bathymetry_GNS_depth_filtered <- bathymetry_GNS %>% filter(Depth1m %in% model_data_depth_range)

grid_pred <- grid_pred %>% 
              left_join(bathymetry_GNS_depth_filtered, by = c("XkmRounded" = "Xkm", "YkmRounded" = "Ykm")) %>%
              dplyr::rename(Depth1m_Unexploited_Immature = Depth1m) %>% 
              mutate(Depth1m_Exploited_Immature = Depth1m_Unexploited_Immature, Depth1m_Exploited_Mature = Depth1m_Unexploited_Immature) 

grid_pred <- grid_pred %>% filter(!is.na(Depth1m_Unexploited_Immature)) # remove grids with NA depth

# combine original data with the newly made grid prediction points that we need for forecasting/predicting
bothnames <- intersect(names(model_data), names(grid_pred)) ; bothnames
comb_data <- bind_rows(subset(model_data, select = c("HaulCode", bothnames)), subset(grid_pred, select = c(bothnames, "prediction"))) %>% mutate(prediction = replace_na(prediction, FALSE))

#remove some objects to seave mem
rm("grid_in_water_indx")
rm("bathymetry_GNS")

# relevel Exploit_Maturity so to decide Unexploited_Immature as intercept
comb_data <- comb_data %>% 
  mutate(Exploit_Maturity = factor(Exploit_Maturity, levels = c("unexploited_immature", "exploited_immature", "exploited_mature")))

comb_data <- comb_data %>% mutate(Year = Year_Unexploited_Immature)

head(comb_data)
dim(comb_data)

fit IM_ST10 no priors pred

init_hyper come from IM_ST10

model_name <- "IM_ST10_no_priors_3RWsYear_9tmpknots_mesh_nonconvex_landboundary_outerbuffer_UTM_km_1407_1988_2023"
model_setup_IM_ST10 <- read_rds(file.path(model_output_path, paste("model_setup", model_name,".rds", sep="_")))

model_name_pred <- paste(c(strsplit(model_name, "_")[[1]][1:2],"no_priors_pred"), collapse="_")
model_name_pred_short <- paste(c(strsplit(model_name, "_")[[1]][1:2],"pred"), collapse="_")

Year_Unexploited_Immature  <- comb_data$Year_Unexploited_Immature
Year_Exploited_Immature    <- comb_data$Year_Exploited_Immature
Year_Exploited_Mature      <- comb_data$Year_Exploited_Mature
Unexploited_Immature       <- comb_data$Unexploited_Immature
Exploited_Immature         <- comb_data$Exploited_Immature
Exploited_Mature           <- comb_data$Exploited_Mature
Depth_Unexploited_Immature <- comb_data$Depth1m_Unexploited_Immature
Depth_Exploited_Immature   <- comb_data$Depth1m_Exploited_Immature
Depth_Exploited_Mature     <- comb_data$Depth1m_Exploited_Mature
Year_Quarter               <- comb_data$Year_Quarter

# design matrix
X <- model.matrix( ~ fSurvey * Exploit_Maturity, data=comb_data)
X <- as.data.frame(X[,-1]) # remove Intercept
names(X) <- c(gsub("[// |//:]","",names(X))) # remove special characters from levels

# for continuous temporal component
y_mesh <- model_setup_IM_ST10$y_mesh 
  
# define mesh
mesh <- model_setup_IM_ST10$mesh 

# projection matrix A requires new Loc of predictions
# define projection matrix
A <- inla.spde.make.A(mesh = mesh,
                      loc = cbind(comb_data$XkmRounded , comb_data$YkmRounded),
                      group = comb_data$Year_Quarter, # use year_quarter resolution
                      group.mesh = y_mesh) # include for continuous temporal component
dim(A)

# define spde
spde <- inla.spde2.matern(mesh)

# define index vector for spde
w.st <- inla.spde.make.index('w', n.spde = spde$n.spde, n.group = y_mesh$n)

# create stack
mystack <-  inla.stack(data    = list(NoPerHaul = comb_data$NoPerHaul),
                        tag     ="Fit",
                        A       = list(1, 
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature), 
                                       Diagonal(length(Exploited_Immature), Exploited_Immature),
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       Diagonal(length(Unexploited_Immature), Unexploited_Immature),
                                       Diagonal(length(Exploited_Immature), Exploited_Immature),
                                       Diagonal(length(Exploited_Mature), Exploited_Mature),
                                       1, 
                                       1,
                                       A),
                        effects = list(Intercept       = rep(1, nrow(comb_data)),
                                       Year_Unexploited_Immature   = Year_Unexploited_Immature, 
                                       Year_Exploited_Immature = Year_Exploited_Immature,
                                       Year_Exploited_Mature     = Year_Exploited_Mature, 
                                       Depth_Unexploited_Immature  = Depth_Unexploited_Immature,
                                       Depth_Exploited_Immature = Depth_Exploited_Immature,
                                       Depth_Exploited_Mature    = Depth_Exploited_Mature,
                                       X               = X,
                                       logSweptAreakm2 = comb_data$logSweptAreakm2,
                                       w = w.st
                        )
)

formula  <-  model_setup_IM_ST10$Formula

init_hyper <- c(log(IM_ST10$summary.hyperpar[c("size for the nbinomial observations (1/overdispersion)", 
                                              "Precision for Year_Unexploited_Immature", 
                                              "Precision for Year_Exploited_Immature", 
                                              "Precision for Year_Exploited_Mature"), "mean"]),
                log(IM_ST10$summary.hyperpar[c("Precision for Depth_Unexploited_Immature",
                                               "Precision for Depth_Exploited_Immature",
                                               "Precision for Depth_Exploited_Mature"), "mean"]),
                IM_ST10$summary.hyperpar[c("Theta1 for w", "Theta2 for w"), "mean"], exp(IM_ST10$summary.hyperpar["GroupRho for w", "mean"]))

# store model setup details in a list, take previous model setup and change relevant parts
model_setup <- model_setup_IM_ST10
model_setup$description <- "ST10 & predictions"
model_setup$data        <- comb_data 
model_setup$init_hyper  <- init_hyper 
model_setup$A           <- A  
model_setup$spde        <- spde 
model_setup$w.st        <- w.st 
model_setup$X           <- X 
model_setup$Stack       <- mystack

# fit INLA model once
## set links for predictions which needs to be added with link = link in control.predictor(..., link = link)
link = rep(NA, nrow(comb_data))
link[which(is.na(comb_data$NoPerHaul))] = 1
threads <- 8 # too many threads/cores may result in exceeding available memory

IM <- inla(formula,
               data              = inla.stack.data(mystack),
               family            = "nbinomial",
               control.predictor = list(A=inla.stack.A(mystack), compute=TRUE, link = link),
               control.compute   = list(dic = TRUE, waic = TRUE, cpo = TRUE, config = TRUE), #openmp.strategy = "pardiso",
               control.inla      = list(strategy = "adaptive", int.strategy = "eb"),
               control.mode      = list(restart=T,  theta = init_hyper), # set initial values for hyperparameters
               verbose           = TRUE,
               inla.mode         = "experimental",
               num.threads = threads) #inla uses all cores by default, this may result in exceeding available memory, down-scaling num core prevents this

i = 1
while (sum(abs(IM$misc$cov.intern[upper.tri(IM$misc$cov.intern)])) == 0 ){ 
  (i = i+1)
  IM <- inla.rerun(IM)
}
##################
# save model results and model setup
IM_ST10_pred <- IM

IM_filename <- paste(model_name_pred, "3RWsYear_9tmpknots", mesh_name, mesh$n,"1988",last_dat_yr,".rds", sep = "_"); IM_filename
write_rds(IM, file = file.path(model_output_path, IM_filename))
write_rds(model_setup, file = file.path(model_output_path,  paste("model_setup", IM_filename, sep="_")))

##################
summary(IM)

IM$dic$dic

# set y scale for Random walks
y_scale_year <- range(cbind(IM$summary.random$Year_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Year_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_depth <- range(cbind(IM$summary.random$Depth_Unexploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Immature[,c("0.025quant", "0.975quant")], IM$summary.random$Depth_Exploited_Mature[,c("0.025quant", "0.975quant")]))
y_scale_w <- range(cbind(IM$summary.random$w[,c("0.025quant", "0.975quant")], IM$summary.random$w[,c("0.025quant", "0.975quant")]))

plots_RandomEffects <- list()
sumrandom <- IM$summary.random
for(i in 1:length(sumrandom)){
  covariate <- names(sumrandom )[i]
  if(str_detect(covariate, pattern = "Year")){y_scale <- y_scale_year}
  if(str_detect(covariate, pattern = "Depth")){y_scale <- y_scale_depth}
  if(str_detect(covariate, pattern = "w")){y_scale <- y_scale_w}
  
  
  plots_RandomEffects[[covariate]] <- sumrandom[[covariate]] %>%
    ggplot(aes(x = ID, y = mean, ymin = `0.025quant`, ymax = `0.975quant`)) +
    geom_line() +
    geom_ribbon(alpha = 0.2) +
    geom_rug(sides="b") +
    ylim(y_scale) +
    theme_bw() +
    ggtitle(covariate)
    # fix scale across RW Mature and Immature for year and depth seperately
}

patchwork::wrap_plots(plots_RandomEffects[1:6]) + patchwork::plot_layout(nrow = 3, ncol =3)

calculate biomass for prediction model by accounting for gear Exploit_Maturity catchability

IM_ST10_pred Results model predictions

# IM <- read_rds(file = "/home/WUR/staud005/RJC_3Lclass_3RW_INLA_2022/JJ_final_RJC_3Lclass_3RW_INLA/RJC_3Lclass_3RW_INLA/final_paper_repo/output/INLA_models/INLAv23.04.24_nonconvex_landboundary_outbuffer_mesh/IM_ST10_no_priors_pred_3RWsYear_9tmpknots_mesh_nonconvex_landboundary_outerbuffer_UTM_km_1407_1988_2023_.rds")
# model_setup <- read_rds(file = "/home/WUR/staud005/RJC_3Lclass_3RW_INLA_2022/JJ_final_RJC_3Lclass_3RW_INLA/RJC_3Lclass_3RW_INLA/final_paper_repo/output/INLA_models/INLAv23.04.24_nonconvex_landboundary_outbuffer_mesh/model_setup_IM_ST10_no_priors_pred_3RWsYear_9tmpknots_mesh_nonconvex_landboundary_outerbuffer_UTM_km_1407_1988_2023_.rds")
# comb_data <- model_setup$data
# IM_ST10_pred <- IM

IM_pred <- IM_ST10_pred

idx1 <- substr(rownames(IM_pred$summary.fitted.values), 1, 10) == "fitted.APr"

comb_data$predmean <- IM_pred$summary.fitted.values[idx1, "mean"]
comb_data$predsd <- IM_pred$summary.fitted.values[idx1, "sd"]
comb_data$`0.025quant` <- IM_pred$summary.fitted.values[idx1, "0.025quant"]
comb_data$`0.5quant` <- IM_pred$summary.fitted.values[idx1, "0.5quant"]
comb_data$`0.975quant` <- IM_pred$summary.fitted.values[idx1, "0.975quant"]

# Exploit_Maturity
comb_data <- within(comb_data,{ Exploit_Maturity = NA
                  Exploit_Maturity[ Unexploited_Immature == 1] <- "Unexploited_Immature"
                  Exploit_Maturity[ Exploited_Immature == 1] <- "Exploited_Immature"
                   Exploit_Maturity[ Exploited_Mature == 1]  <-  "Exploited_Mature"
})

# identify predicted hauls
pred_data <- comb_data %>% filter(prediction == TRUE)

pred_data <- pred_data %>% mutate(predmean_surface = predmean * exp(logSweptAreakm2),
                                  predsd_surface = predsd * logSweptAreakm2,
                                  predmedian_surface = `0.5quant` * exp(logSweptAreakm2)) # multiply predmean with surface

head(pred_data)
# add gear to pred_data
pred_data <- pred_data %>% mutate(Gear = case_when(fSurvey %in% c("FRCGFS", "NSIBTS") ~ "GOV",
                                                   fSurvey %in% c("BTS") ~ "BEAM"))

# add BEAM efficiency to pred_data
pred_data <- pred_data %>% left_join(exploit_maturity_length_class_stats, by = c("Exploit_Maturity", "Gear"))

# adjust ray predmean_surface & predsd_surface by 2 efficiency multiplier:
# - gear Efficiency & group 3 age class 10 efficiency from Walker et al. 2017
# - and estimate biomass per haul by multiplying predicted catch and sd with mean_weight_maturity, we further estimate
walker_group_3_efficiency_largest <- 0.95

pred_data <- pred_data %>%
  mutate(#BEAM_efficiency = 1,
         predicted_catch_mean = predmean_surface * 1/Efficiency * 1/walker_group_3_efficiency_largest,
         predicted_catch_sd = predsd_surface * 1/Efficiency * 1/walker_group_3_efficiency_largest,
         predicted_catch_median = predmedian_surface * 1/Efficiency * 1/walker_group_3_efficiency_largest,
         predicted_biomass_mean = predicted_catch_mean * mean_weight,
         predicted_biomass_sd = predicted_catch_sd * mean_weight,
         predicted_biomass_mean = predicted_catch_median * mean_weight)

head(pred_data) %>% as.data.frame()

pred_results_path <- file.path(results_path, model_name_pred_short)
dir.create(pred_results_path, recursive = TRUE)

write_rds(pred_data, file = file.path(pred_results_path, paste0(model_name_pred_short, "_data.rds")))