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.
rm(list = ls()); gc()
getwd()library(tidyverse)
library(INLA)
library(inlabru)
library(sf)
library(ggplot2)
library(rnaturalearth)Timestamp when saving dowloaded exchange data.
dl_date <- "14_06_2024"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)
}# 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)# 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')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 kmLoc <- 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 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() 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_Maturemodel_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)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)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)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)# 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)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)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)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)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)# 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)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)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)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)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)# 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) ))
}# 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"))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 = "_")))#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)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)# 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")))