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.

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.

#read saved sim object from previous example
# sim <- readRDS("optim_para_sim.RDS")
init <- HTM2_sim_optim3
params<-init@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 = init@n[dim(init@n)[1],,], initial_npp = init@n[dim(init@n)[1],])

plotSummary(sim)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: Removed 35 row(s) containing missing values (geom_path).
## Warning: Removed 272 row(s) containing missing values (geom_path).

plotGrowthCurves2(sim, species_panel = T)

plotDiet2(sim)
## Warning: Removed 6151 rows containing missing values (position_stack).

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

# 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 = sim@params@species_params$species,
               gear = sim@params@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 = sim@params@species_params$catchability)

gear_params(sim@params) <- gear_params

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


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

plotSummary(sim2)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: Removed 30 row(s) containing missing values (geom_path).
## Warning: Removed 299 row(s) containing missing values (geom_path).

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).

#read saved data

# f_history2<-as(read.csv("data/fmat.csv", row.names=1), "matrix")[as.character(1967:2019),]
# f_history <- readRDS("data/FmatWeightedInterpolated.rds")
f_history <- FmatWeightedInterpolated
f_history <- f_history[21:73,]
# why is N.pout 2019 mssing??!
f_history$N.pout[53] <- 0.310

# filling the NAs
# assuming that the fisheries effort was relatively constant instead of NAs

for(iSpecies in 2:13)
{
  if(sum(is.na(f_history[,iSpecies])))
  {
 missing_val <- which(is.na(f_history[,iSpecies]))
new_val <- abs(rnorm(length(missing_val),
                 mean(f_history[(missing_val[length(missing_val)]+1):(missing_val[length(missing_val)]+6),iSpecies]),
                 sd(f_history[(missing_val[length(missing_val)]+1):(missing_val[length(missing_val)]+6),iSpecies])))
 f_history[missing_val,iSpecies] <- new_val
  }
}

f_history <- as.matrix(f_history)
rownames(f_history) <- f_history[,1]
f_history <- f_history[,-1]

#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)
##          Sprat    Sandeel    N.pout Herring       Dab   Whiting  Sole  Gurnard
## 1967 0.9700800 0.04183382 0.9905612   0.642 1.3054427 0.6524305 0.403 1.299340
## 1968 0.9734009 0.20421960 1.3014048   0.968 1.1041024 0.6757865 0.510 1.283472
## 1969 0.7089274 0.05930692 1.1535269   0.872 0.9405752 0.6599328 0.552 1.361376
## 1970 0.9278167 0.26006034 0.9434100   0.934 1.1409250 0.6525347 0.512 1.319668
## 1971 0.9266226 0.18326945 0.9567193   1.276 1.3521582 0.6593338 0.488 1.439509
## 1972 0.8522984 0.16528814 0.7193825   0.664 1.1057730 0.6709367 0.522 1.216895
##      Plaice   Haddock  Cod Saithe
## 1967  0.348 0.9099618 0.62  0.352
## 1968  0.350 0.8000865 0.66  0.295
## 1969  0.363 0.8165137 0.62  0.321
## 1970  0.370 0.7698469 0.66  0.345
## 1971  0.372 0.7662074 0.75  0.373
## 1972  0.391 0.8957193 0.81  0.404
# 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(sim2@params)$catchability,"/")
# check
relative_effort[as.character(1988:1992),]
##          Sprat  Sandeel   N.pout   Herring      Dab  Whiting     Sole  Gurnard
## 1988 0.9580546 3.107062 2.198853 0.3312883 5.439474 1.697276 1.337494 6.165474
## 1989 0.2331446 4.203543 2.291268 0.3169734 6.950737 1.801059 1.270620 6.542471
## 1990 1.1332990 3.396453 2.176546 0.2617587 5.575373 1.619440 1.214445 5.882726
## 1991 0.5288214 2.966564 2.068196 0.2852761 5.561428 1.357821 1.206420 4.932379
## 1992 0.6176016 3.242138 1.934353 0.3159509 6.543379 1.279984 1.273295 4.649632
##        Plaice  Haddock      Cod   Saithe
## 1988 3.338142 3.024876 2.569975 2.669992
## 1989 3.165014 3.180763 2.595420 2.591242
## 1990 3.213706 3.014925 2.442748 2.467492
## 1991 3.381425 2.998342 2.417303 2.358743
## 1992 3.365194 2.951907 2.417303 2.291243
# 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),]
##             [,1]        [,2]       [,3]        [,4]       [,5]       [,6]
## 1867 0.000000000 0.000000000 0.00000000 0.000000000 0.00000000 0.00000000
## 1868 0.007564684 0.006490451 0.03188550 0.006630724 0.07265173 0.01424892
## 1869 0.015129369 0.012980903 0.06377100 0.013261449 0.14530347 0.02849784
## 1870 0.022694053 0.019471354 0.09565649 0.019892173 0.21795520 0.04274676
## 1871 0.030258737 0.025961806 0.12754199 0.026522898 0.29060693 0.05699568
## 1872 0.037823421 0.032452257 0.15942749 0.033153622 0.36325866 0.07124460
## 1873 0.045388106 0.038942709 0.19131299 0.039784347 0.43591040 0.08549352
## 1874 0.052952790 0.045433160 0.22319848 0.046415071 0.50856213 0.09974244
## 1875 0.060517474 0.051923612 0.25508398 0.053045795 0.58121386 0.11399136
## 1876 0.068082159 0.058414063 0.28696948 0.059676520 0.65386560 0.12824028
## 1877 0.075646843 0.064904515 0.31885498 0.066307244 0.72651733 0.14248920
## 1878 0.083211527 0.071394966 0.35074047 0.072937969 0.79916906 0.15673812
## 1879 0.090776212 0.077885417 0.38262597 0.079568693 0.87182079 0.17098704
## 1880 0.098340896 0.084375869 0.41451147 0.086199417 0.94447253 0.18523596
## 1881 0.105905580 0.090866320 0.44639697 0.092830142 1.01712426 0.19948488
## 1882 0.113470264 0.097356772 0.47828246 0.099460866 1.08977599 0.21373381
## 1883 0.121034949 0.103847223 0.51016796 0.106091591 1.16242773 0.22798273
## 1884 0.128599633 0.110337675 0.54205346 0.112722315 1.23507946 0.24223165
## 1885 0.136164317 0.116828126 0.57393896 0.119353040 1.30773119 0.25648057
## 1886 0.143729002 0.123318578 0.60582445 0.125983764 1.38038292 0.27072949
## 1887 0.151293686 0.129809029 0.63770995 0.132614488 1.45303466 0.28497841
## 1888 0.158858370 0.136299481 0.66959545 0.139245213 1.52568639 0.29922733
## 1889 0.166423055 0.142789932 0.70148095 0.145875937 1.59833812 0.31347625
## 1890 0.173987739 0.149280383 0.73336645 0.152506662 1.67098986 0.32772517
## 1891 0.181552423 0.155770835 0.76525194 0.159137386 1.74364159 0.34197409
## 1892 0.189117107 0.162261286 0.79713744 0.165768111 1.81629332 0.35622301
## 1893 0.196681792 0.168751738 0.82902294 0.172398835 1.88894505 0.37047193
## 1894 0.204246476 0.175242189 0.86090844 0.179029559 1.96159679 0.38472085
## 1895 0.211811160 0.181732641 0.89279393 0.185660284 2.03424852 0.39896977
## 1896 0.219375845 0.188223092 0.92467943 0.192291008 2.10690025 0.41321869
## 1897 0.226940529 0.194713544 0.95656493 0.198921733 2.17955199 0.42746761
## 1898 0.234505213 0.201203995 0.98845043 0.205552457 2.25220372 0.44171653
## 1899 0.242069897 0.207694447 1.02033592 0.212183182 2.32485545 0.45596545
## 1900 0.249634582 0.214184898 1.05222142 0.218813906 2.39750718 0.47021437
##            [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
## 1867 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## 1868 0.01088910 0.01424892 0.01901795 0.03048602 0.01593544 0.01333329
## 1869 0.02177819 0.02849784 0.03803590 0.06097203 0.03187087 0.02666658
## 1870 0.03266729 0.04274676 0.05705386 0.09145805 0.04780631 0.03999988
## 1871 0.04355638 0.05699568 0.07607181 0.12194406 0.06374174 0.05333317
## 1872 0.05444548 0.07124460 0.09508976 0.15243008 0.07967718 0.06666646
## 1873 0.06533458 0.08549352 0.11410771 0.18291609 0.09561261 0.07999975
## 1874 0.07622367 0.09974244 0.13312567 0.21340211 0.11154805 0.09333304
## 1875 0.08711277 0.11399136 0.15214362 0.24388812 0.12748349 0.10666634
## 1876 0.09800186 0.12824028 0.17116157 0.27437414 0.14341892 0.11999963
## 1877 0.10889096 0.14248920 0.19017952 0.30486015 0.15935436 0.13333292
## 1878 0.11978006 0.15673812 0.20919748 0.33534617 0.17528979 0.14666621
## 1879 0.13066915 0.17098704 0.22821543 0.36583218 0.19122523 0.15999950
## 1880 0.14155825 0.18523596 0.24723338 0.39631820 0.20716067 0.17333280
## 1881 0.15244734 0.19948488 0.26625133 0.42680421 0.22309610 0.18666609
## 1882 0.16333644 0.21373381 0.28526928 0.45729023 0.23903154 0.19999938
## 1883 0.17422554 0.22798273 0.30428724 0.48777624 0.25496697 0.21333267
## 1884 0.18511463 0.24223165 0.32330519 0.51826226 0.27090241 0.22666597
## 1885 0.19600373 0.25648057 0.34232314 0.54874827 0.28683784 0.23999926
## 1886 0.20689282 0.27072949 0.36134109 0.57923429 0.30277328 0.25333255
## 1887 0.21778192 0.28497841 0.38035905 0.60972030 0.31870872 0.26666584
## 1888 0.22867102 0.29922733 0.39937700 0.64020632 0.33464415 0.27999913
## 1889 0.23956011 0.31347625 0.41839495 0.67069233 0.35057959 0.29333243
## 1890 0.25044921 0.32772517 0.43741290 0.70117835 0.36651502 0.30666572
## 1891 0.26133831 0.34197409 0.45643086 0.73166436 0.38245046 0.31999901
## 1892 0.27222740 0.35622301 0.47544881 0.76215038 0.39838589 0.33333230
## 1893 0.28311650 0.37047193 0.49446676 0.79263639 0.41432133 0.34666559
## 1894 0.29400559 0.38472085 0.51348471 0.82312241 0.43025677 0.35999889
## 1895 0.30489469 0.39896977 0.53250267 0.85360842 0.44619220 0.37333218
## 1896 0.31578379 0.41321869 0.55152062 0.88409444 0.46212764 0.38666547
## 1897 0.32667288 0.42746761 0.57053857 0.91458046 0.47806307 0.39999876
## 1898 0.33756198 0.44171653 0.58955652 0.94506647 0.49399851 0.41333205
## 1899 0.34845107 0.45596545 0.60857447 0.97555249 0.50993395 0.42666535
## 1900 0.35934017 0.47021437 0.62759243 1.00603850 0.52586938 0.43999864
relative_effort <- rbind(initial_effort, relative_effort)

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

sim3 <- project(sim2@params, effort = relative_effort, dt = 0.25, t_save = 1)

plotSummary(sim3)
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: Removed 39 row(s) containing missing values (geom_path).
## Warning: Removed 357 row(s) containing missing values (geom_path).

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

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.

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

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

obsy <- read.csv("../../data-raw/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)
## Warning: Removed 1176 row(s) containing missing values (geom_path).
## Warning: Removed 122 rows containing missing values (geom_point).

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.

# 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
## Warning: Removed 106 row(s) containing missing values (geom_path).
## Warning: Removed 9 rows containing missing values (geom_point).

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).

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

params3 <- setParams(sim3@params)

#re-run model
sim4 <- project(params3, 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(sim4)
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
## Warning: Removed 106 row(s) containing missing values (geom_path).
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 106 row(s) containing missing values (geom_path).

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


params3 <- setParams(params3)
#re-run model again...
sim5 <- project(params3, 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(sim5)
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
## Warning: Removed 106 row(s) containing missing values (geom_path).
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 106 row(s) containing missing values (geom_path).

## Warning: Removed 106 row(s) containing missing values (geom_path).

# # 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:

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.