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" covs <- gsub("\\+", "-", covariates) covs <- "-tv-depth"
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) # 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) do_data <- scale_predictors(data, predictors = c(quo(depth), quo(temp), quo(do_mlpl))) data <- do_data %>% mutate(depth = raw_depth) %>% # filter(year != 2016) %>% filter(ssid %in% model_ssid) %>% filter(!is.na(raw_depth)) %>% filter(!is.na(temp)) #%>% select(X, Y, X10, Y10, adult_density, imm_density, depth_scaled, depth_scaled2, temp_scaled, temp_scaled2, depth_mean, depth_sd, temp_mean, temp_sd, trawled, any_rock, muddy, sandy, year, ssid, fishing_event_id, depth, temp)
if (region == "All synoptic surveys") { spde <- sdmTMB::make_spde(data$X, data$Y, n_knots = 800) } 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) } sdmTMB::plot_spde(spde)
Run temperature sdmTMB model
if (any(names(data) == "adult_density")) { adult_formula <- as.formula(paste( "adult_density ~ 0 + as.factor(year) + temp_scaled + temp_scaled2", covariates, "" )) starttime1 <- Sys.time() adult_biomass_t <- 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_t, file = paste0("data/", spp, "/mod-mat-biomass-", spp,"-fixed-temp", covs, "-", ssid_string, "-prior-", priors, ".rds" )) imm_formula <- as.formula(paste( "imm_density ~ 0 + as.factor(year) + temp_scaled + temp_scaled2", #+ depth_scaled + depth_scaled2", covariates, "" )) starttime2 <- Sys.time() imm_biomass_t <- 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_t, file = paste0("data/", spp, "/mod-imm-biomass-", spp,"-fixed-temp", covs, "-", ssid_string, "-prior-", priors, ".rds" )) } else { starttime3 <- Sys.time() dens_formula <- as.formula(paste("density ~ 0 + as.factor(year) + temp_scaled + temp_scaled2", #+ depth_scaled + depth_scaled2", covariates, "")) total_biomass_t <- 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_t, file = paste0("data/", spp, "/model-total-biomass-", spp, "-fixed-temp", covs, "-", ssid_string, "-prior-", priors, ".rds" )) }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.