code for revision

code/03_fit_nonstationary_models_total_cpue_parallel.R → code/02_fit_nonstationary_models_total_cpue_parallel.R

@@ -1,7 +1,7 @@
 ## Fit Stationary and Non-Stationary Models to zero and positive responses ##
 ## In parallel ##
 
-#remotes::install_github("pbs-assess/sdmTMB@priors-experimental")
+remotes::install_github("pbs-assess/sdmTMB", "nonstationary") # v. 0.0.19.9001
 library(sdmTMB)
 library(dplyr)
 library(future)
@@ -31,8 +31,9 @@ dat$year <- as.factor(dat$year)
 temp <- (dat$depth_m - mean(-grid$depth))
 dat$depth_scaled <- temp / sd(grid$depth)
 dat$depth_scaled2 <- dat$depth_scaled^2
-dat$time <- as.numeric(dat$year) - floor(mean(unique(as.numeric(dat$year))))
+dat$time <- as.numeric(as.character(dat$year)) - floor(mean(unique(as.numeric(as.character(dat$year)))))
 
+#dat$time <- dat$time - mean(dat$time)
 # # spatial + spatiotemporal:
 # fit_models <- function(sub) {
 # spde <- make_mesh(sub, c("lon", "lat"), cutoff = n_cutoff, type = "cutoff")
@@ -61,61 +62,65 @@ dat$time <- as.numeric(dat$year) - floor(mean(unique(as.numeric(dat$year))))
 # AR1 spatiotemporal:
 fit_models_ar1 <- function(sub) {
  spde <- make_mesh(sub, c("lon", "lat"), cutoff = n_cutoff, type = "cutoff")
+ sub$fyear <- as.factor(sub$year)
  ad_fit <- tryCatch({
- sdmTMB(cpue_kg_km2 ~ 0 + depth_scaled + depth_scaled2 + year,
- fields = "AR1",
- include_spatial = FALSE,
+ sdmTMB(cpue_kg_km2 ~ 0 + depth_scaled + depth_scaled2 + fyear,
+ spatiotemporal = "ar1",
+ spatial = "off",
  data = sub,
  time = "year",
- spde = spde,
+ mesh = spde,
  family = tweedie(link = "log"),
  control = sdmTMBcontrol(nlminb_loops = 2, newton_loops = 1),
  priors = sdmTMBpriors(
- phi = halfnormal(0, 10),
- tweedie_p = normal(1.5, 2),
- ar1_rho = normal(0.5, 0.1),
+ #phi = halfnormal(0, 10),
+ #tweedie_p = normal(1.5, 2),
+ #ar1_rho = normal(0.5, 0.1),
  matern_st = pc_matern(range_gt = 10, sigma_lt = 5))
  )}, error = function(e) NA)
  #ad_fit <- refit_model_if_needed(ad_fit)
  ad_fit_ll <- tryCatch({
- sdmTMB(cpue_kg_km2 ~ 0 + depth_scaled + depth_scaled2 + year,
-  fields = "AR1",
-  include_spatial = FALSE,
-  data = sub,
-  time = "year",
-  spde = spde,
-  family = tweedie(link = "log"),
-  epsilon_predictor = "time",
-  control = sdmTMBcontrol(nlminb_loops = 2, newton_loops = 1,
-  lower = list(b_epsilon = -1),
-  upper = list(b_epsilon = 1)),
-  priors = sdmTMBpriors(
-  phi = halfnormal(0, 10),
-  tweedie_p = normal(1.5, 2),
-  ar1_rho = normal(0.5, 0.1),
-  matern_st = pc_matern(range_gt = 10, sigma_lt = 5))
+ sdmTMB(cpue_kg_km2 ~ 0 + depth_scaled + depth_scaled2 + fyear,
+ spatiotemporal = "ar1",
+ spatial = "off",
+ data = sub,
+ time = "year",
+ mesh = spde,
+ family = tweedie(link = "log"),
+ experimental = list(epsilon_predictor = "time", epsilon_model = "trend"),
+ control = sdmTMBcontrol(nlminb_loops = 2, newton_loops = 1,
+ lower = list(b_epsilon = -1),
+ upper = list(b_epsilon = 1)),
+ priors = sdmTMBpriors(
+ #phi = halfnormal(0, 10),
+ #tweedie_p = normal(1.5, 2),
+ #ar1_rho = normal(0.5, 0.1),
+ matern_st = pc_matern(range_gt = 10, sigma_lt = 5))
  )}, error = function(e) NA)
  #ad_fit_ll <- refit_model_if_needed(ad_fit_ll)
  save(ad_fit, ad_fit_ll,
  file = paste0("output/", sub(" ", "_", sub$common_name[[1]]), "_ar1_priors.RData")
  )
 }
 
-refit_model_if_needed <- function(m) {
- if (class(m)!="try-error") {
- if (max(m$gradients) > 0.01) {
- m <- tryCatch({
- sdmTMB::run_extra_optimization(m,
- nlminb_loops = 1L,
- newton_steps = 1L
- )
- }, error = function(e) m)
- }
- }
- m
-}
+# refit_model_if_needed <- function(m) {
+#  if (class(m)!="try-error") {
+#  if (max(m$gradients) > 0.01) {
+#  m <- tryCatch({
+#  sdmTMB::run_extra_optimization(m,
+#  nlminb_loops = 1L,
+#  newton_steps = 1L
+#  )
+#  }, error = function(e) m)
+#  }
+#  }
+#  m
+# }
 
+df = dplyr::group_by(dat, scientific_name) %>%
+ dplyr::summarise(common_name = common_name[1])
 
+species_names = saveRDS(df, file="species.rds")
 split(dat, dat$scientific_name) %>%
  furrr::future_walk(fit_models_ar1)
 # purrr::walk(fit_models_ar1) # serial for testing

code/04_process_model_output.R

@@ -68,11 +68,13 @@ species = dplyr::rename(species,
  common_name = common.name,
  scientific_name = scientific.name)
 
+#species_names = saveRDS(data.frame("species"=unique(dat$scientific_name)), file="species.rds")
+
 for(i in 1:nrow(species)){
 
- comm_name = species$common_name[i]
+ comm_name = species$species[i]
 
- load(file=paste0("output/", sub(" ", "_", comm_name),"_ar1_priors.RData"))
+ load(file=paste0("output/", sub(" ", "_", comm_name), "_ar1_priors.RData"))
 
  df = data.frame(name = comm_name,
  model = c("adult", "adult"),
@@ -81,15 +83,18 @@ for(i in 1:nrow(species)){
  df$trend = NA
  df$trend_se = NA
  df$sigma = NA
-
+ df$common_name <- NA
+ df$scientific_name <- NA
  #if(class(ad_fit)!="try-error") {
  # df$pred_dens[1] = ad_fit$sum_loglik
  #}
- if(class(ad_fit_ll)!="try-error") {
+ if(class(ad_fit_ll)!="logical") {
  #df$pred_dens[2] = ad_fit_ll$sum_loglik
  df$trend[2] = ad_fit_ll$sd_report$value[which(names(ad_fit_ll$sd_report$value) == "b_epsilon")]
  df$trend_se[2] = ad_fit_ll$sd_report$sd[which(names(ad_fit_ll$sd_report$value) == "b_epsilon")]
- df$sigma[2] = ad_fit_ll$sd_report$sd[which(names(ad_fit_ll$sd_report$value) == "sigma_E")[1]]
+ df$sigma[2] = ad_fit_ll$sd_report$value[which(names(ad_fit_ll$sd_report$value) == "log_sigma_E")[1]]
+ df$common_name <- ad_fit_ll$data$common_name[1]
+ df$scientific_name <- ad_fit_ll$data$scientific_name[1]
  }
 
  if(i==1) {

code/05_calc_index.R

@@ -23,8 +23,10 @@ grid = dplyr::mutate(grid,
 grid$cell = seq(1,nrow(grid))
 pred_grid = expand.grid(cell = grid$cell, year = 2003:2018)
 pred_grid = dplyr::left_join(pred_grid, grid)
+pred_grid$fyear = as.factor(pred_grid$year)
+pred_grid$time = as.numeric(as.character(pred_grid$year)) - floor(mean(unique(as.numeric(as.character(pred_grid$year)))))
+#as.numeric(pred_grid$year) - mean(pred_grid$year)
 pred_grid$year = as.factor(pred_grid$year)
-pred_grid$time = as.numeric(pred_grid$year) - floor(mean(unique(as.numeric(pred_grid$year))))
 
 null_predictions_summ = list()
 ll_predictions_summ = list()
@@ -33,24 +35,26 @@ ll_index = list()
 
 for(i in 1:nrow(species)){
 
- load(file=paste0("output/", sub(" ", "_", species$common_name[[i]]),"_ar1_priors.RData"))
+ load(file=paste0("output/", sub(" ", "_", species$species[[i]]),"_ar1_priors.RData"))
 
  est_index = TRUE
- if(class(ad_fit)=="try-error") est_index = FALSE
- if(class(ad_fit_ll)=="try-error") est_index = FALSE
+ if(class(ad_fit)=="logical") est_index = FALSE
+ if(class(ad_fit_ll)=="logical") est_index = FALSE
 
  if(est_index==TRUE) {
- null_predictions <- predict(ad_fit, newdata = pred_grid, sims = 500)
+ null_predictions <- predict(ad_fit, newdata = pred_grid, nsim = 500)
  mean_null <- apply(null_predictions, 1, mean)
  sd_null <- apply(null_predictions, 1, sd)
  null_predictions_summ[[i]] <- cbind(mean_null, sd_null)
  null_index[[i]] <- get_index_sims(null_predictions)
+ null_index[[i]]$common_name <- ad_fit$data$common_name[1]
 
- ll_predictions <- predict(ad_fit_ll, newdata = pred_grid, sims = 500)
+ ll_predictions <- predict(ad_fit_ll, newdata = pred_grid, nsim = 500)
  mean_ll <- apply(ll_predictions, 1, mean)
  sd_ll <- apply(ll_predictions, 1, sd)
  ll_predictions_summ[[i]] <- cbind(mean_ll, sd_ll)
  ll_index[[i]] <- get_index_sims(ll_predictions)
+ ll_index[[i]]$common_name <- ad_fit_ll$data$common_name[1]
 
  # also calculate epsilon_st
  pred_null <- predict(ad_fit, newdata = pred_grid)

example_fig4.Rmd

@@ -43,30 +43,67 @@ ll_predictions = readRDS("output/ll_predictions_summary.rds")
 ```
 
 
-This is just using darkblotched rockfish, and years 2003 - 2018 for comparison
+This is just using lingcod, and years 2003 - 2018 for comparison
+
 ```{r}
-indx = which(species$species=="Darkblotched rockfish")
+indx = which(species$species=="Lingcod")
 # get predictions for each model
- pred_grid$pred_ll = ll_predictions[[indx]][,1]
- pred_grid$pred_null = null_predictions[[indx]][,1]
- # subset years to 2003 and 2018
- pred_grid = dplyr::filter(pred_grid, year %in%c(2003,2018))
-```
+pred_grid$pred_ll_mean = ll_predictions[[indx]][,1]
+pred_grid$pred_null_mean = null_predictions[[indx]][,1]
+pred_grid$pred_ll_se = ll_predictions[[indx]][,2]
+pred_grid$pred_null_se = null_predictions[[indx]][,2]
+# subset years to 2003 and 2018
+pred_grid = dplyr::filter(pred_grid, year %in%c(2003,2018))
 
-First we can make a plot of the absolute difference
+# 2 data frames, one for diff in means, one for diff in ses
+diff_mean = dplyr::group_by(pred_grid, year, cell) %>% 
+ dplyr::summarise(diff = pred_ll_mean - pred_null_mean,
+ lat = lat, lon = lon, metric="mean")
+diff_se = dplyr::group_by(pred_grid, year, cell) %>% 
+ dplyr::summarise(diff = pred_ll_se - pred_null_se,
+ lat = lat, lon = lon, metric="se")
+diff_df = rbind(diff_mean, diff_se)
+```
 
 ```{r}
  # make plots of difference
- g0 = pred_grid %>%
- ggplot(aes(lon,lat,fill=pred_ll - pred_null)) +
+ g1 = diff_mean %>%
+ ggplot(aes(lon,lat,fill=exp(diff))) +
  geom_tile() +
- scale_fill_gradient2() +
- facet_grid(~year) +
+ scale_fill_gradient2(low="red",high="blue",midpoint=1, name="Ratio of means") +
+ facet_wrap(~year) +
  coord_fixed() +
- theme_bw()
+ theme_bw() + 
+ xlab("") + ylab("") + 
+ theme(strip.background =element_rect(fill="white"))
+
+ g2 = diff_se %>%
+ ggplot(aes(lon,lat,fill=exp(diff))) +
+ geom_tile() +
+ scale_fill_gradient2(low="red",high="blue",midpoint=1, name="Ratio of SEs") +
+ facet_wrap(~year) +
+ coord_fixed() +
+ theme_bw() + 
+ xlab("") + ylab("") + 
+ theme(strip.background =element_rect(fill="white"))
+ 
+ library(grid)
+ library(gridExtra)
+ library(cowplot)
+ library(egg)
+ library(ggpubr)
+ #g0 = cowplot::plot_grid(g1,g2,nrow=2)
+ g0=egg::ggarrange(g1,g2,nrow=2)
+ pdf("test_fig_4.pdf", height = 7, width=7)
+ h0 = annotate_figure(g0, 
+ left = textGrob("Latitude", rot = 90, vjust = 13, gp = gpar(cex = 1)),
+ bottom = textGrob("Longitude", hjust=1,vjust = -1,gp = gpar(cex =1)))
+ h0
+ dev.off()
 print(g0)
 ```
 
+
 ```{r}
  # calculate de-meaned data frame, subtracting off spatial mean for each cell
  pred_grid <- dplyr::group_by(pred_grid, cell) %>%