example_fig4.Rmd

---
title: "Untitled"
author: "EW"
date: "11/5/2021"
output: pdf_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(sdmTMB)
library(dplyr)
library(ggplot2)
library(viridis)
```

```{r}
# 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)

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$year = as.factor(pred_grid$year)
pred_grid$time = as.numeric(pred_grid$year) - floor(mean(unique(as.numeric(pred_grid$year))))

null_predictions = readRDS("output/null_predictions_summary.rds")
ll_predictions = readRDS("output/ll_predictions_summary.rds")

```


This is just using lingcod, and years 2003 - 2018 for comparison

```{r}
indx = which(species$species=="Lingcod")
# get predictions for each model
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))

# 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
  g1 = diff_mean %>%
    ggplot(aes(lon,lat,fill=exp(diff))) +
    geom_tile() +
    scale_fill_gradient2(low="red",high="blue",midpoint=1, name="Ratio of means") +
    facet_wrap(~year) +
    coord_fixed() +
    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) %>%
    dplyr::mutate(ll_mean = mean(pred_ll),
                  null_mean = mean(pred_null)) %>%
    dplyr::ungroup()
  pred_grid$ll_dev <- pred_grid$pred_ll - pred_grid$ll_mean
  pred_grid$null_dev <- pred_grid$pred_null - pred_grid$null_mean
  pred_grid = dplyr::select(pred_grid, -ll_mean, -null_mean)


  pred_grid2 <- pred_grid
  pred_grid2$est = pred_grid2$null_dev
  pred_grid2$type = "Null"
  pred_grid$est = pred_grid$ll_dev
  pred_grid$type = "Non-stationary"
  pred_grid = rbind(pred_grid, pred_grid2)
  pred_grid$model = paste0(pred_grid$year, pred_grid$type)

# first version: this plots fields that are spatially de-meaned
  g2 =  ggplot(pred_grid, aes(lon,lat,fill=est)) +
    geom_tile() +
    scale_fill_gradient2() +
    facet_grid(year~type) +
    coord_fixed() +
    theme_bw()
print(g2)
```