SkillAssessInit.Rmd

---
title: "Comparing assessment results with truth: first pass"
author: "Sarah Gaichas and Christine Stawitz"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## So your SS model finally ran!

Now we get to plot outputs using `r4ss`

```{r}
library(here)
library(r4ss)
library(atlantisom)
library(dplyr)
library(ggplot2)
library(ggthemes)
```

```{r species-setup}

#Use the right Atlantis output
source(here("config/CC3Config.R"))
       
#Set species name
species_ss <- "Pacific_sardine" #change to species_ss if we want both 

#Directory with SS files
model_dir <- "Sardine_SS_files/"

#Name of SS data file
datfile_name <- "sardEM_3_3.dat"

#ssinput.data <- r4ss::SS_readdat_3.30(paste0("./inst/extdata/", model_dir, datfile_name))

#Name of SS control file
ctlfile_name <-  "sardEM_3_3.ctl"

# this does not appear to read correctly; has 20 ages for maturity 
# but the model reads it correcty, so don't use for now

#ssinput.ctl <-r4ss::SS_readctl_3.30(paste0("./inst/extdata/", model_dir, ctlfile_name))

```

## Basic SS Plots

```{r ssplots}

replist <- r4ss::SS_output(dir = here(paste0("./inst/extdata/", model_dir)), verbose=TRUE, printstats=TRUE,covar=FALSE)

SS_plots(replist)

replist$timeseries
replist$wtatage
```

## Initial skill assessment

This code is from the multispecies model skill assessment work using Hydra as an operating model and presented in 2015. Original file is hydra_sim_wrapper_modtest.R

```{r skill-functions}

#skill assessment metrics to evaluate model fit

# following from http://heuristically.wordpress.com/2013/07/12/calculate-rmse-and-mae-in-r-and-sas/

# where error = modeled-true or relative error  = modeled/true

# Function that returns Root Mean Squared Error
rmse <- function(error)
{
  sqrt(mean(error^2, na.rm=T))
}

# Function that returns Mean Absolute Error (SAME AS AAE FOR US)
aae <- function(error)
{
  mean(abs(error), na.rm=T)
}

# Function that returns Mean Absolute Error (SAME AS AE FOR US)
ae <- function(error)
{
  mean(error, na.rm=T)
}

# Function that returns Modeling Efficiency MEF
mef <- function(obs, error)
{
  obsavg <- mean(obs, na.rm=T)
  obserr <- obs - obsavg
  (sum(obserr^2)-sum(error^2, na.rm=T))/sum(obserr^2)
}


```

Reading in the truth--these files were generated and saved previously [here](https://sgaichas.github.io/poseidon-dev/SardinesHakeatlantisom2SStest.html). Recruitment and mortality are detailed [here](https://sgaichas.github.io/poseidon-dev/TestCalcZ2.html). 

```{r truth}

# generalized timesteps all models
runpar <- atlantisom::load_runprm(d.name, run.prm.file)
noutsteps <- runpar$tstop/runpar$outputstep
stepperyr <- if(runpar$outputstepunit=="days") 365/runpar$toutinc
midptyr <- round(median(seq(0,stepperyr)))

#Load functional groups
funct.groups <- load_fgs(dir=d.name,
                         file_fgs = functional.groups.file)
#Get just the names of active functional groups
funct.group.names <- funct.groups %>% 
  filter(IsTurnedOn == 1) %>%
  select(Name) %>%
  .$Name

trueB <- readRDS(file.path(d.name, paste0(scenario.name,"surveyBcensus.rds")))

trueB_ss <- trueB %>%
  filter(species==species_ss)

trueN <- readRDS(file.path(d.name, paste0(scenario.name,"surveyNcensus.rds")))

trueN_ss <- trueN %>%
  filter(species==species_ss)

Natage <- readRDS(file.path(d.name, paste0(scenario.name, "natage_census_sard_hake.rds"))) 

Natage_ss <- Natage %>%
  filter(species==species_ss)

lengthwt_census_ss <- readRDS(file.path(d.name, paste0(scenario.name, "lengthwt_census_sard_hake.rds")))

len_ss <- lengthwt_census_ss$natlength %>%
  filter(species==species_ss)

YOY <- atlantisom::load_yoy(d.name, paste0("output", scenario.name, "YOY.txt"))

# load biolprm in some initialize file?
biol <- atlantisom::load_biolprm(d.name, biol.prm.file)

# get code matching species name to split YOY file
code_ss <- funct.groups$Code[which(funct.groups$Name == species_ss)]

# cut to a single species in YOY file
YOY_ss <- YOY %>%
  select(Time, paste0(code_ss, ".0"))

# mg carbon converted to wet weight in tonnes
k_wetdry <- biol$kgw2d / 1000000000

# WARNING only works for CCA because YOY.txt rows 2:end are already in numbers
# and we are merging out the incorrect and irrelevant YOY row 1 (Time=0)
# also hardcoded for sardine example
recnums <- YOY_ss %>%
  mutate(yr = as.integer(round(YOY_ss$Time)/365)) %>%
  mutate(recnums = SAR.0/k_wetdry/biol$redfieldcn) %>%
  filter(yr>0)

file.mort <- file.path(d.name, paste0("output", scenario.name, "Mort.txt"))

mortish <- read.table(file.mort, header = TRUE)

relF_ss <- mortish %>%
  select(Time, relF = paste0(code_ss, ".F"))

```

```{r ests}

replist$timeseries

names(replist$timeseries)[19]<-"Fmort"

```


How does SS3 estimated biomass compare with true?

```{r compareB}

plotB <-ggplot() +
  geom_line(data=trueB_ss, aes(x=time/stepperyr,y=atoutput, color="True B"), 
            alpha = 10/10) +
  geom_point(data=replist$timeseries, aes(x=Yr,y=Bio_all, color="SS3 Est B"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

plotB +
  facet_wrap(~species, scales="free") 


```

What about recruitment?

```{r compareR}

plotR <-ggplot() +
  geom_line(data=recnums, aes(x=yr,y=recnums, color="True R"), 
            alpha = 10/10) +
  geom_point(data=replist$timeseries, aes(x=Yr,y=(Recruit_0*1000), color="SS3 Est R"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

plotR 

```
And F (understanding that F in the mort.txt file from Atlantis is relative, just looking for patterns).

```{r compare-F}

plotF <-ggplot() +
  geom_line(data=relF_ss, aes(x=Time/365, y=relF, color="True F"), 
            alpha = 10/10) +
  geom_point(data=replist$timeseries, aes(x=Yr,y=Fmort, color="SS3 Est F"),
             alpha = 10/10) + 
  theme_tufte() +
  theme(legend.position = "top") +
  labs(colour=scenario.name)

plotF 

```