In pbs-assess/gfranges: Groundfish Ranges

knitr::opts_chunk$set(
  fig.width = 11, fig.height = 8.5,
  # fig.path=paste0("figs/", spp, "/"),
  echo = FALSE, warning = FALSE, message = FALSE
)
library(dplyr)
library(ggplot2)
library(gfplot)
library(gfdata)
library(sdmTMB)
library(gfranges)

species <- params$species
region <- params$region
covariates <- params$covariates
covs <- params$covs
knots <- params$knots
priors <- FALSE

paste("region =", region)
paste("model covariates =", covariates)
paste("model label =", covs)
paste("priors =", priors)
paste("knots =", knots)

spp <- gsub(" ", "-", gsub("\\/", "-", tolower(species)))

# folder to hold figs for this species
dir.create(file.path("figs", spp))
dir.create(file.path("data", spp))

if (region == "Both odd year surveys") {
  survey <- c("SYN QCS", "SYN HS")
  model_ssid <- c(1, 3)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "West Coast Vancouver Island") {
  survey <- c("SYN WCVI")
  model_ssid <- c(4)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "West Coast Haida Gwaii") {
  survey <- c("SYN WCHG")
  model_ssid <- c(16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

if (region == "All synoptic surveys") {
  survey <- c("SYN QCS", "SYN HS", "SYN WCVI", "SYN WCHG")
  model_ssid <- c(1, 3, 4, 16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  years <- NULL
}

Build spatiotemporal density model

Load and filter data

biomass <- readRDS(paste0(
  "data/", spp, "/data-by-maturity-", spp, "-", ssid_string, ".rds"
))

# if(spp=="redbanded-rockfish"){ 
# 1 redbanded fish found at too shallow a depth... probably leftover in net?
# slender sole also has an outlier here, so maybe should throw out this event for everyone?
# could mean there's a location or depth error  
biomass <- biomass[biomass$fishing_event_id!=2536031,]
# }

if(spp=="pacific-ocean-perch"){
# fish sizes are impossible given recorded catch weight
biomass <- biomass[biomass$fishing_event_id!=1506954,]
}

if(nrow(biomass)<4000) {stop("Need to recalculate split by maturity!")}

covars <- readRDS("data/event-covariates.rds") %>% select(-geometry)
data <- dplyr::left_join(biomass, covars)

# scale predictors before filtering to ensure mean and SD are global
data <- data %>% mutate(raw_depth = depth, depth = log(raw_depth))
data <- gfranges::scale_predictors(data, # predictors = c(quo(depth)))
  predictors = c(quo(depth), quo(mixed), quo(muddy), quo(sandy), quo(rocky), quo(any_rock))
)
data <- data %>%
  mutate(depth = raw_depth) %>%
  filter(ssid %in% model_ssid) %>%
  filter(year < 2019)

Make mesh

if (region == "All synoptic surveys") {
  # spde <- sdmTMB::make_spde(data$X, data$Y, n_knots = knots)
  spde <- sdmTMB::make_mesh(data, c("X", "Y"), n_knots = knots, type = "kmeans")
}

sdmTMB::plot_spde(spde)

Run climate independent sdmTMB model

if (params$update_model) {
  tictoc::tic()
  if (any(names(data) == "adult_density")) {
    # adult_formula <- as.formula(paste(
    #   "adult_density ~ 0 + as.factor(year)", covariates, ""
    # ))
    # 
    # adult_biomass <- sdmTMB::sdmTMB(
    #   data = data,
    #   adult_formula,
    #   time_varying = ~ 0 + depth_scaled + depth_scaled2, #+ trawled
    #   time = "year", spde = spde,
    #   family = tweedie(link = "log"),
    #   #ar1_fields = TRUE,
    #   include_spatial = TRUE,
    #   # reml = TRUE,
    #   nlminb_loops = 2,
    #   newton_steps = 2,
    #   enable_priors = priors,
    #   control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
    #   silent = FALSE
    # )
    # 
    # saveRDS(adult_biomass,
    #   file = paste0(
    #     "data/", spp,
    #     "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-new.rds"
    #   )
    # )


   try({ imm_formula <- as.formula(paste(
      "imm_density ~ 0 + as.factor(year)", covariates, ""
    ))

    imm_biomass <- sdmTMB::sdmTMB(
      data = data,
      imm_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, #+ trawled
      time = "year", spde = spde,
      family = tweedie(link = "log"),
      # ar1_fields = TRUE, # changed to TRUE on Oct 17 2019
      # reml = TRUE,
      include_spatial = TRUE,
      nlminb_loops = 2,
      newton_steps = 2,
      enable_priors = priors,
      control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
      silent = FALSE
    )

    saveRDS(imm_biomass,
      file = paste0("data/", spp, "/mod-imm-biomass-",
        spp, covs, "-", ssid_string, ".rds"))
   })

  } else {
    dens_formula <- as.formula(paste("density ~ 0 + as.factor(year)", covariates, ""))

    total_biomass <- sdmTMB::sdmTMB(
      data = data,
      dens_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, # + trawled
      time = "year", spde = spde,
      family = tweedie(link = "log"),
      # ar1_fields = TRUE, # all off for final version
      include_spatial = TRUE,
      nlminb_loops = 2, 
      newton_steps = 2,
      enable_priors = priors,
      control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
      silent = FALSE
    )

    saveRDS(total_biomass, file = paste0(
      "data/", spp,
      "/model-total-biomass-", spp, covs, "-", ssid_string, ".rds"
    ))
  }
  tictoc::toc()
}

Make prediction grids for all surveys and years, if not done before.

try(nd_all <- readRDS("data/nd_all_synoptic.rds"))

if (!exists("nd_all")) {

  # choose base year(s) to create grid from
  dummy_year <- c(2005, 2006)

  bath <- readRDS("data/bathymetry-data")
  # dat <- bath$data %>% scale_survey_predictors()
  dat <- bath$data %>% mutate(raw_depth = depth, depth = log(raw_depth))
  dat <- gfranges::scale_predictors(dat, predictors = c(quo(depth)))

  # create grids for each ssid separately
  ssid <- 1
  survey_abbrev <- "SYN QCS"
  nd_1 <- spatiotemporal_grid(dat, ssid, survey_abbrev, dummy_year)


  unique(nd_1$year)

  ssid <- 3
  survey_abbrev <- "SYN HS"
  nd_3 <- spatiotemporal_grid(dat, ssid, survey_abbrev, dummy_year)
  unique(nd_3$year)

  ssid <- 4
  survey_abbrev <- "SYN WCVI"
  nd_4 <- spatiotemporal_grid(dat, ssid, survey_abbrev, dummy_year)
  unique(nd_4$year)

  ssid <- 16
  survey_abbrev <- "SYN WCHG"
  nd_16 <- spatiotemporal_grid(dat, ssid, survey_abbrev, dummy_year)
  unique(nd_16$year)

  nd_all <- rbind(nd_1, nd_3, nd_4, nd_16) %>%
    mutate(year = as.numeric(year))

  newcov <- readRDS("data/new_covariates.rds") %>% dplyr::select(X, Y, trawled, muddy, rocky, mixed, any_rock)
  nd_all <- left_join(nd_all, newcov, by = c("X", "Y"))
  nd_all$trawled[is.na(nd_all$trawled)] <- 0 # correct missing values for small area on edge of QCS grid
  # nd <- nd %>% dplyr::mutate(shallow = ifelse(depth > 35, 0, 1))

  saveRDS(nd_all, file = paste0("data/nd_all_synoptic.rds"))
}

Save predictions for full spatial grid

# fish <- readRDS(paste0("raw/bio-data-", spp, ""))

rm(adult_predictions)
rm(imm_predictions)
rm(predictions)
if(params$update_predictions) {
if (any(names(data) == "adult_density")) {
  # if (length(unique(fish$maturity_code)) > 2) {
  nd_all <- readRDS(paste0("data/nd_all_synoptic.rds"))

  # try({
  #   adult_biomass <- readRDS(paste0(
  #     "data/", spp,
  #     "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
  #   ))
  # 
  #   nd <- nd_all %>%
  #     filter(ssid %in% model_ssid) %>%
  #     filter(year %in% unique(adult_biomass$data$year))
  #   nd <- na.omit(nd)
  #   nd$year <- as.integer(nd$year)
  # 
  # 
  #   adult_predictions <- predict(adult_biomass, newdata = nd)
  #   saveRDS(adult_predictions, file = paste0(
  #     "data/", spp,
  #     "/predictions-", spp, covs, "-", ssid_string, "-mature-biomass.rds"
  #   ))
  # })

  try({
    imm_biomass <- readRDS(paste0(
      "data/", spp,
      "/mod-imm-biomass-", spp, covs, "-", ssid_string, ".rds"
    ))

    nd <- nd_all %>%
      filter(ssid %in% model_ssid) %>%
      filter(year %in% unique(imm_biomass$data$year))
    nd <- na.omit(nd)
    nd$year <- as.integer(nd$year)

    imm_predictions <- predict(imm_biomass, newdata = nd)
    saveRDS(imm_predictions, file = paste0(
      "data/", spp,
      "/predictions-", spp, covs, "-", ssid_string, "-imm-biomass.rds"
    ))
  })
} else {
  total_biomass <- readRDS(paste0(
    "data/", spp,
    "/model-total-biomass-", spp, covs, "-", ssid_string, ".rds"
  ))

  nd_all <- readRDS(paste0("data/nd_all_synoptic.rds"))
  nd <- nd_all %>%
    filter(ssid %in% model_ssid) %>%
    filter(year %in% unique(total_biomass$data$year))
  nd <- na.omit(nd)
  nd$year <- as.integer(nd$year)

  predictions <- predict(total_biomass, newdata = nd)
  saveRDS(predictions, file = paste0(
    "data/", spp,
    "/predictions-", spp, covs, "-", ssid_string, "-total-biomass.rds"
  ))
}
}

if(params$update_model_check) {
try({
  biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-", ssid_string, ".rds"
  ))
})

try({
  biomass <- readRDS(paste0(
    "data/", spp,
    "/model-total-biomass-", spp, covs, "-", ssid_string, ".rds"
  ))
})

point_predictions <- predict(biomass)
point_predictions$residuals <- residuals(biomass)

saveRDS(point_predictions, file = paste0(
  "data/", spp,
  "/check-mod-predictions-", spp, covs, "-", ssid_string, ".rds"
))
}

depth_only_predictions <- readRDS(paste0(
  "data/", spp,
  "/check-mod-predictions-", spp, covs, "-", ssid_string, ".rds"
))

depth_only_predictions <- filter(
  depth_only_predictions,
  # ssid %in% c(1,3)) %>% filter (
  year > 2007
)

try({
  g <- ggplot(depth_only_predictions, aes(adult_density, residuals, colour = adult_density)) +
    geom_point() + scale_x_continuous(trans = "log10") + facet_wrap(~year)
})

try({
  g <- ggplot(depth_only_predictions, aes(density, residuals, colour = density)) +
    geom_point() + scale_x_continuous(trans = "log10") + facet_wrap(~year)
})

g

try({
  g2 <- ggplot(depth_only_predictions, aes((est), residuals, colour = adult_density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
})

try({
  g2 <- ggplot(depth_only_predictions, aes((est), residuals, colour = density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
})

g2

ggplot(depth_only_predictions, aes(X, Y, colour = (residuals))) +
  geom_point() +
  scale_colour_gradient2() +
  scale_x_continuous(trans = "log10") +
  facet_wrap(~year)

fits <- readRDS("data/climate-fits.rds")
check_fits <- left_join(depth_only_predictions,fits)

try({
  g3 <- ggplot(check_fits, aes(t_resids, residuals, colour = density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
  g4 <- ggplot(check_fits, aes(d_resids, residuals, colour = density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
})
try({
  g3 <- ggplot(check_fits, aes(t_resids, residuals, colour = adult_density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
  g4 <- ggplot(check_fits, aes(d_resids, residuals, colour = adult_density)) +
    geom_point() + # scale_x_continuous(trans = 'log10') +
    facet_wrap(~year)
})
g3
g4

Build time-varying depth plots

Model without environmental covariates

max_depth_found <- max(data[data$present == 1, ]$raw_depth, na.rm = TRUE)
rm(depth_model_list)
rm(adult_biomass)
rm(imm_biomass)
rm(total_biomass)
rm(depth_plots)
#max_depth_found <- 800

# try({
#   adult_biomass <- readRDS(paste0(
#     "data/", spp,
#     "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
#   ))
# })
# adult_biomass<-sdmTMB:::update_model(adult_biomass)

try({
  imm_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, covs, "-1n3n4n16.rds"
  ))
  # imm_biomass<-sdmTMB:::update_model(imm_biomass)
  # depth_model_list <- list(adult = adult_biomass, imm = imm_biomass)
})

if (!exists("adult_biomass")) {
  try({
    total_biomass <- readRDS(paste0(
      "data/", spp,
      "/model-total-biomass-", spp, covs, "-", ssid_string, ".rds"
    ))
  })
}


if (!exists("adult_biomass")) adult_biomass <- NULL
if (!exists("imm_biomass")) imm_biomass <- NULL
if (!exists("total_biomass")) total_biomass <- NULL

depth_model_list <- list(adult = adult_biomass, imm = imm_biomass, total = total_biomass)
depth_model_list <- depth_model_list[!sapply(depth_model_list, is.null)]


d <- list()
depth_plots <- list()
for (i in seq_len(length(depth_model_list))) {
  #if (depth_model_list[[i]]$model$convergence == 0) {
    d[[i]] <- time_varying_density(depth_model_list[[i]], predictor = "depth")

    if (length(d[[i]]) == 0) {
      depth_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
    } else {
      d[[i]]$x <- exp(d[[i]]$x)
      depth_plots[[i]] <- plot_mountains(d[[i]], variable_label = "Depth (without environmental variables)", xlimits = c(0, max_depth_found)) +
        ggtitle(paste(species, names(depth_model_list[i])))
    }
  # } else {
  #   depth_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
  # }
}
print(depth_plots)

png(
  file = paste0(
    "figs/", spp,
    "/depth-", spp, covs, "-no-AR1.png"
  ),
  res = 600,
  units = "in",
  width = 8.5,
  height = 6
)
gridExtra::grid.arrange(
  grobs = c(depth_plots),
  nrow = 2,
  top = grid::textGrob(paste(species, "(", covs, ")"))
)
dev.off()

depth_model_list

write_grads <- function(x, macro, ...) {
  paste0("\\newcommand{\\", macro, "}{", x, "}") %>%
    readr::write_lines("model-grads.tex", append = TRUE)
}

write_grads(signif(max(depth_model_list$adult$gradients), digits = 8), paste0(spp, "-mat"))
write_grads(signif(max(depth_model_list$imm$gradients), digits = 8), paste0(spp, "-imm-redo"))
write_grads(signif(max(depth_model_list$total$gradients), digits = 8), paste0(spp, "-total-redo"))

try({
  if(max(depth_model_list$adult$gradients)>0.001) {
    adult_biomass <- sdmTMB:::run_extra_optimization(depth_model_list$adult, nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(adult_biomass$gradients), digits = 8), paste0(spp, "-mat2"))
  }
})

try({
  if(max(depth_model_list$imm$gradients)>0.001) {
    imm_biomass <- sdmTMB::run_extra_optimization(depth_model_list$imm, nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(imm_biomass$gradients), digits = 8), paste0(spp, "-imm3"))
  }
})

try({
  if(max(depth_model_list$total$gradients)>0.001) {
    total_biomass <- sdmTMB::run_extra_optimization(depth_model_list$total, nlminb_loops = 2, newton_steps = 2)
    write_grads(signif(max(total_biomass$gradients), digits = 8), paste0(spp, "-total3"))
  }
})

 if (any(names(data) == "adult_density")) {

    # saveRDS(adult_biomass,
    #   file = paste0(
    #     "data/", spp,
    #     "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-prior-", priors, ".rds"
    #   )
    # )

    saveRDS(imm_biomass,
      file = paste0("data/", spp, "/mod-imm-biomass-", 
        spp, covs, "-", ssid_string, ".rds"))

  } else {
    saveRDS(total_biomass, file = paste0(
      "data/", spp,
      "/model-total-biomass-", spp, covs, "-", ssid_string, ".rds"
    ))
  }