In pbs-assess/gfranges: Groundfish Ranges

knitr::opts_chunk$set(dpi = 300,
  fig.width = 5, fig.height = 2.5,
  echo=FALSE, warning=FALSE, message=FALSE)
library(dplyr)
library(ggplot2)
library(gfplot)
library(gfdata)
library(sdmTMB)
library(gfranges)

species <- params$species
priors <- params$priors
covs <- params$covs
log_temp <- params$log_temp
num_params <- params$cov_number
model_w_do <- params$model_w_do
fix_depth <- params$fix_depth

paste("priors =", priors)
paste("model =", covs)
paste("model_w_do =", model_w_do)
paste("log(temperature) =", log_temp)
paste("depth fixed =", fix_depth)

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

Load models with environmental covariates

rm(adult, imm,
  adult_survey1, imm_survey1,
  adult_survey4, imm_survey4,
  adult_survey16, imm_survey16,
  total, total1, total4, total16
  )

try({
  # loads global models with 2016 excluded due to strangely high DO 
  add2016 <- ""
  x2016 <- "x2016"
  # for models using 2016 
  # add2016 <- "-do"
  #  x2016 <- ""
  survey <- c("SYN QCS", "SYN HS", "SYN WCVI", "SYN WCHG")
  model_ssid <- c(1, 3, 4, 16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  try({adult <- readRDS(paste0("data/", spp,
    "/mod-mat-biomass-", spp, "-fixed", add2016, covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  adult<-sdmTMB:::update_model(adult) 
  imm <- readRDS(paste0("data/", spp,
    "/mod-imm-biomass-", spp, "-fixed", add2016, covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  imm <-sdmTMB:::update_model(imm) 
  })
  try({total <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, "-fixed", add2016, covs, "-", ssid_string, "-prior-", priors, ".rds"))
  })
})
if (!exists("adult")) adult <- NULL
if (!exists("imm")) imm <- NULL
if (!exists("total")) total <- NULL

try({
  survey <- c("SYN QCS", "SYN HS")
  model_ssid <- c(1, 3)
  ssid_string <- paste0(model_ssid, collapse = "n")
  try({adult_survey1 <- readRDS(paste0("data/", spp,
    "/mod-mat-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  imm_survey1 <- readRDS(paste0("data/", spp,
    "/mod-imm-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  })
  try({total1 <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"))
  })
  })
if (!exists("adult_survey1")) adult_survey1 <- NULL
if (!exists("imm_survey1")) imm_survey1 <- NULL
if (!exists("total1")) total1 <- NULL

try({
  survey <- c("SYN WCVI")
  model_ssid <- c(4)
  ssid_string <- paste0(model_ssid, collapse = "n")
  try({adult_survey4 <- readRDS(paste0("data/", spp,
    "/mod-mat-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  imm_survey4 <- readRDS(paste0("data/", spp,
    "/mod-imm-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
  ))
  })
  try({total4 <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"))
  })
})
if (!exists("adult_survey4")) adult_survey4 <- NULL
if (!exists("imm_survey4")) imm_survey4 <- NULL
if (!exists("total4")) total4 <- NULL

try({
  survey <- c("SYN WCHG")
  model_ssid <- c(16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  try({adult_survey16 <- readRDS(paste0("data/", spp,
    "/mod-mat-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
   ))
   imm_survey16 <- readRDS(paste0("data/", spp,
    "/mod-imm-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"
   ))
  })
  try({total16 <- readRDS(paste0("data/", spp,
    "/model-total-biomass-", spp, "-fixed-do", covs, "-", ssid_string, "-prior-", priors, ".rds"))
  })
})
if (!exists("adult_survey16")) adult_survey16 <- NULL
if (!exists("imm_survey16")) imm_survey16 <- NULL
if (!exists("total16")) total16 <- NULL


model_list <- list(total1 = total1, total4 = total4, total16 = total16, 
  adult_survey1 = adult_survey1, imm_survey1 = imm_survey1,
  adult_survey4 = adult_survey4, imm_survey4 = imm_survey4,
  adult_survey16 = adult_survey16, imm_survey16 = imm_survey16,
  total = total, adult = adult, imm = imm
)
model_list <- model_list[!sapply(model_list, is.null)]

model_list

Density curves from model coefficients

Calculate DO curves

do <- list()
do_plots <- list()
if (model_w_do) {
for (i in seq_len(length(model_list))) {

 # check if models converged
 if(model_list[[i]]$model$convergence==0){
  do[[i]] <- fixed_density(model_list[[i]], predictor = "do_mlpl", quadratics_only = FALSE)

  if (length(do[[i]])>0) {
   do[[i]]$species <- species
   do[[i]]$age <- ifelse(grepl("imm", names(model_list[i])), "imm", "adult")
   do[[i]]$ssid <- paste0(sort(unique(as.numeric(model_list[[i]]$data$ssid))), collapse = "n")
   do[[i]]$model <- names(model_list[i])
   do[[i]]$covs <- paste0(x2016, covs, priors) 
   roots <- get_quadratic_roots(model_list[[i]], predictor = "do_mlpl", fixed_param = 1, threshold = 0.5)
   do[[i]]$optimal_do <- exp(roots[["optimal"]])
   do[[i]]$lower_do_50 <- exp(roots[["lower_threshold"]])
   do[[i]]$upper_do_50 <- exp(roots[["upper_threshold"]])
   roots75 <- get_quadratic_roots(model_list[[i]], predictor = "do_mlpl", fixed_param = 1, threshold = 0.75)
   do[[i]]$lower_do_75 <- exp(roots75[["lower_threshold"]])
   do[[i]]$upper_do_75 <- exp(roots75[["upper_threshold"]])
   roots25 <- get_quadratic_roots(model_list[[i]], predictor = "do_mlpl", fixed_param = 1, threshold = 0.25)
   do[[i]]$lower_do_25 <- exp(roots25[["lower_threshold"]])
   do[[i]]$upper_do_25 <- exp(roots25[["upper_threshold"]])
   do[[i]]$x <- exp(do[[i]]$x)
   do_plots[[i]] <- 
     plot_mountains(do[[i]], variable_label = "Dissolved O2", xlimits = c(0.2, 8), x_breaks = c(1,2,3,4,5,6,7,8)) +
        ggtitle(paste(species, names(model_list[i]))) +
        geom_vline(xintercept=(do[[i]]$lower_do_50), linetype="longdash", alpha = 0.75) +
        #geom_vline(xintercept=(do[[i]]$upper_do_50), linetype="longdash", alpha = 0.75) +
        geom_vline(xintercept=(do[[i]]$lower_do_75), linetype="dotted", alpha = 0.5) +
        #geom_vline(xintercept=(do[[i]]$upper_do_75), linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=(do[[i]]$lower_do_25), linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=(do[[i]]$optimal_do), size= 1.2, colour = "darkgrey") 
        #geom_vline(xintercept=(do[[i]]$upper_do_25), linetype="dotted", alpha = 0.5) 
   } else {
    do_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
   }
 } else {
  do_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
 }
}

odat <- do.call("rbind", do) 
if (length(odat)>0) {
  odat <- odat %>% group_by(model) %>%
  mutate(do_maxy = max(y_hat)) %>% 
  select(-x, -y_hat, -year) %>% 
  distinct() 
} else {
  odat <- NULL
}
} else {  odat <- NULL }
print(do_plots)

Calculate temperature curves

td <- list()
temp_plots <- list()
for (i in seq_len(length(model_list))) {
 if(model_list[[i]]$model$convergence==0){
  td[[i]] <- fixed_density(model_list[[i]], predictor = "temp", quadratics_only = FALSE, fixed_params = num_params)
  if (length(td[[i]])>0) {
    td[[i]]$species <- species
    td[[i]]$age <- ifelse(grepl("imm", names(model_list[i])), "imm", "adult")
    td[[i]]$ssid <- paste0(sort(unique(as.numeric(model_list[[i]]$data$ssid))), collapse = "n")
    td[[i]]$model <- names(model_list[i])
    td[[i]]$covs <- paste0(x2016, covs, priors)
   roots <- get_quadratic_roots(model_list[[i]], predictor = "temp", fixed_param = num_params, threshold = 0.5)

   if (log_temp) {

   td[[i]]$optimal_temp <- exp(roots[["optimal"]])
   td[[i]]$lower_t_50 <- exp(roots[["lower_threshold"]])
   td[[i]]$upper_t_50 <- exp(roots[["upper_threshold"]])
   roots75 <- get_quadratic_roots(model_list[[i]], predictor = "temp", fixed_param = num_params, threshold = 0.75)
   td[[i]]$lower_t_75 <- exp(roots75[["lower_threshold"]])
   td[[i]]$upper_t_75 <- exp(roots75[["upper_threshold"]])
   roots25 <- get_quadratic_roots(model_list[[i]], predictor = "temp", fixed_param = num_params, threshold = 0.25)
   td[[i]]$lower_t_25 <- exp(roots25[["lower_threshold"]])
   td[[i]]$upper_t_25 <- exp(roots25[["upper_threshold"]])
   td[[i]]$x <- exp(td[[i]]$x)

   } else {

   td[[i]]$optimal_temp <-  (roots[["optimal"]])
   td[[i]]$lower_t_50 <-  (roots[["lower_threshold"]])
   td[[i]]$upper_t_50 <-  (roots[["upper_threshold"]])
   roots75 <- get_quadratic_roots(model_list[[i]], predictor = "temp", fixed_param = num_params, threshold = 0.75)
   td[[i]]$lower_t_75 <-  (roots75[["lower_threshold"]])
   td[[i]]$upper_t_75 <-  (roots75[["upper_threshold"]])
   roots25 <- get_quadratic_roots(model_list[[i]], predictor = "temp", fixed_param = num_params, threshold = 0.25)
   td[[i]]$lower_t_25 <-  (roots25[["lower_threshold"]])
   td[[i]]$upper_t_25 <-  (roots25[["upper_threshold"]])
   td[[i]]$x <-  (td[[i]]$x)

   }

   temp_plots[[i]] <- 
     plot_mountains(td[[i]], variable_label = "Temperature", xlimits = c(2, 12), x_breaks = c(3,4,5,6,7,8,9,10,11)) +
        ggtitle(paste(species, names(model_list[i]))) + 
        geom_vline(xintercept=td[[i]]$lower_t_50, linetype="longdash", alpha = 0.75) +
        geom_vline(xintercept=td[[i]]$upper_t_50, linetype="longdash", alpha = 0.75) +
        geom_vline(xintercept=td[[i]]$lower_t_75, linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=td[[i]]$upper_t_75, linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=td[[i]]$lower_t_25, linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=td[[i]]$upper_t_25, linetype="dotted", alpha = 0.5) +
        geom_vline(xintercept=td[[i]]$optimal_t, size= 1.2, colour = "darkgrey") 
  } else {
    temp_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
  } 
 } else {
    temp_plots[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
  } 
}

tdat <- do.call("rbind", td) 
if (length(tdat)>0) {
  tdat <- tdat %>% group_by(model) %>%
  mutate(t_maxy = max(y_hat)) %>% 
  select(-x, -y_hat, -year) %>% 
  distinct() 
} else {
  tdat <- NULL
}
print(temp_plots)

Check optimals against raw capture denstiies

Combine biomass and sensor data

  model_ssid <- c(1, 3, 4, 16)
  ssid_string <- paste0(model_ssid, collapse = "n")
  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) 

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)
depth_only <- T

if(depth_only) {

try({  adult_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-mat-biomass-", spp, "-trawled-ssid-reml-1n3n4n16-prior-FALSE.rds"
  ))
adult_biomass<-sdmTMB:::update_model(adult_biomass) 
  imm_biomass <- readRDS(paste0(
    "data/", spp,
    "/mod-imm-biomass-", spp, "-trawled-ssid-reml-1n3n4n16-prior-FALSE.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,
    "/mod-total-biomass-", spp, "-trawled-ssid-reml-1n3n4n16-prior-FALSE.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)]

}  else {
  depth_model_list <- model_list
}

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 (in model 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)

Model with environmental covariates

max_depth_found <- max(data[data$present==1,]$raw_depth, na.rm = TRUE)

depth_model_list <- model_list
d <- list()
depth_plots2 <- list()
for (i in seq_len(length(depth_model_list))) {
 if(depth_model_list[[i]]$model$convergence==0){

   if(fix_depth) {
     d[[i]] <- fixed_density(model_list[[i]], predictor = "depth", 
       quadratics_only = FALSE, fixed_params = 3)
   } else {
   d[[i]] <- time_varying_density(depth_model_list[[i]], predictor = "depth")
   }

  if (length(d[[i]]) == 0) {
    depth_plots2[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
  } else {
    d[[i]]$x <- exp(d[[i]]$x)
    depth_plots2[[i]] <- plot_mountains(d[[i]], variable_label = "Depth (controling for DO and temperature)", xlimits = c(0, max_depth_found)) +
      ggtitle(paste(species, names(depth_model_list[i])))
  }
 } else {
    depth_plots2[[i]] <- grid::grid.rect(gp = grid::gpar(col = "white"))
  } 
}
print(depth_plots2)

Plot raw DO

if (model_w_do) {
if(!is.null(data$adult_density)){
odata <- data %>% filter(!is.na(do_mlpl)) %>% 
  mutate(#do_mlpl = raw_do, 
    ad_density_class = ifelse(adult_density==0, 0, round(log10((adult_density*100000)+1))),
    im_density_class = ifelse(imm_density==0, 0, round(log10((imm_density*10000000)+1)))) %>% 
  group_by(ad_density_class) %>% mutate(ad_mean_density = mean(adult_density, na.rm = TRUE)) %>% ungroup() %>%
  group_by(im_density_class) %>% mutate(im_mean_density = mean(imm_density, na.rm = TRUE)) %>% ungroup()

ad_density_class_labs <- unique(signif(odata$ad_mean_density*10000),2) # convert to kg/ha
ad_density_class_labs <- ad_density_class_labs[order(ad_density_class_labs)]
ad_density_class_labs <- formatC(ad_density_class_labs, digits = 1, format = "fg")
ad_density_class_labs

im_density_class_labs <- unique(signif(odata$im_mean_density*10000),2) # convert to kg/ha
im_density_class_labs <- im_density_class_labs[order(im_density_class_labs)]
im_density_class_labs <- formatC(im_density_class_labs, digits = 1, format = "fg")
im_density_class_labs

do_count_ad <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(ad_density_class))) +
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=ad_density_class_labs) + 
  geom_vline(xintercept=odat[odat$model =="adult", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=odat[odat$age =="adult", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Count") +
  ggtitle(paste("Adult", species))

do_prop_ad <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(ad_density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=ad_density_class_labs) + 
  geom_vline(xintercept=odat[odat$model =="adult", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=odat[odat$age =="adult", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Proportion") +
  ggtitle(paste("Adult", species))

do_count_im <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(im_density_class))) +
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=im_density_class_labs) + 
  geom_vline( xintercept = odat[odat$model =="imm", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline( xintercept = odat[odat$age =="imm", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Count") +
  ggtitle(paste("Immature", species))

do_prop_im <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(im_density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=im_density_class_labs) + 
  geom_vline( xintercept = odat[odat$model =="imm", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline( xintercept = odat[odat$age =="imm", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Proportion") +
  ggtitle(paste("Immature", species))

raw_do_plots <- list(do_count_ad, do_prop_ad, do_count_im, do_prop_im )
print(raw_do_plots)

} else {

odata <- data %>% filter(!is.na(do_mlpl)) %>% 
  mutate(density_class = ifelse(density==0, 0, (round(log10((density*10000)+1)*2)+1))) %>% 
  group_by(density_class) %>% mutate(mean_density = mean(density, na.rm = TRUE)) %>% ungroup()

density_class_labs <- unique(signif(odata$mean_density*10000),2) # convert to kg/ha
density_class_labs <- density_class_labs[order(density_class_labs)]
density_class_labs <- formatC(density_class_labs, digits = 1, format = "fg")
density_class_labs


do_count <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(density_class))) +
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=density_class_labs) + 
  geom_vline(xintercept=odat[odat$model =="total", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=odat[odat$age =="adult", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Count") +
  ggtitle(paste("Total", species))

do_prop <- ggplot(odata, aes(x = do_mlpl, stat(count), fill = as.factor(density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=density_class_labs) + 
  geom_vline(xintercept=odat[odat$model =="total", ]$optimal_do, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=odat[odat$age =="adult", ]$optimal_do, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Proportion") +
  ggtitle(paste("Total", species))

raw_do_plots <- list(do_count,do_prop)

print(raw_do_plots)
}
} else {
  raw_do_plots <- NULL
}

Plot raw temperature

if(!is.null(data$adult_density)){
tdata <- data %>% filter(!is.na(temp)) %>% 
  mutate(ad_density_class = ifelse(adult_density==0, 0, round(log10((adult_density*100000)+1))),
    im_density_class = ifelse(imm_density==0, 0, round(log10((imm_density*10000000)+1)))) %>% 
  group_by(ad_density_class) %>% mutate(ad_mean_density = mean(adult_density, na.rm = TRUE)) %>% ungroup() %>%
  group_by(im_density_class) %>% mutate(im_mean_density = mean(imm_density, na.rm = TRUE)) %>% ungroup()

ad_density_class_labs <- unique(signif(tdata$ad_mean_density*10000),2) # convert to kg/ha
ad_density_class_labs <- ad_density_class_labs[order(ad_density_class_labs)]
ad_density_class_labs <- formatC(ad_density_class_labs, digits = 1, format = "fg")
ad_density_class_labs

im_density_class_labs <- unique(signif(tdata$im_mean_density*10000),2) # convert to kg/ha
im_density_class_labs <- im_density_class_labs[order(im_density_class_labs)]
im_density_class_labs <- formatC(im_density_class_labs, digits = 1, format = "fg")
im_density_class_labs

t_count_ad <- ggplot(tdata, aes(x = temp, stat(count),fill = as.factor(ad_density_class))) + 
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=ad_density_class_labs, begin=0, end=1) + 
  geom_vline(xintercept=tdat[tdat$model =="adult", ]$optimal_temp, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=tdat[tdat$age =="adult", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Count") +
  ggtitle(paste("Adult", species))

t_prop_ad <- ggplot(tdata, aes(x = temp, stat(count), fill = as.factor(ad_density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=ad_density_class_labs, begin=0, end=1) + 
  geom_vline(xintercept=tdat[tdat$model =="adult", ]$optimal_temp, color="black", size=1, na.rm = TRUE) +
  geom_vline(xintercept=tdat[tdat$age =="adult", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Proportion") +
  ggtitle(paste("Adult", species))

t_count_im <- ggplot(tdata, aes(x = temp, stat(count), fill = as.factor(im_density_class))) + 
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=im_density_class_labs) + 
  geom_vline( xintercept = tdat[tdat$age =="imm", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Count") +
  ggtitle(paste("Immature", species))

t_prop_im <- ggplot(tdata, aes(x = temp, stat(count), fill = as.factor(im_density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=im_density_class_labs) + 
  geom_vline( xintercept = tdat[tdat$model =="imm", ]$optimal_temp, color="black", size=1, na.rm = TRUE) +
  geom_vline( xintercept = tdat[tdat$age =="imm", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Proportion") +
  ggtitle(paste("Immature", species))

raw_t_plots <- list(t_count_ad, t_prop_ad, t_count_im, t_prop_im )
print(raw_t_plots)

} else {

tdata <- data %>% filter(!is.na(temp)) %>% 
  mutate(density_class = ifelse(density==0, 0, (round(log10((density*10000)+1)*2)+1))) %>% 
  group_by(density_class) %>% mutate(mean_density = mean(density, na.rm = TRUE)) %>% ungroup()

density_class_labs <- unique(signif(tdata$mean_density*10000),2) # convert to kg/ha
density_class_labs <- density_class_labs[order(density_class_labs)]
density_class_labs <- formatC(density_class_labs, digits = 1, format = "fg")
density_class_labs

t_count <- ggplot(tdata, aes(x = temp, stat(count), fill = as.factor(density_class))) +
  geom_density(position = "stack") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=density_class_labs) + 
  geom_vline(xintercept=tdat[tdat$model =="total", ]$optimal_temp, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=tdat[tdat$age =="adult", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Count") +
  ggtitle(paste("Total", species))

t_prop <- ggplot(tdata, aes(x = temp, stat(count), fill = as.factor(density_class))) +
  geom_density(position = "fill") +
  scale_fill_viridis_d(name = "Mean Density \n(kg/ha)", label=density_class_labs) + 
  geom_vline(xintercept=tdat[tdat$model =="total", ]$optimal_temp, color="black",  size=1, na.rm = TRUE) +
  geom_vline(xintercept=tdat[tdat$age =="adult", ]$optimal_temp, color="black", linetype="dashed", size=0.5, na.rm = TRUE) +
  theme_pbs() +
  xlab("Temperature from CTD") +
  ylab("Proportion") +
  ggtitle(paste("Total", species))

raw_t_plots <- list(t_count,t_prop)
print(raw_t_plots)
}

Save DO plots

if(model_w_do) {
png(
  file = paste0("figs/", spp, "/do-", spp, x2016, covs, "-priors-", priors, ".png"),
  res = 600,
  units = "in",
  width = 8.5,
  height = 11
)
gridExtra::grid.arrange(
  grobs = do_plots,
  nrow = 4,
  top = grid::textGrob(paste(species))
)
dev.off()
}

Save temperature plots

png(
  file = paste0("figs/", spp, "/temp-", spp, x2016, covs, "-priors-", priors, ".png"),
  res = 600,
  units = "in",
  width = 8.5,
  height = 11
)
gridExtra::grid.arrange(
  grobs = temp_plots,
  nrow = 4,
  top = grid::textGrob(paste(species))
)
dev.off()

Save plots of raw densities vs model optimals

png(
  file = paste0("figs/", spp, "/check-optimals-", spp, x2016, covs, "-priors-", priors, ".png"),
  res = 600,
  units = "in",
  width = 15,
  height = 6
)
gridExtra::grid.arrange(grobs = c(raw_do_plots, raw_t_plots),
  nrow = 2,
  top = grid::textGrob(paste("Raw densities of", species, "with modeled optimals in vertical lines (if solid includes all surveys)"))
)
dev.off()

Save plots of time-varying depth from climate model

png(
  file = paste0(
    "figs/", spp,
    "/depth-", spp, x2016, covs, "-priors-", priors, "-revised.png"
  ),
  res = 600,
  units = "in",
  width = 8.5,
  height = 11
)
gridExtra::grid.arrange(
  grobs = c(depth_plots, depth_plots2),
  nrow = 4,
  top = grid::textGrob(paste(species))
)
dev.off()

Save optimal values

other_params <- readRDS(paste0("data/root-values-by-species-by-age.rds"))
# # If data far a particular species needs redoing
# other_params <- subset(other_params, species!="Redbanded Rockfish")

if(!is.null(odat)) {
  if(!is.null(tdat)) {
    new_params <- full_join(tdat, odat, by=c("species", "age", "model", "covs", "ssid"))
    all_params <- full_join(other_params, new_params)
  } else {
  try({all_params <- full_join(other_params, odat)})
  }
} else {
  try({all_params <- full_join(other_params, tdat)})
}

all_params <- all_params %>% distinct()
# # or certain rows need removal
# all_params <- all_params[-c(#row number#), ] 
saveRDS(all_params, file = "data/root-values-by-species-by-age.rds")

# spp_params <- all_params[all_params$species==species,]
# View(spp_params)

glimpse(new_params)

Alternative presense-absense code

odata <- data %>% filter(!is.na(do_mlpl)) #%>% filter(do_mlpl < 1.5)

odata1 <- filter(odata, adult_density > 0)
# ymax <- which.max(density(odata1$do_mlpl)$y)
# do_optimal_raw <- density(odata1$do_mlpl)$x[ymax]

oplot <- ggplot(odata1) + 
  geom_density(aes(x = do_mlpl), colour = "blue") +
  # geom_vline(aes(xintercept=do_optimal_raw),
  #           color="blue", linetype="dashed", size=1) +
  geom_density(data=filter(odata, adult_density == 0), aes(x = do_mlpl), colour = "red", inherit.aes = FALSE) + theme_bw() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Frequency") +
  ggtitle(paste("Adult"))
oplot

oploti <- ggplot(filter(odata, imm_density > 0)) + geom_density(aes(x = do_mlpl), colour = "blue") + 
  geom_density(data=filter(odata, imm_density == 0), aes(x = do_mlpl), colour = "red", inherit.aes = FALSE) + theme_bw() +
  xlab("Dissolved O2 from CTD (ml/L)") +
  ylab("Frequency") +
  ggtitle(paste("Immature"))
oploti


tdata <- data %>% filter(!is.na(temp)) #%>% filter(do_mlpl < 1.5)
tplot <-ggplot(filter(tdata, adult_density > 0)) + geom_density(aes(x = temp), colour = "blue") + 
  geom_density(data=filter(tdata, adult_density == 0), aes(x = temp), colour = "red", inherit.aes = FALSE) + theme_bw() +
  xlab("Temperature from CTD") +
  ylab("Frequency") +
  ggtitle(paste("Adult"))
tplot
tploti <-ggplot(filter(tdata, imm_density > 0)) + geom_density(aes(x = temp), colour = "blue") + 
  geom_density(data=filter(tdata, imm_density == 0), aes(x = temp), colour = "red", inherit.aes = FALSE) + theme_bw() +
  xlab("Temperature from CTD") +
  ylab("Frequency") +
  ggtitle(paste("Immature"))
tploti

Save raw presense-absence density plots

png(
  file = paste0("figs/", spp, "/raw-presense-by-climate-", spp, "-density-plots.png"),
  res = 600,
  units = "in",
  width = 8.5,
  height = 11
)
gridExtra::grid.arrange(
  oplot,  tplot, oploti, tploti,
  ncol = 2,
  top = grid::textGrob(paste("Proportion of tows with", species, "(blue) vs without (red)"), gp=grid::gpar(fontsize=16))
)
dev.off()

TRAWLING

Plot time-varying effect of trawling at mean depth

trawl_plots <- list()
for (i in seq_len(length(depth_model_list))) {
  r <- depth_model_list[[i]]$tmb_obj$report()
  b3 <- r$b_rw_t[, 3]
  yrs <- sort(unique(depth_model_list[[i]]$data$year))
  b_j <- depth_model_list[[i]]$model$par[1:length(yrs)]

  SE <- 0.35 # Place-holder

  trawled <- data.frame(year = yrs, intercept = b_j, trawling = b3)

  # trawl_plots[[i]] <- ggplot(trawled, aes(year, y=exp(trawling) ) ) +
  #   geom_point(colour ="red") +
  #   #geom_pointrange(aes(ymin=exp(intercept-0.35), ymax=exp(intercept+0.35))) +
  #   ylab("Difference in kg/ha within trawled zone") + ggtitle(names(depth_model_list[i]))

  trawl_plots[[i]] <- ggplot(trawled, aes(year, exp(intercept))) + geom_point() +
    geom_point(aes(year, y = exp(intercept + trawling)), colour = "red") +
    geom_pointrange(aes(ymin = exp(intercept - (1.96 * SE)), ymax = exp(intercept + (1.96 * SE)))) + # using mean se on year effects
    ylim(min(exp(trawled$intercept - trawled$trawling - (1.96 * SE))), max(exp(trawled$intercept + trawled$trawling + (1.96 * SE)))) +
    ylab("Predicted kg/ha at mean depth") + ggtitle(names(depth_model_list[i]))
}
print(trawl_plots)

Code for calculation SE on time-varying effects...

# FIXME: need to add to sdmTMB!
r <- m$tmb_obj$report()
summary(m)
yrs <- sort(unique(m$data$year))
.sd <- summary(m$sd_report)
.sd <- as.data.frame(.sd)
.sd$par <- row.names(.sd)
.sd <- .sd[grepl("^b_rw_t", .sd$par), ]
varnames <- colnames(model.matrix(m$time_varying, m$data))

rep(varnames, each = length(yrs))
.sd$parameter <- rep(varnames, each = length(yrs))
.sd$year <- rep(seq_along(yrs), each = length(varnames))
row.names(.sd) <- NULL
.sd$par <- NULL

names(.sd$Estimate) <- "estimate"
names(.sd[2]) <- "se"
.sd$ci <- .sd$se * 1.96
m$time_varying


nrow(.sd)

b3 <- r$b_rw_t[, 3]
yrs <- sort(unique(m$data$year))
b_j <- m$tmb_params$b_j[1:8]

trawled <- data.frame(year = yrs, intercept = b_j, trawling = b3)

pbs-assess/gfranges documentation built on Dec. 13, 2021, 4:50 p.m.