In pbs-assess/gfranges: Groundfish Ranges

Choose one value for each parameter:

Species

species <- params$species

# # Species run so far...
 species <- "Arrowtooth Flounder"
# species <- "Pacific Cod"
# species <- "Sablefish"
# species <- "Silvergray Rockfish"
# species <- "Lingcod"
# species <- "North Pacific Spiny Dogfish" # note: using all data for maturity thresholds
# species <- "Quillback Rockfish"
# species <- "Pacific Ocean Perch"
# species <- "Yelloweye Rockfish"

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 <- ""

# covariates <- "+muddy+mixed+rocky"
# covariates <- "+mixed+rocky"
# covariates <- "+trawled+muddy+rocky+mixed"
# covariates <- "+trawled+mixed+rocky"
# covariates <- "+trawled+mixed"

Run all subsequent code...

covs <- gsub("\\+", "-", covariates)
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)

fish <- readRDS(paste0("raw/bio-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
}
#update_spatial_model <- TRUE

Calculate maturity weighted biomass if not done before

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

if (!exists("biomass")) {
  survey_sets <- readRDS(paste0("raw/event-data-", spp, ""))
  fish <- readRDS(paste0("raw/bio-data-", spp, ""))
  bath <- readRDS("data/bathymetry-data")
  # for POP: exclude sample with unusually large fish 
  # sizes are impossible given recorded catch weight
  # fish <- fish[fish$fishing_event_id!=1506954,]
 rm(maturity)
 try( maturity <- split_catch_maturity(survey_sets, fish, bath,
    # survey = c("SYN QCS","SYN HS","SYN WCVI","SYN WCHG"),
    # survey = c("SYN HS", "SYN WCHG"),
    survey = survey,
    years = years,
    cutoff_quantile = c(.9995),
    plot = TRUE
  ))

  if(!exists("maturity")) {
  maturity <- split_catch_maturity(survey_sets, fish, bath,
    survey = c("SYN QCS","SYN HS","SYN WCVI","SYN WCHG"),
    years = years,
    cutoff_quantile = c(.9995),
    plot = TRUE
  )}

  saveRDS(maturity$data, file = paste0("data/", spp, 
    "/data-by-maturity-", spp, "-", ssid_string, ".rds"
    ))
}

if (exists("maturity")) {
  print(maturity$mass_model)
  ggsave(
    file = paste0("figs/", spp, "/", spp, "-mass-model-", ssid_string, ".png"),
    dpi = 600,
    width = 8, 
    height = 10
  )
}

if (exists("maturity")) {
  png(
    file = paste0("figs/", spp, "/", spp, "-maturity-", ssid_string, ".png"), 
    res = 600,
    units = "in",
    width = 8, 
    height = 10
  )
  print(maturity$maturity)
  dev.off()
}

try(print(maturity$maturity))

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, depth = log(raw_depth), temp = temperature_c) 
  do_data <- gfranges::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(do_mlpl))

opt_data <- data %>% filter(depth > 150) %>% filter(temp < 7.5) 
ggplot(filter(opt_data, adult_density < 0.02), aes(do_mlpl, (adult_density), colour = as.factor(ssid))) + 
  geom_point() +
  facet_wrap(~year, scales="free_y")

#filter(data, adult_density < log(0.002))
ggplot(filter(data, ssid ==16), aes(do_mlpl, (adult_density), colour = as.factor(ssid))) + geom_point() + facet_wrap(~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)", covariates, ""
      ))

    starttime1 <- Sys.time()
    adult_biomass_d <- sdmTMB::sdmTMB(data,
      adult_formula,
      time_varying = ~ 0 + do_mlpl_scaled + do_mlpl_scaled2, # + trawled
      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,"-do3", covs, "-", ssid_string, "-prior-", priors, ".rds"
        ))


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

    starttime2 <- Sys.time()
    imm_biomass_d <- sdmTMB::sdmTMB(data, 
      imm_formula,
      time_varying = ~ 0 + do_mlpl_scaled + do_mlpl_scaled2, # + trawled
      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,"-do3", covs, "-", ssid_string, "-prior-", priors, ".rds"
        ))

  } else {

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

    total_biomass_d <- sdmTMB::sdmTMB(data,
      dens_formula,
      time_varying = ~ 0 + do_mlpl_scaled + do_mlpl_scaled2 + do_mlpl_scaled3, # + trawled
      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, "-do3", covs, "-", ssid_string, "-prior-", priors, ".rds"
      ))
  }

Test gbm

library(gbm)
#glimpse(data)
data$year.f <- as.factor(data$year)
data$log_density_ad <- log(data$adult_density*100)
data$log_do <- log(data$do_mlpl)

 adult_formula <- as.formula(paste(
       "adult_density ~ 0 + year.f * raw_depth + do_mlpl * temperature_c + muddy + sandy + rocky + mixed + X * Y * year.f + raw_depth:do_mlpl + year.f:do_mlpl"
      ))

m <- gbm( adult_formula ,
  # log_density_ad ~ 0 + year.f + raw_depth + muddy + rocky + mixed + # sandy +
  # do_mlpl + temperature_c + X * Y, 
  data = data, 
  distribution = "gaussian",
  n.trees = 2000, interaction.depth = 3, 
  shrinkage = 0.02)

hist(data$log_density_ad)
hist(log(data$do_mlpl))

m

plot(m,i.var=1)
plot(m,i.var=2)
plot(m,i.var=3)
plot(m,i.var=4)
plot(m,i.var=5)
plot(m,i.var=6)
plot(m,i.var=7)
plot(m,i.var=8)
plot(m,i.var=c(1,2))

plot(m,i.var=c(1,3))
plot(m,i.var=c(2,3))
plot(m,i.var=c(3,4))

plot(m,i.var=c(9,10))

data$resid <- predict(m, n.trees = 2000) - data$adult_density

ggplot(data, aes(raw_depth, resid)) +
  geom_point(alpha=0.4) +
  ylim(-0.0025,0.0025) +
  geom_smooth()

ggplot(data, aes(do_mlpl, resid)) +
  geom_point(alpha=0.4) +
  ylim(-0.0025,0.0025) +
  geom_smooth()

ggplot(data, aes(temperature_c, resid)) +
  geom_point(alpha=0.4) +
  ylim(-0.0025,0.0025) +
  geom_smooth()

Run DO sdmTMB model for all regions together

if (region=="All synoptic surveys")
  if (any(names(data) == "adult_density")) {

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

    starttime1 <- Sys.time()
    adult_biomass_d <- sdmTMB::sdmTMB(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,
      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,"-do3", #covs, 
        "-", ssid_string, "-prior-", priors, ".rds"
        ))


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

    starttime2 <- Sys.time()
    imm_biomass_d <- sdmTMB::sdmTMB(data, 
      imm_formula,
      time_varying = ~ 0 + depth_scaled + depth_scaled2, # + trawled
      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,"-do3", #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 + do_mlpl_scaled3"#, covariates, ""
      ))

    total_biomass_d <- sdmTMB::sdmTMB(data,
      dens_formula,
      #time_varying = ~ 0 + do_mlpl_scaled + do_mlpl_scaled2 + do_mlpl_scaled3, # + trawled
      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, "-do3", #covs, 
      "-", ssid_string, "-prior-", priors, ".rds"
      ))
  }

#d <- time_varying_density3(adult_biomass_d, predictor = "do_mlpl")

d <- fixed_density(adult_biomass_d, predictor = "do_mlpl")

do_plot <- plot_mountains(d, variable_label = "Dissolved O2", xlimits = c(0,2)) + 
    ggtitle(" ")

do_plot