In pbs-assess/gfranges: Groundfish Ranges

Choose one value for each parameter:

Species

# species <- params$species

# # Species run so far...
# species <- "Arrowtooth Flounder"
# species <- "Petrale Sole"
# species <- "English Sole"
# species <- "Dover Sole" # maturity data error for WCHG
# species <- "Southern Rock Sole"
# species <- "Flathead Sole"



# species <- "Silvergray Rockfish"
# species <- "Quillback Rockfish" # summer-birthing
# species <- "Yelloweye Rockfish" # summer-birthing, overfished, 
# species <- "Bocaccio" # winter-birthing, overfished
# species <- "Sharpchin Rockfish"
# species <- "Splitnose Rockfish"

# species <- "Rougheye/Blackspotted Rockfish Complex"
# species <- "Redbanded Rockfish"
# species <- "Greenstriped Rockfish"
# # species <- "Copper Rockfish" # small sample
# species <- "Darkblotched Rockfish"
# # species <- "Shortbelly Rockfish" # small sample



# species <- "Pacific Ocean Perch" # schooling
# species <- "Widow Rockfish" # schooling
# species <- "Yellowtail Rockfish" # schooling
# species <- "Canary Rockfish" # schooling, winter-birthing

# species <- "Longspine Thornyhead"
# species <- "Shortspine Thornyhead"


# species <- "Walleye Pollock"
# species <- "Pacific Cod"
# species <- "Sablefish"
# species <- "Lingcod"
# species <- "Pacific Hake"

# species <- "North Pacific Spiny Dogfish" # note: using all data for maturity thresholds
# species <- "Longnose Skate"
# species <- "Big Skate"
# species <- "Spotted Ratfish"
# species <- "Sandpaper Skate"
# species <- "Brown Cat Shark"

Region

 region <- "All synoptic surveys"
  # region <- "Both odd year surveys"
# region <- "West Coast Vancouver Island"
 #  region <- "West Coast Haida Gwaii"

Choose model details

# priors <- FALSE
 priors <- TRUE
# covariates <- "+muddy+any_rock"
 covariates <- "+as.factor(ssid)"
# covariates <- "+trawled+as.factor(ssid)"
# covariates <- "+muddy+mixed+rocky"
# covariates <- "+mixed+rocky"
# covariates <- "+trawled+muddy+rocky+mixed"
# covariates <- "+trawled+mixed+rocky"
# covariates <- "+trawled+mixed"
# covs <- gsub("\\+", "-", covariates)
# covs <- "-both-tv-depth-ssid"
# covs <- "-both-tv-depth"
covs <- "-log-do"
covs <- "-log-do-ssid"

Run all subsequent code...

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

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

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)

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
}

Combine biomass and sensor data

  biomass <- readRDS(paste0(
    "data/", spp, "/data-by-maturity-", spp, "-", ssid_string, ".rds"
    ))
  covars <- readRDS("data/event-covariates.rds")
  data <- dplyr::left_join(biomass, covars) %>% mutate(depth_akima = depth) %>% select(-depth)

  CTD <- readRDS("../tmb-sensor-explore/data/all-sensor-data-processed.rds")
  data <- dplyr::left_join(data, CTD) %>% 
    mutate(exclude = if_else(do_mlpl>8, 1, 0)) %>% 
    #mutate(exclude = if_else(do_mlpl>6, 1, 0)) %>% #try excluding high values?
    filter(exclude != 1) %>% # added on july 24 after first round of models run...
    filter(year > 2007) # 2007 DO data is flawed

  # scale predictors before filtering to ensure mean and SD are global
  data <- data %>% mutate(raw_depth = depth_bath, depth = log(raw_depth), temp = temperature_c, raw_do = do_mlpl, do_mlpl = log(do_mlpl)) 
  do_data <- scale_predictors(data, predictors = c(quo(depth), quo(temp), quo(do_mlpl))) 
       # c(quo(log_depth), quo(mixed), quo(muddy), quo(sandy), quo(rocky), quo(any_rock)))
  data <- do_data %>% mutate(depth = raw_depth) %>% #filter(year != 2016) %>%
    filter(ssid %in% model_ssid) %>% 
    filter(!is.na(depth_scaled)) %>% filter(!is.na(do_mlpl)) %>% filter(!is.na(temp)) 

if (any(names(data) == "adult_density")) { 
  data <- data %>% select(X, Y, X10, Y10,
      adult_density, imm_density,
      depth_scaled, depth_scaled2, depth_mean, depth_sd, 
      do_mlpl_scaled, do_mlpl_scaled2,  do_mlpl_mean, do_mlpl_sd, 
      temp_scaled, temp_scaled2, temp_mean, temp_sd, 
      trawled, any_rock, muddy, sandy, year, ssid, fishing_event_id, 
      depth, do_mlpl, raw_do) 
 } else {
    data <- data %>% select(X, Y, X10, Y10,
      density,
      depth_scaled, depth_scaled2, depth_mean, depth_sd, 
      do_mlpl_scaled, do_mlpl_scaled2,  do_mlpl_mean, do_mlpl_sd, 
      temp_scaled, temp_scaled2, temp_mean, temp_sd, 
      trawled, any_rock, muddy, sandy, year, ssid, fishing_event_id, 
      depth, do_mlpl, raw_do) 
 }

 data <- filter(data, !(year == 2016 & ssid == 4) ) # 2016 do an outlier for survey 4
# unique(data$year)
# table(data$year,data$ssid)

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

  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 = 500)
  }

  sdmTMB::plot_spde(spde)

Run DO sdmTMB model

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

    adult_formula <- as.formula(paste(
      "adult_density ~ 0 + as.factor(year) + do_mlpl_scaled + do_mlpl_scaled2 + temp_scaled + temp_scaled2", covariates, ""
      ))

    starttime1 <- Sys.time()
    adult_biomass_d <- sdmTMB::sdmTMB(data,
      adult_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, 
      time = "year", spde = spde,
      family = tweedie(link = "log"),
      ar1_fields = TRUE,
      include_spatial = TRUE,
      enable_priors = priors,
      control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
      silent = FALSE
    )
    endtime1 <- Sys.time()
    time1 <- round(starttime1 - endtime1)

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


    imm_formula <- as.formula(paste(
      "imm_density ~ 0 + as.factor(year) + do_mlpl_scaled + do_mlpl_scaled2 + temp_scaled + temp_scaled2", covariates, ""
      ))

    starttime2 <- Sys.time()
    imm_biomass_d <- sdmTMB::sdmTMB(data, 
      imm_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, 
      time = "year", spde = spde,
      family = tweedie(link = "log"),
      ar1_fields = FALSE,
      include_spatial = TRUE,
      enable_priors = priors,
      control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
      silent = FALSE
    )
    endtime2 <- Sys.time()
    time2 <- round(starttime2 - endtime2)

    saveRDS(imm_biomass_d, 
      file = paste0("data/", spp, 
        "/mod-imm-biomass-", spp,"-fixed", covs, "-", ssid_string, "-prior-", priors, ".rds"
        ))

  } else {

    starttime3 <- Sys.time()
    dens_formula <- as.formula(paste("density ~ 0 + as.factor(year) + do_mlpl_scaled + do_mlpl_scaled2 + temp_scaled + temp_scaled2", covariates, ""))

    total_biomass_d <- sdmTMB::sdmTMB(data,
      dens_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, 
      time = "year", spde = spde,
      family = tweedie(link = "log"),
      ar1_fields = TRUE,
      include_spatial = TRUE,
      enable_priors = priors,
      control = sdmTMBcontrol(step.min = 0.01, step.max = 1),
      silent = FALSE
    )
    endtime3 <- Sys.time()
    time3 <- round(starttime3 - endtime3)

    saveRDS(total_biomass_d, file = paste0("data/", spp, 
      "/model-total-biomass-", spp, "-fixed", covs, "-", ssid_string, "-prior-", priors, ".rds"
      ))
  }