This is a simple example of how to use \textbf{'MixFishSim'} to generate simulations of the dynamics in a mixed fishery. We describe how to calibrate the habitat fields, the population models, the fishery model and implement a simple fixed spatial closure. \

First, load the packages and set a seed for reproducibility.

Load MixFishSim

library(MixFishSim)
library(knitr)
opts_chunk$set(tidy = TRUE)

set.seed(123)

Initialise the simulation

This vignette is a paired down example of how to construct a simulation using MixFishSim. We include only a basic example and encourage users to explore the other features of the package. \

Base parameters

First we specify the basic parameters of the simulation. This includes the dimensions of the spatial domain, the number of years to simulate, the number of fleets and vessels per fleet and the number of species and how often (in weeks) the fish move.

The object returned is used internally by MixFishSim a list with two levels:

sim <- init_sim(nrows = 10, ncols = 10, n_years = 10, n_tows_day = 4,
        n_days_wk_fished = 5, n_fleets = 2, n_vessels = 20, n_species =
            2, move_freq = 2)

class(sim)
sim$idx
names(sim$brk.idx)

Habitat setup

This function creates the spatial fields which support the fish populations and determine their spatial distributions. You define the parameters for the matern covariance function for each population and optionally the location of any spawning closure areas.

It returns a list of suitable habitat for each species (hab), the habitat as adjusted during the spawning period (spwn_hab) and the binary location of spawning areas (spwn_loc). It also returns the locations as x1,x2,y1,y2 and the multiplier of attractiveness to the spawning area during spawning periods (spwn_mult).

If plot.dist = TRUE, it returns the plots to a file.

hab <- create_hab(sim_init = sim, 
          spp.ctrl = list(
                  "spp.1" = list('nu' = 1/0.015,
                        'var' = 1,
                        'scale' = 1,
                        'Aniso' = matrix(nc = 2, c(1.5, 3, -3, 4))),
                  "spp.2" = list('nu' = 1/0.05,
                        'var'  = 2,
                        'scale' = 12,
                        'Aniso' = matrix(nc = 2, c(1, 2, -1, 2)))
                  ),
          spawn_areas = list(
                     "spp1" = list(
                           'area1' = c(2,5,2,5),
                           'area2' = c(6,8,6,8)
                           ),
                     "spp2" = list(
                           'area1' = c(5,6,6,6)
                           )),
                     spwn_mult = 10, 
                     plot.dist = FALSE)

print(hab)

## Plot the unadjusted habitat fields
plot_habitat(hab$hab)

## Plot the adjusted habitat fields
plot_habitat(hab$spwn_hab)

Population models

Now we need to set up the population models for the simulations. We do this with the init_pop function. We set the initial population biomasses, movement rates, recruitment parameter and growth and natural mortality rates.

The object created stores all the starting conditions and containers for recording the changes in the populations during the simulations.

We can plot the starting distributions for each population as a check.

Pop <- init_pop(sim_init = sim, Bio = c("spp1" = 1e5, "spp2" = 1e5),
        hab = hab[["hab"]], start_cell = c(5,5),
        lambda = c("spp1" = 0.1, "spp2" = 0.1),
        init_move_steps = 20,
        rec_params = list("spp1" = c("model" = "BH", "a" = 54, "b" = 2, "cv" = 0.7),
                  "spp2" = c("model" = "BH", "a" = 27, "b" = 4,"cv" = 0.3)),
        rec_wk = list("spp1" = 3:6, "spp2" = 4:8),
        spwn_wk = list("spp1" = 4:8, "spp2" = 4:8),
      M = c("spp1" = 0.2/365, "spp2" = 0.2/365),
                wt = c("spp1"= 1, "spp2" = 1),
                wtm1 = c("spp1" = 0.1, "spp2" = 0.1),
                K = c("spp1" = 0.3, "spp2" = 0.3)
        )

names(Pop)

Pop$dem_params

image(Pop$Start_pop[[1]], main = "spp1 starting biomass")
image(Pop$Start_pop[[2]], main = "spp2 starting biomass")

Population movement

Now we set up the population tolerance to different temperatures which determines how the populations move during the course of a year. We can then plot the combined spatiotemporal suitable habitat to examine how these interact.

moveCov <- init_moveCov(sim_init = sim, steps = 52,
            spp_tol = list("spp1" = list("mu" = 12, "va" = 8),
                       "spp2" = list("mu" = 15, "va" = 7)
                       )
            )

plot(norm_fun(x = 0:25, mu = 12, va = 8)/max(norm_fun(0:25, 12, 8)), 
     type = "l", xlab = "Temperature", ylab = "Tolerance", lwd = 2)
lines(norm_fun(x = 0:25, mu = 15, va = 7)/ max(norm_fun(0:25, 15, 7)),
      type = "l", col = "blue", lwd = 2)
legend(x = 2, y = 0.9, legend = c("spp1", "spp2"), lwd = 2, col = c("black", "blue"))

plot_spatiotemp_hab(hab = hab, moveCov = moveCov, spwn_wk = list("spp1" = 4:8, "spp2" = 4:8))

Fleet models

Here we initialise the fleet with fish landings price per tonne, catchability coefficients per population, fuel cost, the coefficients for the step function and fleet behaviour.

We can plot the behaviour of the step function to check its suitable for our simulations. This determines the relationship between the monetary value gained from a fishing tow and the next move by the vessel when using the correlated random walk function.

fleets <- init_fleet(sim_init = sim, VPT = list("spp1" = 4, "spp2" = 3),
             Qs = list("fleet 1" = c("spp1" = 1e-5, "spp2" = 3e-5),
                   "fleet 2" = c("spp1" = 5e-5, "spp2" = 1e-5)
                   ),
             fuelC = list("fleet1" = 3, "fleet 2" = 8),
             step_params = list("fleet 1" = c("rate" = 3, "B1" = 1, "B2" = 2, "B3" = 3),
                    "fleet 2" = c("rate" = 3, "B1" = 2, "B2" = 4, "B3" = 4)
                    ),              
             past_knowledge = TRUE,
             past_year_month = TRUE,
             past_trip = TRUE,
             threshold = 0.7
             )

test_step(step_params = fleets$fleet_params[[1]]$step_params, rev.max = 1e2)
test_step(step_params = fleets$fleet_params[[2]]$step_params, rev.max = 1e2)

Spatial closure

We set up a spatial closure. There are multiple options in defining this, but we simply define a static fixed site closure for demonstration purposes.

closure <- init_closure(input_coords = data.frame(x = c(9,10), y = c(6,10)),
            spp1 = "spp1", year_start = 5) 

Survey

Its also possible to define a survey design using the init_survey function, but we do not do so for this demonstration. Please refer to the function help file if this is required.

Run simulation

Finally we run the simulation. The output is a list of objects containing all the information on fisheries catches, the population dynamics and population distributions. These can be examined with some inbuilt plotting functions.

res <- run_sim(sim_init = sim,
           pop_init = Pop,
           move_cov = moveCov,
           fleets_init = fleets,
           hab_init = hab,
           save_pop_bio = TRUE,
           survey = NULL,
           closure = closure)

Summary plots

There are a series of input plotting functions to visualise the results of the simulation. For example, we can explore:

Users will wish to define their own plots, depending on the issues of interest and all the results are saved in the output from the run_sim function.

## Biological
p1 <- plot_pop_summary(results = res,
         timestep = "annual",
         save = FALSE
         )
p1

p2 <- plot_daily_fdyn(res)
p2

## Fishery

logs <- combine_logs(res[["fleets_catches"]])

p3 <- plot_vessel_move(sim_init = sim, logs = logs, fleet_no = 1, vessel_no = 5,
         year_trip = 5, trip_no = 10)
p3

p4 <- plot_realised_stepF(logs = logs, fleet_no = 1, vessel_no = 1)
p4      

Note in our example how the fishing mortality rate for species 2 changes following the spatial closure, which was set to cover some of the core distribution of the population.



pdolder/MixFishSim documentation built on Oct. 17, 2023, 4:25 p.m.