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)

# options(scipen = 999)
theme_set(
    gfplot::theme_pbs(base_size = 14) 
)

species <- params$species
region <- params$region
covariates <- params$covariates
covs <- params$covs
knots <- params$knots
priors <- FALSE
update_model <- params$update_model
update_model_check <- params$update_model_check
update_predictions <- params$update_predictions
update_index <- params$update_index

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

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, "-1n3n4n16.rds"
))

  positive_sets <- filter(biomass, density != 0)
  prop_pos <- round(nrow(positive_sets)/nrow(biomass), digits = 3)
      if (prop_pos < 0.1) {
        # update_model <- T
        # update_model_check <- T
        # update_predictions <- T
        # update_index <- T
        knots <- 150
        }
    if (prop_pos < 0.05) {
        # update_model <- F
        # update_model_check <- F
        # update_predictions <- F
        # update_index <- T
        knots <- 150
        }

# View(biomass)
if (nrow(biomass)<2000) {
  # update_model <- T
  # update_model_check <- T
  # update_predictions <- T
  # update_index <- T
  knots <- 100
  if(prop_pos < 0.05) { knots <- 50 }
  }
if(nrow(biomass)<1000) {stop("Need to recalculate split by maturity!")}

# covars <- readRDS("data/event-covariates.rds")
# data <- dplyr::left_join(biomass, covars)
data <- biomass
# 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))
)
data <- data %>%
  mutate(depth = raw_depth) %>% 
  filter(ssid %in% model_ssid) %>%
  filter(year >= params$start_year)

rm(adult_biomass)
rm(imm_biomass)
rm(total_biomass)

Make mesh

if (region == "Both odd year surveys") {
  spde <- sdmTMB::make_spde(data$X, data$Y, n_knots = knots)
}

if (region == "West Coast Vancouver Island") {
  spde <- sdmTMB::make_spde(data$X, data$Y, n_knots = 200)
}

if (region == "West Coast Haida Gwaii") {
  spde <- sdmTMB::make_spde(data$X, data$Y, n_knots = 200)
}

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

sdmTMB::plot_spde(spde)

paste("knots =", knots)

Run climate independent sdmTMB model

if (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 = params$AR1,
      include_spatial = params$fixed_spatial,
      reml = TRUE,
      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, "-ar1-", params$AR1, "-reml.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 = params$AR1, # changed to TRUE on Oct 17 2019
      reml = TRUE,
      include_spatial = params$fixed_spatial,
      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, "-ar1-", params$AR1, "-reml.rds"))
   })

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

    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 = params$AR1,
      reml = TRUE,
      include_spatial = params$fixed_spatial,
      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, "-ar1-", params$AR1, "-reml.rds"
    ))
  }
  tictoc::toc()
}

Save predictions for spatial grid

rm(ad_predictions)
rm(im_predictions)
rm(predicted)

if(update_predictions) {

  # nd_all <- readRDS(paste0("data/nd_just_depth.rds"))
  # nd_all <- readRDS(paste0("data/nd_whole_coast_index.rds"))
  nd_all <- readRDS(paste0("data/nd_odd.rds"))

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

try({

if (!exists("adult_biomass")) {    
  adult_biomass <- readRDS(paste0("data/", spp, 
    "/mod-mat-biomass-", spp, covs, "-", ssid_string, 
    "-ar1-", params$AR1, "-reml.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)


    ad_predictions <- predict(adult_biomass, newdata = nd, return_tmb_object = TRUE)
    # saveRDS(ad_predictions, file = paste0("data/", spp,
    #   "/predictions-", spp, covs, "-", ssid_string, 
    #   "-mature-biomass-ar1-", params$AR1, "-reml.rds"))
    saveRDS(ad_predictions,
      file = paste0("data/", spp,
        "/mod-mat-biomass-", spp, covs, "-", ssid_string, "-ar1-", params$AR1, "-reml.rds"))
})

  try({

    if (!exists("imm_biomass")) {
    imm_biomass <- readRDS(paste0("data/", spp,
      "/mod-imm-biomass-", spp, covs, "-", ssid_string, 
      "-ar1-", params$AR1, "-reml.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)
    im_predictions <- predict(imm_biomass, newdata = nd, return_tmb_object = TRUE)

    # saveRDS(im_predictions, file = paste0("data/", spp,
    #   "/predictions-", spp, covs, "-", ssid_string, "-imm-biomass-ar1-", params$AR1, "-reml.rds"
    # ))

     saveRDS(im_predictions, file = paste0("data/", spp, "/mod-imm-biomass-", 
        spp, covs, "-", ssid_string, "-ar1-", params$AR1, "-reml.rds"))

  })

} else {

  if (!exists("total_biomass")) {
  total_biomass <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, covs, "-", ssid_string, 
    "-ar1-", params$AR1, "-reml.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)

  predicted <- predict(total_biomass, newdata = nd, return_tmb_object = TRUE)

  # saveRDS(predicted, file = paste0("data/", spp,
  #   "/predictions-", spp, covs, "-", ssid_string, 
  #   "-total-biomass-ar1-", params$AR1,"-reml.rds"))

  saveRDS(predicted, file = paste0("data/", spp,
    "/model-total-biomass-", spp, covs, "-", ssid_string, 
    "-ar1-", params$AR1, "-reml.rds"))

}
}

Check residuals

if(update_model_check) {
try({

if (!exists("ad_predictions")) {    
  ad_predictions <- readRDS(paste0("data/", spp, 
    "/mod-mat-biomass-", spp, covs, "-", ssid_string, 
    "-ar1-", params$AR1, "-reml.rds"))
}
  point_predictions <- predict(ad_predictions$fit_obj)
  point_predictions$residuals <- residuals(ad_predictions$fit_obj)
})

try({

if (!exists("total_biomass")) {    
  predicted <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, covs, "-", ssid_string, 
    "-ar1-", params$AR1, "-reml.rds"
  ))
}
  point_predictions <- predict(predicted$fit_obj)
  point_predictions$residuals <- residuals(predicted$fit_obj)
})

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

if(update_model_check) {
check_predictions <- readRDS(paste0(
  "data/", spp,
  "/check-mod-predictions-", spp, covs, "-", ssid_string, "-ar1-", params$AR1, "-reml.rds"
))


check_predictions <- filter(
  check_predictions,
  # ssid %in% c(1,3)) %>% filter (
  year > 2004
)

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

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

g
}

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

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

try ({g2})

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

Load predictions

Transform estimates to kg/ha

rm(adult_predictions)
rm(imm_predictions)
rm(predictions)

if (!exists("ad_predictions")) {
  try({
    ad_predictions <- readRDS(paste0("data/", spp, 
      "/mod-mat-biomass-", spp, covs, "-", ssid_string, 
      "-ar1-", params$AR1, "-reml.rds"))
  })

  if (!exists("ad_predictions")) {
     try({ ad_predictions <- readRDS(paste0("data/", spp,
      "/mod-mat-biomass-", spp, "-no-covs-", ssid_string,
      "-ar1-", params$AR1, "-reml.rds")) })
  }

  try({
    im_predictions <- readRDS(paste0("data/", spp, 
      "/mod-imm-biomass-", spp, covs, "-", ssid_string, 
      "-ar1-", params$AR1, "-reml.rds"))
  })
    if (!exists("im_predictions")) {
     try({ im_predictions <- readRDS(paste0("data/", spp,
      "/mod-imm-biomass-", spp, "-no-covs-", ssid_string,
      "-ar1-", params$AR1, "-reml.rds")) })
  }


  if (!exists("ad_predictions")) {
    try({
      predicted <-  readRDS(paste0("data/", spp,
        "/model-total-biomass-", spp, covs, "-", ssid_string, 
        "-ar1-", params$AR1, "-reml.rds"))
    })
    if (!exists("predicted")) {
     try({
       predicted <- readRDS(paste0("data/", spp,
      "/model-total-biomass-", spp, "-no-covs-", ssid_string,
      "-ar1-", params$AR1, "-reml.rds"))
       })
  }
  }
}

if (exists("ad_predictions")) {
  # to convert from kg/m2 to kg/hectare multiply by 10000
  adult_predictions <- ad_predictions$data
  adult_predictions$est_exp <- exp(adult_predictions$est) * 10000


  try({
    imm_predictions <- im_predictions$data
    imm_predictions$est_exp <- exp(imm_predictions$est) * 10000

    adult_predictions$total_bio <- imm_predictions$est_exp +   
        adult_predictions$est_exp
    imm_predictions$prop_imm <- imm_predictions$est_exp /
        (adult_predictions$est_exp + imm_predictions$est_exp)

  plot_imm <- filter(imm_predictions, year>2004)
  max_imm <-  signif(max(plot_imm$est_exp), digits = 2)
  max_raster_imm <- quantile(plot_imm$est_exp, 0.999)


  saveRDS(imm_predictions, file = paste0("data/", spp,
    "/sopo-predictions-", spp, covs, "-imm-biomass-ar1-", 
    params$AR1, "-reml.rds"))
  })

  saveRDS(adult_predictions, file = paste0("data/", spp,
    "/sopo-predictions-", spp, covs, "-mat-biomass-ar1-", 
    params$AR1, "-reml.rds"))

  plot_ad <- filter(adult_predictions, year>2004)

  max_raster_ad <- quantile(plot_ad$est_exp, 0.999)
  max_adult <- signif(max(plot_ad$est_exp), digits = 2)

  # model_ssid <- unique(adult_predictions$ssid)
} else {

  predictions <- predicted$data
  predictions$est_exp <- exp(predictions$est) * 10000
  predictions$total_bio <- exp(predictions$est) * 10000
  max_raster <- quantile(predictions$total_bio, .99)
  max_bio <- signif(quantile(predictions$total_bio, .999), digits = 2)

  saveRDS(predictions, file = paste0(
    "data/", spp,
    "/sopo-predictions-", spp, covs, "-total-biomass-ar1-", params$AR1, "-reml.rds"
  ))

}

Model summary

 try({ ad_predictions$fit_obj})
 try({ im_predictions$fit_obj})
 try({ predicted$fit_obj})

Plot biomass predictions

legend_coords <- c(0.9, 0.17) #"none"


if (exists("adult_predictions")) {

 p_adult_all <- adult_predictions %>% filter(!((year<2005) & (ssid == 3)))%>%
   mutate(x = X, y = Y, X = 2 * round(X/2), Y = 2 * round(Y/2)) %>%  
   plot_facet_map("est_exp", 
      raster_limits = c(0, max_raster_ad),
      legend_position =  legend_coords, 
      transform_col = fourth_root_power
    ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " mature biomass \n(max = ", max_adult, " kg/ha)"))

 print(p_adult_all)

try({
 p_imm_all <- imm_predictions %>% filter(!((year<2005) & (ssid == 3)))%>%
    mutate(x = X, y = Y, X = 2 * round(X/2), Y = 2 * round(Y/2)) %>%
    plot_facet_map("est_exp",
      raster_limits = c(0, max_raster_imm),
      legend_position =  legend_coords, 
      transform_col = fourth_root_power
    ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " immature biomass \n(max = ", max_imm, " kg/ha)"))

print(p_imm_all)
   })

} else {

p_adult_all <- predictions  %>% filter(!((year<2005) & (ssid == 3))) %>% #filter(year>2003) %>%
  mutate(x = X, y = Y, X = 2 * round(X/2), Y = 2 * round(Y/2)) %>% 
  plot_facet_map("total_bio",
    raster_limits = c(0, max_raster),
    legend_position =  legend_coords, 
    transform_col = fourth_root_power
  ) +
    labs(fill = "kg/ha") +
    ggtitle(paste0("", species, " total biomass \n(max = ", max_bio, " kg/ha)"))

  print(p_adult_all)
}

Calculate index of total biomass for full survey grid

rm(ind)
rm(ind_imm)
library(data.table)
if(update_index) {

if (exists("ad_predictions")) {

  ind <- get_index(ad_predictions, bias_correct = FALSE)
  ind$species <- species
  ind$maturity <- "mature"
  ind$ssid <- ssid_string
  ind$covs <- covs
  ind <- as.data.frame(ind)
  write.csv(ind, file = paste0("data/_indices/sopo-index-", 
    spp,"-", ssid_string, "-", covs, "-reml.csv"), row.names = FALSE)

  try ({
    ind_imm <- get_index(im_predictions, bias_correct = FALSE)
    ind_imm$species <- species
    ind_imm$maturity <- "immature"
    ind_imm$ssid <- ssid_string
    ind_imm$covs <- covs
    ind_imm <- as.data.frame(ind_imm)
    write.csv(ind_imm, file = paste0("data/_indices/sopo-index-imm-", 
      spp,"-", ssid_string, "-", covs, "-reml.csv"), row.names = FALSE)
  })

} else {
  ind <- get_index(predicted, bias_correct = FALSE)
  ind$species <- species
  ind$maturity <- "all"
  ind$ssid <- ssid_string
  ind$covs <- covs

  ind <- as.data.frame(ind)

  write.csv(ind, file = paste0("data/_indices/sopo-index-", 
    spp,"-", ssid_string, "-", covs, "-reml.csv"), row.names = FALSE)
}
try({ glimpse(ind_imm) })
try({ glimpse(ind) })
}

# scale <- 2 * 2 / 1000 # if density was in kg/km2: 2 x 2 km grid and converted from kg to tonnes 
scale <- 1000000 * 4 /1000 # if density was in kg/m2: m2 to km2 * 4 km2 grid size and kg to tonnes
ad_col <- "deepskyblue4"

ind <- read.csv(paste0("data/_indices/sopo-index-", 
  spp,"-", ssid_string, "-", covs, "-reml.csv"))

ind %>% filter(year > 2004) %>%
ggplot(aes(year, est*scale)) + geom_line(col=ad_col) +
  geom_ribbon(aes(ymin = lwr*scale, ymax = upr*scale), fill=ad_col, alpha = 0.4) +
  xlab('Year') +  
  ylab('Mature biomass estimate (metric tonnes)') +
  gfplot::theme_pbs() +
    ggtitle(paste0("", species, " "))

Add immature trend if available

try({
  ind_imm <- read.csv(paste0("data/_indices/sopo-index-imm-", 
    spp,"-", ssid_string, "-", covs, "-reml.csv")) %>% filter(year > 2004)
})

if (exists("ind_imm")) {
ratio <- mean(ind$est)/max(ind_imm$est)
imm_col <- "orangered"

ind %>% filter(year > 2004) %>%
ggplot(aes(year, est*scale)) + geom_line(col = ad_col, ) +
  geom_ribbon(aes(ymin = lwr*scale, ymax = upr*scale), fill = ad_col, alpha = 0.4) +
# adding the relative humidity data, transformed to match roughly the range of the temperature
  geom_ribbon(aes(ymin = ind_imm$lwr*scale*ratio, ymax = ind_imm$upr*scale*ratio), 
    fill = "orangered", alpha = 0.4) + 
  geom_line(aes(ind_imm$year, ind_imm$est*scale*ratio), col = "orangered") +
# now adding the secondary axis, following the example in the help file ?scale_y_continuous and, reverting the above transformation
  scale_y_continuous(sec.axis = sec_axis(~./ratio, name = "Immature biomass estimate (metric tonnes)")) +
  xlab('Year') + 
  ylab('Mature biomass estimate (metric tonnes)')+
  gfplot::theme_pbs(base_size = 16)
}