diff --git a/code/03_fit_nonstationary_models_total_cpue_parallel.R b/code/02_fit_nonstationary_models_total_cpue_parallel.R
similarity index 63%
rename from code/03_fit_nonstationary_models_total_cpue_parallel.R
rename to code/02_fit_nonstationary_models_total_cpue_parallel.R
index d0fc9cc..301ae2c 100644
--- a/code/03_fit_nonstationary_models_total_cpue_parallel.R
+++ b/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,39 +62,40 @@ 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,
@@ -101,21 +103,24 @@ fit_models_ar1 <- function(sub) {
   )
 }
 
-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
diff --git a/code/04_process_model_output.R b/code/04_process_model_output.R
index 0de7dd2..411d513 100644
--- a/code/04_process_model_output.R
+++ b/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) {
diff --git a/code/05_calc_index.R b/code/05_calc_index.R
index 1aaff97..8d9f5c0 100644
--- a/code/05_calc_index.R
+++ b/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)
diff --git a/code/02_fit_nonstationary_models_total_cpue_parallel_cv.R b/code/archive/02_fit_nonstationary_models_total_cpue_parallel_cv.R
similarity index 100%
rename from code/02_fit_nonstationary_models_total_cpue_parallel_cv.R
rename to code/archive/02_fit_nonstationary_models_total_cpue_parallel_cv.R
diff --git a/example_fig4.Rmd b/example_fig4.Rmd
index 82d9036..841f8ce 100644
--- a/example_fig4.Rmd
+++ b/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) %>%
diff --git a/output/pred_dens_summary.csv b/output/pred_dens_summary.csv
index e5a404c..77b7ece 100644
--- a/output/pred_dens_summary.csv
+++ b/output/pred_dens_summary.csv
@@ -1,77 +1,77 @@
-"","name","model","loglinear","pred_dens","trend","trend_se","sigma"
-"1","sablefish","adult",FALSE,NA,NA,NA,NA
-"2","sablefish","adult",TRUE,NA,-0.00156998867935203,0.00757772066033416,0.122708185936352
-"3","flounder, arrowtooth","adult",FALSE,NA,NA,NA,NA
-"4","flounder, arrowtooth","adult",TRUE,NA,-0.0172068410002163,0.0146720995063858,0.397834211658528
-"5","sanddab, Pacific","adult",FALSE,NA,NA,NA,NA
-"6","sanddab, Pacific","adult",TRUE,NA,-0.0518359210933006,0.0135687797577054,0.431001693597088
-"7","Pacific grenadiers","adult",FALSE,NA,NA,NA,NA
-"8","Pacific grenadiers","adult",TRUE,NA,0.0045798643104747,0.0199109767863361,0.244286977508086
-"9","sole, deepsea","adult",FALSE,NA,NA,NA,NA
-"10","sole, deepsea","adult",TRUE,NA,-0.0834276883622552,0.110351611233545,0.22976819555884
-"11","sole, petrale","adult",FALSE,NA,NA,NA,NA
-"12","sole, petrale","adult",TRUE,NA,0.0275132580709065,0.0130667705594635,0.219445591965961
-"13","sole, rex","adult",FALSE,NA,NA,NA,NA
-"14","sole, rex","adult",TRUE,NA,-0.0549971384627835,0.0184155503565747,0.324098956119735
-"15","sole, flathead","adult",FALSE,NA,NA,NA,NA
-"16","sole, flathead","adult",TRUE,NA,-0.0679363471129145,0.0307502497722832,1.18572730732805
-"17","sole, Dover","adult",FALSE,NA,NA,NA,NA
-"18","sole, Dover","adult",TRUE,NA,-0.0498698237551584,0.02071579408609,0.319046309533818
-"19","lingcod","adult",FALSE,NA,NA,NA,NA
-"20","lingcod","adult",TRUE,NA,0.0406623464966729,0.00818616414251104,0.185522544530894
-"21","sole, English","adult",FALSE,NA,NA,NA,NA
-"22","sole, English","adult",TRUE,NA,0.00127684475152133,0.0139293810621791,0.336632504757078
-"23","turbot, curlfin","adult",FALSE,NA,NA,NA,NA
-"24","turbot, curlfin","adult",TRUE,NA,-0.0369364331645313,0.0370503883120801,1.01559041780341
-"25","skate, longnose","adult",FALSE,NA,NA,NA,NA
-"26","skate, longnose","adult",TRUE,NA,-0.0521500516729773,0.025603683023489,0.24339505880582
-"27","Pacific ocean perch","adult",FALSE,NA,NA,NA,NA
-"28","Pacific ocean perch","adult",TRUE,NA,-0.0578007139096175,0.0183384758382203,0.758487552916459
-"29","rockfish, aurora","adult",FALSE,NA,NA,NA,NA
-"30","rockfish, aurora","adult",TRUE,NA,-0.010504427627849,0.0212629126382578,0.656528880729546
-"31","rockfish, redbanded","adult",FALSE,NA,NA,NA,NA
-"32","rockfish, redbanded","adult",TRUE,NA,-0.0134427445579375,0.0143424481948757,0.378398654491473
-"33","rockfish, greenspotted","adult",FALSE,NA,NA,NA,NA
-"34","rockfish, greenspotted","adult",TRUE,NA,0.00827203439575741,0.0249164964533957,0.63771671536518
-"35","rockfish, darkblotched","adult",FALSE,NA,NA,NA,NA
-"36","rockfish, darkblotched","adult",TRUE,NA,0.0434828215470383,0.0136678195480753,0.415189961873599
-"37","rockfish, splitnose","adult",FALSE,NA,NA,NA,NA
-"38","rockfish, splitnose","adult",TRUE,NA,0.118062842592973,0.0110163025558626,0.308306242778169
-"39","rockfish, greenstriped","adult",FALSE,NA,NA,NA,NA
-"40","rockfish, greenstriped","adult",TRUE,NA,-0.0178260281197252,0.0162809324577883,0.810985539237843
-"41","rockfish, chilipepper","adult",FALSE,NA,NA,NA,NA
-"42","rockfish, chilipepper","adult",TRUE,NA,-0.0175386050346645,0.012537980143781,0.636431647259849
-"43","rockfish, rosethorn","adult",FALSE,NA,NA,NA,NA
-"44","rockfish, rosethorn","adult",TRUE,NA,-0.0574898862378932,0.0363271485581029,0.760826888193704
-"45","rockfish, shortbelly","adult",FALSE,NA,NA,NA,NA
-"46","rockfish, shortbelly","adult",TRUE,NA,-0.0245502082334868,0.0117842915992335,0.667277458152544
-"47","rockfish, melanostomus","adult",FALSE,NA,NA,NA,NA
-"48","rockfish, melanostomus","adult",TRUE,NA,-0.0205232242729431,0.0189428560315607,0.464191090176458
-"49","bocaccio","adult",FALSE,NA,NA,NA,NA
-"50","bocaccio","adult",TRUE,NA,-0.0186605708339269,0.0138072633337549,0.456647096045443
-"51","rockfish, stripetail","adult",FALSE,NA,NA,NA,NA
-"52","rockfish, stripetail","adult",TRUE,NA,0.0028209317775982,0.0143499461965302,0.454180723038511
-"53","rockfish, sharpchin","adult",FALSE,NA,NA,NA,NA
-"54","rockfish, sharpchin","adult",TRUE,NA,-0.0456562288932143,0.0177666112317523,0.865876051787256
-"55","thornyhead, shortspine ","adult",FALSE,NA,NA,NA,NA
-"56","thornyhead, shortspine ","adult",TRUE,NA,0.0403541847270536,0.0304869939046843,0.754381394968426
-"57","thornyhead, longspine","adult",FALSE,NA,NA,NA,NA
-"58","thornyhead, longspine","adult",TRUE,NA,-0.0209190024732427,0.0162396154723142,0.153958851227133
-"59","dogfish, North Pacific","adult",FALSE,NA,NA,NA,NA
-"60","dogfish, North Pacific","adult",TRUE,NA,0.00634737904679518,0.00857870107503903,0.305764627307779
-"61","thornyhead, longspine","adult",FALSE,NA,NA,NA,NA
-"62","thornyhead, longspine","adult",TRUE,NA,-0.0209190024732427,0.0162396154723142,0.153958851227133
-"63","crab, Dungeness","adult",FALSE,NA,NA,NA,NA
-"64","crab, Dungeness","adult",TRUE,NA,-0.0283238887806066,0.0113980053092378,0.340267492998066
-"65","ratfish, spotted","adult",FALSE,NA,NA,NA,NA
-"66","ratfish, spotted","adult",TRUE,NA,-0.0522515492816547,0.00950875923572334,0.318844578377859
-"67","grenadier, giant","adult",FALSE,NA,NA,NA,NA
-"68","grenadier, giant","adult",TRUE,NA,-0.0377487774936798,0.0312263614566731,0.310321801786068
-"69","rockfish, halfbanded","adult",FALSE,NA,NA,NA,NA
-"70","rockfish, halfbanded","adult",TRUE,NA,-0.0436820590316302,0.016103354672798,1.01852972398691
-"71","sole, slender","adult",FALSE,NA,NA,NA,NA
-"72","sole, slender","adult",TRUE,NA,0.0694549062423454,0.0160234208832948,0.424068965537814
-"73","crab, tanner","adult",FALSE,NA,NA,NA,NA
-"74","crab, tanner","adult",TRUE,NA,0.00352188712708748,0.00950742683723024,0.187111579095462
-"75","skate, sandpaper","adult",FALSE,NA,NA,NA,NA
-"76","skate, sandpaper","adult",TRUE,NA,-0.128748323883596,0.0264414923575475,0.419540008093663
+"","name","model","loglinear","pred_dens","trend","trend_se","sigma","common_name","scientific_name"
+"1","Sablefish","adult",FALSE,NA,NA,NA,NA,"Sablefish","Anoplopoma fimbria"
+"2","Sablefish","adult",TRUE,NA,-0.000729247847161701,0.00777757358436007,0.691839647753991,"Sablefish","Anoplopoma fimbria"
+"3","Arrowtooth flounder","adult",FALSE,NA,NA,NA,NA,"Arrowtooth flounder","Atheresthes stomias"
+"4","Arrowtooth flounder","adult",TRUE,NA,-0.0312040845716178,0.0144917998402774,1.73264057378115,"Arrowtooth flounder","Atheresthes stomias"
+"5","Pacific sanddab","adult",FALSE,NA,NA,NA,NA,"Pacific sanddab","Citharichthys sordidus"
+"6","Pacific sanddab","adult",TRUE,NA,-0.0450166829766253,0.0135760354087097,1.48685378255112,"Pacific sanddab","Citharichthys sordidus"
+"7","Pacific grenadiers","adult",FALSE,NA,NA,NA,NA,"Pacific grenadiers","Coryphaenoides acrolepis"
+"8","Pacific grenadiers","adult",TRUE,NA,0.00720577090361807,0.0210489169963226,1.03085454342087,"Pacific grenadiers","Coryphaenoides acrolepis"
+"9","Deepsea sole","adult",FALSE,NA,NA,NA,NA,"Deepsea sole","Embassichthys bathybius"
+"10","Deepsea sole","adult",TRUE,NA,-0.06930589593086,0.0912039666429259,0.836637041358958,"Deepsea sole","Embassichthys bathybius"
+"11","Petrale sole","adult",FALSE,NA,NA,NA,NA,"Petrale sole","Eopsetta jordani"
+"12","Petrale sole","adult",TRUE,NA,0.0410460037455776,0.0137973393078627,0.826615008858185,"Petrale sole","Eopsetta jordani"
+"13","Rex sole","adult",FALSE,NA,NA,NA,NA,"Rex sole","Glyptocephalus zachirus"
+"14","Rex sole","adult",TRUE,NA,-0.0413594578532196,0.01670384325415,1.14157268786767,"Rex sole","Glyptocephalus zachirus"
+"15","Flathead sole","adult",FALSE,NA,NA,NA,NA,"Flathead sole","Hippoglossoides elassodon"
+"16","Flathead sole","adult",TRUE,NA,-0.0627717357740916,0.0313909418301552,2.14904092271672,"Flathead sole","Hippoglossoides elassodon"
+"17","Dover sole","adult",FALSE,NA,NA,NA,NA,"Dover sole","Microstomus pacificus"
+"18","Dover sole","adult",TRUE,NA,-0.0413480266508006,0.0176905840903171,1.22171232373613,"Dover sole","Microstomus pacificus"
+"19","Lingcod","adult",FALSE,NA,NA,NA,NA,"Lingcod","Ophiodon elongatus"
+"20","Lingcod","adult",TRUE,NA,0.0455561853230786,0.00883212379721078,0.79852747444074,"Lingcod","Ophiodon elongatus"
+"21","English sole","adult",FALSE,NA,NA,NA,NA,"English sole","Parophrys vetulus"
+"22","English sole","adult",TRUE,NA,-0.0131122605830874,0.0118681860991206,1.09154767940991,"English sole","Parophrys vetulus"
+"23","Curlfin sole","adult",FALSE,NA,NA,NA,NA,"Curlfin sole","Pleuronichthys decurrens"
+"24","Curlfin sole","adult",TRUE,NA,-0.0181984973835078,0.033497765437781,1.90136287546855,"Curlfin sole","Pleuronichthys decurrens"
+"25","Longnose skate","adult",FALSE,NA,NA,NA,NA,"Longnose skate","Raja rhina"
+"26","Longnose skate","adult",TRUE,NA,-0.049901896751081,0.0240312672568576,0.686235283341633,"Longnose skate","Raja rhina"
+"27","Pacific ocean perch","adult",FALSE,NA,NA,NA,NA,"Pacific ocean perch","Sebastes alutus"
+"28","Pacific ocean perch","adult",TRUE,NA,-0.0624786356759652,0.0209594173923592,2.2164247950147,"Pacific ocean perch","Sebastes alutus"
+"29","Aurora rockfish","adult",FALSE,NA,NA,NA,NA,"Aurora rockfish","Sebastes aurora"
+"30","Aurora rockfish","adult",TRUE,NA,0.0161103705058428,0.0201106287262072,1.30907550908841,"Aurora rockfish","Sebastes aurora"
+"31","Redbanded rockfish","adult",FALSE,NA,NA,NA,NA,"Redbanded rockfish","Sebastes babcocki"
+"32","Redbanded rockfish","adult",TRUE,NA,-0.0297709873875709,0.0170679449642739,1.20696543507631,"Redbanded rockfish","Sebastes babcocki"
+"33","Greenspotted rockfish","adult",FALSE,NA,NA,NA,NA,NA,NA
+"34","Greenspotted rockfish","adult",TRUE,NA,NA,NA,NA,NA,NA
+"35","Darkblotched rockfish","adult",FALSE,NA,NA,NA,NA,"Darkblotched rockfish","Sebastes crameri"
+"36","Darkblotched rockfish","adult",TRUE,NA,0.0580053686562304,0.0144672073204298,1.41518134604214,"Darkblotched rockfish","Sebastes crameri"
+"37","Splitnose rockfish ","adult",FALSE,NA,NA,NA,NA,"Splitnose rockfish ","Sebastes diploproa"
+"38","Splitnose rockfish ","adult",TRUE,NA,0.0947641170527645,0.0105070931176857,0.990079316398321,"Splitnose rockfish ","Sebastes diploproa"
+"39","Greenstriped rockfish","adult",FALSE,NA,NA,NA,NA,NA,NA
+"40","Greenstriped rockfish","adult",TRUE,NA,NA,NA,NA,NA,NA
+"41","Chilipepper rockfish","adult",FALSE,NA,NA,NA,NA,"Chilipepper rockfish","Sebastes goodei"
+"42","Chilipepper rockfish","adult",TRUE,NA,-0.0170489324132962,0.0128632362296588,1.92589740941759,"Chilipepper rockfish","Sebastes goodei"
+"43","Rosethorn rockfish","adult",FALSE,NA,NA,NA,NA,"Rosethorn rockfish","Sebastes helvomaculatus"
+"44","Rosethorn rockfish","adult",TRUE,NA,-0.0742310318614993,0.0433795588191206,1.90575812500007,"Rosethorn rockfish","Sebastes helvomaculatus"
+"45","Shortbelly rockfish","adult",FALSE,NA,NA,NA,NA,"Shortbelly rockfish","Sebastes jordani"
+"46","Shortbelly rockfish","adult",TRUE,NA,-0.0293777302565219,0.0117267621831699,1.95174228338942,"Shortbelly rockfish","Sebastes jordani"
+"47","Melanostomus rockfish","adult",FALSE,NA,NA,NA,NA,"Melanostomus rockfish","Sebastes melanostomus"
+"48","Melanostomus rockfish","adult",TRUE,NA,-0.0161124975878593,0.019814982067417,1.33032342149391,"Melanostomus rockfish","Sebastes melanostomus"
+"49","Bocaccio","adult",FALSE,NA,NA,NA,NA,"Bocaccio","Sebastes paucispinis"
+"50","Bocaccio","adult",TRUE,NA,-0.0173876305154365,0.0151853053072607,1.48803626757189,"Bocaccio","Sebastes paucispinis"
+"51","Stripetail rockfish","adult",FALSE,NA,NA,NA,NA,NA,NA
+"52","Stripetail rockfish","adult",TRUE,NA,NA,NA,NA,NA,NA
+"53","Sharpchin rockfish","adult",FALSE,NA,NA,NA,NA,NA,NA
+"54","Sharpchin rockfish","adult",TRUE,NA,NA,NA,NA,NA,NA
+"55","Shortspine thornyhead ","adult",FALSE,NA,NA,NA,NA,"Shortspine thornyhead ","Sebastolobus alascanus"
+"56","Shortspine thornyhead ","adult",TRUE,NA,0.0379739435923417,0.032191833086369,1.64829340333681,"Shortspine thornyhead ","Sebastolobus alascanus"
+"57","Longspine thornyhead","adult",FALSE,NA,NA,NA,NA,"Longspine thornyhead","Sebastolobus altivelis"
+"58","Longspine thornyhead","adult",TRUE,NA,-0.0114641752109057,0.0166825877198747,0.882400802229726,"Longspine thornyhead","Sebastolobus altivelis"
+"59","Dogfish","adult",FALSE,NA,NA,NA,NA,"Dogfish","Squalus suckleyi"
+"60","Dogfish","adult",TRUE,NA,0.0155414591424771,0.00792438436438408,1.07489107877019,"Dogfish","Squalus suckleyi"
+"61","Longspine thornyhead","adult",FALSE,NA,NA,NA,NA,"Longspine thornyhead","Sebastolobus altivelis"
+"62","Longspine thornyhead","adult",TRUE,NA,-0.0114641752109057,0.0166825877198747,0.882400802229726,"Longspine thornyhead","Sebastolobus altivelis"
+"63","Dungeness crab","adult",FALSE,NA,NA,NA,NA,"Dungeness crab","Cancer magister"
+"64","Dungeness crab","adult",TRUE,NA,-0.0215859685454883,0.0105396879490761,1.33548759136981,"Dungeness crab","Cancer magister"
+"65","Spotted ratfish","adult",FALSE,NA,NA,NA,NA,"Spotted ratfish","Hydrolagus colliei"
+"66","Spotted ratfish","adult",TRUE,NA,-0.0513124290355387,0.00877204700868534,1.35467447367653,"Spotted ratfish","Hydrolagus colliei"
+"67","Giant grenadier","adult",FALSE,NA,NA,NA,NA,"Giant grenadier","Albatrossia pectoralis"
+"68","Giant grenadier","adult",TRUE,NA,-0.043337306216974,0.0281545488371632,0.893296405183462,"Giant grenadier","Albatrossia pectoralis"
+"69","Halfbanded rockfish","adult",FALSE,NA,NA,NA,NA,"Halfbanded rockfish","Sebastes semicinctus"
+"70","Halfbanded rockfish","adult",TRUE,NA,-0.0483025286885406,0.0171877177753199,2.2216528342162,"Halfbanded rockfish","Sebastes semicinctus"
+"71","Slender sole","adult",FALSE,NA,NA,NA,NA,"Slender sole","Lyopsetta exilis"
+"72","Slender sole","adult",TRUE,NA,-0.00830428991895639,0.0221041479850011,3.14314842460068,"Slender sole","Lyopsetta exilis"
+"73","Tanner crab","adult",FALSE,NA,NA,NA,NA,"Tanner crab","Chionoecetes tanneri"
+"74","Tanner crab","adult",TRUE,NA,0.00383094858284821,0.00956661249588154,0.880061112888092,"Tanner crab","Chionoecetes tanneri"
+"75","Sandpaper skate","adult",FALSE,NA,NA,NA,NA,"Sandpaper skate","Bathyraja kincaidii"
+"76","Sandpaper skate","adult",TRUE,NA,-0.136576237198,0.0256476741575999,1.41251586963881,"Sandpaper skate","Bathyraja kincaidii"
diff --git a/plots/Figure_02_trends.pdf b/plots/Figure_02_trends.pdf
index 8eb355c..577748d 100644
Binary files a/plots/Figure_02_trends.pdf and b/plots/Figure_02_trends.pdf differ
diff --git a/plots/Figure_3_biomass_index_top5.pdf b/plots/Figure_3_biomass_index_top5.pdf
index 16019b3..9d13dc2 100644
Binary files a/plots/Figure_3_biomass_index_top5.pdf and b/plots/Figure_3_biomass_index_top5.pdf differ
diff --git a/plots/Figure_4.pdf b/plots/Figure_4.pdf
index b82f999..356b46b 100644
Binary files a/plots/Figure_4.pdf and b/plots/Figure_4.pdf differ
diff --git a/plots/Figure_S3_trends.pdf b/plots/Figure_S3_trends.pdf
index d2fff23..72be9a2 100644
Binary files a/plots/Figure_S3_trends.pdf and b/plots/Figure_S3_trends.pdf differ
diff --git a/plots/Figure_S4_biomass_index_log.pdf b/plots/Figure_S4_biomass_index_log.pdf
index 08cdc0d..11bf4da 100644
Binary files a/plots/Figure_S4_biomass_index_log.pdf and b/plots/Figure_S4_biomass_index_log.pdf differ
diff --git a/plots/Figure_S5_ratio.pdf b/plots/Figure_S5_ratio.pdf
index b10f187..5f72cd9 100644
Binary files a/plots/Figure_S5_ratio.pdf and b/plots/Figure_S5_ratio.pdf differ
diff --git a/plots/Figure_S6_ratio_se.pdf b/plots/Figure_S6_ratio_se.pdf
index f4e916b..69a7e2a 100644
Binary files a/plots/Figure_S6_ratio_se.pdf and b/plots/Figure_S6_ratio_se.pdf differ
diff --git a/plots/Figure_S8.pdf b/plots/Figure_S8.pdf
new file mode 100644
index 0000000..2c9228a
Binary files /dev/null and b/plots/Figure_S8.pdf differ
diff --git a/plots/figure_02_S3.R b/plots/figure_02_S3.R
index ecc3e22..d69b8ee 100644
--- a/plots/figure_02_S3.R
+++ b/plots/figure_02_S3.R
@@ -8,22 +8,32 @@ pdf("plots/Figure_02_trends.pdf", height = 5.5, width = 5.5)
 pred = read.csv("output/pred_dens_summary.csv")
 # remove duplicates in pred
 
-names = read.csv("data/name_change.csv") %>%
-  dplyr::rename(name = old_name)
+pred = pred[-which(pred$common_name=="Longspine thornyhead")[1:2],]
 
-pred = dplyr::left_join(pred,names) %>%
+#names = read.csv("data/name_change.csv") %>%
+#  dplyr::rename(name = old_name)
+
+pred <- pred %>%
   dplyr::mutate(trend_sigma = trend/sigma) %>%
-  dplyr::filter(loglinear==TRUE, name != c("sole, deepsea")) %>%
+  dplyr::filter(loglinear==TRUE, !is.na(trend_se), common_name != "Deepsea sole",
+                common_name != "Giant grenadier") %>%
   dplyr::arrange(trend_sigma)
-pred = pred[which(duplicated(pred$new_name)==FALSE),]
-pred$new_name[which(is.na(pred$new_name))] = "Shortspine thornyhead"
+#pred = dplyr::left_join(pred,names) %>%
+#  dplyr::mutate(trend_sigma = trend/sigma) %>%
+#  dplyr::filter(loglinear==TRUE, name != c("sole, deepsea")) %>%
+#  dplyr::arrange(trend_sigma)
+#pred = pred[which(duplicated(pred$new_name)==FALSE),]
+#pred$new_name[which(is.na(pred$new_name))] = "Shortspine thornyhead"
 #x$name <- factor(x$name, levels = x$name[order(x$val)])
 
+pred$common_name[which(pred$common_name == "Shortspine thornyhead ")] = "Shortspine thornyhead"
+pred$common_name[which(pred$common_name == "Splitnose rockfish ")] = "Splitnose rockfish"
+
 top5_bottom5 = pred[c(1:5, (nrow(pred)-4):nrow(pred)),]
 saveRDS(top5_bottom5,"output/top5_bottom5.rds")
 
-pred$new_name <- factor(pred$new_name, levels = pred$new_name)
 # change order of factor levels
+pred$new_name <- factor(pred$common_name, levels = pred$common_name)
 
 my_col = viridis(1, end=0.8)
 p1 = ggplot(pred, aes(new_name,trend_sigma)) +
@@ -46,23 +56,23 @@ dev.off()
 pdf("plots/Figure_S3_trends.pdf", height = 5.5, width = 5.5)
 
 pred = read.csv("output/pred_dens_summary.csv")
-# remove duplicates in pred
 
-names = read.csv("data/name_change.csv") %>%
-  dplyr::rename(name = old_name)
+pred = pred[-which(pred$common_name=="Longspine thornyhead")[1:2],]
 
-pred = dplyr::left_join(pred,names) %>%
-  dplyr::filter(loglinear==TRUE, name != c("sole, deepsea")) %>%
+pred <- pred %>%
+  dplyr::mutate(trend_sigma = trend/sigma) %>%
+  dplyr::filter(loglinear==TRUE, !is.na(trend_se), common_name != "Deepsea sole",
+                common_name != "Giant grenadier") %>%
   dplyr::arrange(trend)
-pred = pred[which(duplicated(pred$new_name)==FALSE),]
-pred$new_name[which(is.na(pred$new_name))] = "Shortspine thornyhead"
-#x$name <- factor(x$name, levels = x$name[order(x$val)])
+
+pred$common_name[which(pred$common_name == "Shortspine thornyhead ")] = "Shortspine thornyhead"
+pred$common_name[which(pred$common_name == "Splitnose rockfish ")] = "Splitnose rockfish"
 
 #top5_bottom5 = pred[c(1:5, (nrow(pred)-4):nrow(pred)),]
 #saveRDS(top5_bottom5,"output/top5_bottom5.rds")
 
-pred$new_name <- factor(pred$new_name, levels = pred$new_name)
 # change order of factor levels
+pred$new_name <- factor(pred$common_name, levels = pred$common_name)
 
 my_col = viridis(1, end=0.8)
 p1 = ggplot(pred, aes(new_name,trend)) +
diff --git a/plots/figure_03_S4_S5_S6.R b/plots/figure_03_S4_S5_S6.R
index 53e917a..9940ad3 100644
--- a/plots/figure_03_S4_S5_S6.R
+++ b/plots/figure_03_S4_S5_S6.R
@@ -4,46 +4,34 @@ library(viridis)
 
 # for plots, drop 2 elements of the list: Deepsea sole (doesn't converge) and
 # Longspine thornyhead because it was run 2x
-spp_to_drop = c(5, 31)
-# Species of interest
-species = read.csv("survey_data/species_list.csv", fileEncoding="UTF-8-BOM")
-names(species) = tolower(names(species))
-species = dplyr::rename(species,
-                        common_name = common.name,
-                        scientific_name = scientific.name)
-species = species[-c(spp_to_drop),]
-
 
 #---- plot results
 null_index = readRDS("output/null_index_range15_sigma10.rds")
 ll_index = readRDS("output/ll_index_range15_sigma10.rds")
-null_index = null_index[-spp_to_drop]
-ll_index = ll_index[-spp_to_drop]
+#species_names = readRDS(file="species.rds")
+#for(i in 1:length(null_index)) {
+  #if(!is.null(null_index[[i]])) {
+  #  null_index[[i]]$species <- species_names$common_name[i]
+  #}
+  #if(!is.null(ll_index[[i]])) {
+  #  ll_index[[i]]$species <- species_names$common_name[i]
+  #}
+#}
 
 null_df = bind_rows(null_index)
-null_df$species = c(t(replicate(species$common_name,n=length(2003:2018))))
+ll_df = bind_rows(ll_index)
 
 null_df$model = "Constant"
-ll_df = bind_rows(ll_index)
-ll_df$species = c(t(replicate(species$common_name,n=length(2003:2018))))
 ll_df$model = "Log-linear"
 joined_df = rbind(ll_df, null_df)
 
-joined_df = dplyr::filter(joined_df, species != "")
-
-# join in names
-new_names = read.csv("data/name_change.csv") %>%
-  dplyr::rename(species = old_name)
-joined_df = dplyr::left_join(joined_df, new_names)
-joined_df$new_name[which(is.na(joined_df$new_name))] = "Shortspine thornyhead"
-
-joined_df$new_name[which(joined_df$new_name == "Pacific grenadiers")] = "Pacific grenadier"
+joined_df = dplyr::filter(joined_df, common_name != "")
 
 pdf("plots/Figure_S4_biomass_index_log.pdf")
 ggplot(joined_df, aes(year, log_est, fill=model, col=model, group=model)) +
   geom_line(alpha=0.5) +
   geom_ribbon(aes(ymin = log(lwr), ymax = log(upr)), alpha=0.4, colour = NA) +
-  facet_wrap(~ new_name, scale="free_y", ncol = 5) +
+  facet_wrap(~ common_name, scale="free_y", ncol = 5) +
   theme_bw() +
   ylab("Ln biomass index (+/- 2SE)") +
   theme(strip.background =element_rect(fill="white")) +
@@ -66,11 +54,16 @@ joined_df$model = as.factor(joined_df$model)
 
 top5_bottom5 = readRDS("output/top5_bottom5.rds")
 
-sub_df = dplyr::filter(joined_df, new_name %in% top5_bottom5$new_name)
-
 # order based on species
-sub_df$new_name = factor(sub_df$new_name, levels = c("Splitnose rockfish","Lingcod","Slender sole","Petrale sole","Darkblotched rockfish","Dover sole",
-                                                     "Spotted ratfish","Rex sole","Longnose skate","Sandpaper skate"))
+sub_df <- joined_df
+sub_df$common_name[which(sub_df$common_name == "Splitnose rockfish ")] = "Splitnose rockfish"
+sub_df$common_name[which(sub_df$common_name == "Shortspine thornyhead ")] = "Shortspine thornyhead"
+
+sub_df = dplyr::filter(sub_df, common_name %in% top5_bottom5$common_name)
+
+sub_df$new_name = factor(sub_df$common_name, levels = c("Splitnose rockfish","Lingcod","Petrale sole","Darkblotched rockfish",
+                                                     "Shortspine thornyhead"
+                                                     ,"Rex sole","Spotted ratfish","Rosethorn rockfish","Longnose skate","Sandpaper skate"))
 
 pdf("plots/Figure_3_biomass_index_top5.pdf",height = 5,width = 7)
 ggplot(sub_df, aes(year, log_est, fill=model, col=model, group=model)) +
@@ -93,23 +86,23 @@ ggplot(sub_df, aes(year, log_est, fill=model, col=model, group=model)) +
   scale_fill_brewer(palette="Dark2")
 dev.off()
 
-df = null_df[,c("year","est","species","se")]
+df = null_df[,c("year","est","common_name","se")]
 df$est_ll = ll_df$est
 df$se_ll = ll_df$se
 df$ratio = df$est_ll / df$est
 df$ratio_se = df$se_ll / df$se
-df = dplyr::filter(df, species!="")
+df = dplyr::filter(df, common_name!="")
 
 # join in names
-new_names = read.csv("data/name_change.csv") %>%
-  dplyr::rename(species = old_name)
-df = dplyr::left_join(df, new_names)
-df$new_name[which(is.na(df$new_name))] = "Shortspine thornyhead"
+# new_names = read.csv("data/name_change.csv") %>%
+#   dplyr::rename(species = old_name)
+# df = dplyr::left_join(df, new_names)
+# df$new_name[which(is.na(df$new_name))] = "Shortspine thornyhead"
 
 pdf("plots/Figure_S5_ratio.pdf")
-ggplot(df, aes(year, ratio,group=new_name)) +
+ggplot(df, aes(year, ratio,group=common_name)) +
   geom_line() +
-  facet_wrap(~ new_name, scale="free_y", ncol = 5) +
+  facet_wrap(~ common_name, scale="free_y", ncol = 5) +
   theme_bw() +
   ylab("Ratio log-linear estimate / null estimated biomass") +
   theme_bw() +
@@ -122,9 +115,9 @@ dev.off()
 
 
 pdf("plots/Figure_S6_ratio_se.pdf")
-ggplot(df, aes(year, ratio_se,group=new_name)) +
+ggplot(df, aes(year, ratio_se,group=common_name)) +
   geom_line() +
-  facet_wrap(~ new_name, scale="free_y", ncol = 5) +
+  facet_wrap(~ common_name, scale="free_y", ncol = 5) +
   theme_bw() +
   ylab("Ratio of log-linear SE / null biomass SE") +
   theme_bw() +
diff --git a/plots/figure_S8.R b/plots/figure_S8.R
new file mode 100644
index 0000000..710ffc3
--- /dev/null
+++ b/plots/figure_S8.R
@@ -0,0 +1,45 @@
+# Species of interest
+species_names = readRDS(file="species.rds")
+grid = readRDS("data/wc_grid.rds")
+grid = dplyr::rename(grid, lon = X, lat = Y)
+grid = dplyr::mutate(grid,
+                     lon = lon*10, # scale to units of km
+                     lat = lat*10,
+                     depth_scaled = as.numeric(scale(-depth)),
+                     depth_scaled2 = depth_scaled^2) %>%
+  dplyr::select(-log_depth_scaled,
+                -log_depth_scaled2)
+
+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)
+
+null_predictions_summ = list()
+ll_predictions_summ = list()
+null_index = list()
+ll_index = list()
+
+for(ii in 28:nrow(species)){
+print(ii)
+  load(file=paste0("output/", sub(" ", "_", species$common_name[[ii]]), "_ar1_priors.RData"))
+
+  df <- data.frame(resid = residuals(ad_fit_ll, type="mvn-laplace"),
+                   common_name = ad_fit_ll$data$common_name[ii])
+  if(ii==1) {
+    all_df = df
+  } else {
+    all_df = rbind(all_df, df)
+  }
+}
+
+pdf("plots/Figure_S8.pdf")
+ggplot(all_df, aes(sample=resid)) +
+  geom_qq(alpha=0.1) + geom_qq_line(col='blue') +
+  facet_wrap(~common_name) +
+  xlab("Theoretical") + ylab("Sample") +
+  theme(strip.text.x = element_text(size =6))
+dev.off()