toyexample2.Rmd

---
title: "Example #2: Modelling Changes Through Time"
author: Julia L. Blanchard
date: July 2020
place: Hobart, Australia
output:
  pdf_document: default
  html_document: default
---

# Introduction 

There are many emerging size spectrum modelling (including mizer) applications that aim to examine changes in time series through time. Depending on your question and the goals you have in mind for your model, it may even be worth fitting models to time series data. We may wish to discuss this later. A first step in exploration of ecosystem models with time series however, often starts by simply varying input or "forcing" parameters through time. 

Here, we begin with the steady state or equilibrium model that has already been calibrated and evaluated. 

Presumably these get the model in the correct ball-park for each species time-average biomass, abundance, catches, growth etc. We then examine how different variables can "force"" the model away from the equilibrium state. Often a goal is being asked whether the forcing alone is enough to capture the trends in time series - e.g. fishing mortality, phytoplankton abundance, temperature  include examples that have been published.

Aims of this practical example: 
1) Learn the main steps involved in forcing a size spectrum model
2) Visually compare some of the model predictions with time-series data 
3) Explore how post-hoc parameter changes can affect model skill through time


We previously forced with fishing mortality time series using the North Sea model and there are examples for this in the mizer vignette. This model compared predictions to observations, but we did not capture directional environmental change (only noise in the realised recruitment). One potential issue with the deterministic version of the  model is related to the stock recruitment dynamics we assumed. First, we assumed an eRepro of 1 (which essentially ignores any losses of eggs, and assumes all eggs enter the size spectrum are available to be eaten and potentially grow). The second assumption was related to our values of Rmax. We calibrated the model to catches and biomass and estimated Rmax values (least known parameter). 

PART A. Here we will explore the calibrated model and apply the dynamical forcing.

Preliminary set up again... if needed.
```{r}
#get required packages
library(mizer)
library(tidyverse)

```

Let's read in the saved calibrated parameters of the North sea model stored in the mizer package. These examples do not use the exact same parameters as in the published papers, so are for illustrative purposes here.

```{r,warnings=F}
#read saved sim object from previous example
sim <- readRDS("optim_para_sim.RDS")
params<-sim@params

# run model to equilibrium and plot results, with fishing.
# here an effort = 1 will equate to a fihsing mortality rate = 1 and uses the default knife-edge selectivity function
sim <- project(params, effort = 1, t_max = 200, dt=0.1,initial_n = sim@n[dim(sim@n)[1],,], initial_npp = sim@n[dim(sim@n)[1],])

plot(sim)
plotSpectra(sim,power=2,total =T)
plotlyGrowthCurves(sim,percentage=T)

plotDiet(params,species="Cod")

```

If we agree the model has reached an equilibrium, we can take these equilibrium values (n form last timestep) and set up a dynamical run through time (a simlar example is also shown in the mizer vignette).

#### Forcing the model with fishing mortality rate (F) time series


```{r}

# note we have sigmoidal trawl slectivity parameters, but we need to set up a gear for each species to force each species separately

# to do this, we need to rebuild the params object 
inter <- sim@params@interaction
species_params <- params@species_params
gear_params<-data.frame(species = species_params$species,
               gear = species_params$species,
               sel_func = "sigmoid_length",
               l25 =  c(7.6, 9.8, 8.7, 10.1, 11.5, 19.8, 16.4, 19.8, 11.5,
                        19.1, 13.2, 35.3),
               l50 = c(8.1, 11.8, 12.2, 20.8, 17.0, 29.0, 25.8, 29.0, 17.0,
                       24.3, 22.9, 43.6),
               catchability = params@species_params$catchability)

# reinitiate the params
params <- newMultispeciesParams(species_params,
                                interaction = inter,
                                kappa = params@resource_params$kappa,
                                gear_params = gear_params)


# re-run
sim <- project(params, effort = 1, t_max = 200, dt=0.1,initial_n = sim@n[200,,], initial_npp = sim@n[200,])

plot(sim)
```

### Forcing changes in species' fishing mortality rates through time

Next, we will read in fishing mortality rate time series (note: this matrix is different to the one I used in the paper, and is out of date. It's just for illustration here for consistency with mizer website, needsto be updated with more recent ICES data).

```{r}
#read saved data

f_history<-as(read.csv("data/fmat.csv", row.names=1), "matrix")[as.character(1967:2019),]

#the below one is stored on the mizer website.
#f_history<-as(read.csv("data/NS_f_history.csv", row.names=1), "matrix")
# for some reason these values are different than on the mizer website - need to look into that that, but use for now.

head(f_history)

# express as "relative effort" - here relative to the time-averaged fishing mortality rates used to calirbate the model.

relative_effort <- sweep(f_history,2,gear_params(params)$catchability,"/")
# check
relative_effort[as.character(1988:1992),]


# add a linear ramping up up period from 1867 -1967 
firstRecordedEffort <- NULL
for(iSpecies in 1:ncol(relative_effort))
firstRecordedEffort <- c(firstRecordedEffort,relative_effort[which(!is.na(relative_effort[,iSpecies]))[1],iSpecies]) # first value of F for each species (different years)
# Gurnard has no data so inputing same value as whiting
firstRecordedEffort[8] <- firstRecordedEffort[6]

initial_effort <- matrix(firstRecordedEffort, byrow = TRUE, nrow = 100,
                         ncol = ncol(relative_effort), dimnames = list(1867:1966))

initial_effort ["1867",] <- 0

# doesn't work cause fmat is full of NAs

for (i in 1:12) initial_effort [as.character(1867:1966),i] <- seq(from = initial_effort ["1867",i], to = initial_effort ["1966",i], length = 1966-1867+1)

# check that has worked
initial_effort[as.character(1867:1900),]

relative_effort <- rbind(initial_effort, relative_effort)

# still a bunch of NAs in there, what do we do about it?
##### project model

simt <- project(params, effort = relative_effort, dt = 0.25, t_save = 1)

plot(simt)

# have a look at the modelled yield
plotlyYield(simt)

```

You can zoom in to get a closer look at these in the forcing stage. Here, we are interested in examining the changes along side observations. Let's read in some observe landings for the North Sea and add these to our plot.

```{r}
# output modelled yields and reshape for plotting
y <- getYield(simt)
y <- reshape2::melt(y)

#read in observed yield values (again need to update these data from ICES)

obsy <- read.csv("data/obslandings.csv")[,-1] 
# plot these
ggplot(y) + geom_line(data=y, aes(x = time, y = (value)/1e6, 
            colour = sp)) +
      geom_point(data=obsy,aes(x = time, y = (value)/1e6, 
            colour = sp),size=0.1) +
    facet_wrap(~sp,scales="free_y") +
    scale_y_continuous(name = "ield [g/year]")  +
    scale_colour_manual(values = sim@params@linecolour) +
    xlim(1957, 2011)

```

The trends look kind of OK for some but really not for others.  Remember we re-calibrated this model with completely different assumptions than the before. 

Are the fits in line with our goals for model? They to pass through the cloud of points for some...but not all. Let's have a closer look at a particular species and make sure we use the (less forgiving) linear scale.

```{r}
# look only at  Sprat and examine on linear 

p<-ggplot(y) + geom_line(data=filter(y,sp=="Sprat"), aes(x = time, y = value, 
            colour = sp)) +
      geom_point(data=filter(obsy,sp=="Sprat"),aes(x = time, y = value, 
            colour = sp),size=0.6) +
    #facet_wrap(~sp) +
    scale_y_continuous(name = "Yield [g/year]")  +
    scale_colour_manual(values = sim@params@linecolour) +
    xlim(1965, 2011)
p
```

As expected some of the trends are captured but not the fluctuations. This isn't really suprising, given that the only driver that is changing is fishing (and also the estimates of the fishing mortality rates come from single species stock assessments). Our goal was to cpature trends, hence the fact that the model passes through atleast some of the data points was satisifying our original expectations. 

But we'd really like much better agreement with data here. One issue could be that the erepro values we just re-calibrated the model make the species much more reactive to fishing. Let's examine how sensitive the time series (and their visual agreement to data look when we change our assumptions about eRepro, and possibly Rmax). 

Remember when erepro is 1 essentially all eggs (after density dependent recruitment, Rmax) are available to be eaten and potentially grow.  Explore the consequences of changing erepro at very high (and perhaps very low values of Rmax). 

```{r}

#increase erepro of Sprat
params@species_params$erepro[params@species_params$species=="Sprat"] <- 10^(-0.5)

params <- setReproduction(params)

#re-run model
simt <- project(params, effort=relative_effort, dt = 0.1, t_save = 1,initial_n = sim@n[200,,],initial_n_pp = sim@n_pp[200,])

# output modelled yields and reshape for plotting
y <- getYield(simt)
y <- reshape2::melt(y)

# again look only at this species and examine on linear not log scale
p2 <- p + geom_line(data=filter(y,sp=="Sprat"), aes(x = time, y = value),colour="red") 
p2


#reduce Rmax of Sprat
params@species_params$R_max[params@species_params$species=="Sprat"] <- 1e+12


params <- setParams(params)
#re-run model again...
simt <- project(params, effort=relative_effort, dt = 0.1, t_save = 1,initial_n = sim@n[200,,],initial_n_pp = sim@n_pp[200,])

# output modelled yields and reshape for plotting
y <- getYield(simt)
y <- reshape2::melt(y)

# again look only at Sprat and examine on linear not log scale
p3 <- p2 + geom_line(data=filter(y,sp=="Sprat"), aes(x = time, y = value),colour="blue") 
p3

# # reset Rmax to our calibrated values
# params@species_params$R_max[params@species_params$species=="Sprat"] <- species_params$R_max[params@species_params$species=="Sprat"]
# 
# # reset erepro to our calibrated values
# params@species_params$erepro[params@species_params$species=="Sprat"] <- species_params$erepro[params@species_params$species=="Sprat"]

```


You could use RShiny to do this interactively instead:

```{r}
library(shiny)
runApp("shiny_timeseries")
# need to work out how to input the values
# is there a way to save the final chosen values?
```


Questions: How do these trends differ in terms of the catch time series when these parameters are varied for each species? What happens to the visual "fit"? What does this tell us about density dependence for different species? How would these different values of erepro and Rmax influence the goodness of fit to data if we were to next use a time series fitting approach? Do changes in these values affect the other species dynamics? 

Feel free to explore a little further with your own questions. Then let's come back to the "tips" and our discussion session.



######### JB notes to do: need to update with more recent north sea stock assessment data, here are the catches:

Overall database:
https://standardgraphs.ices.dk/stockList.aspx


Sprat: https://standardgraphs.ices.dk/ViewSourceData.aspx?key=13322

Sandeel : https://standardgraphs.ices.dk/ViewCharts.aspx?key=13303
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13301
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13298
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13304

N.pout - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13166

Herring - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13422

Dab - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13184
 
Whiting - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13525
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13506

Sole - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13743
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13495
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13828

Gurnard - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13493

Plaice - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13484
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13744

Haddock - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13204

Cod - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13740
https://standardgraphs.ices.dk/ViewCharts.aspx?key=13838

Saithe - https://standardgraphs.ices.dk/ViewCharts.aspx?key=13511