tastyfish: a package to account for a non-linear functional response when estimating prey dynamics using predator diet data
Matthew Robertson 10/20/2023
tastyfish is a package designed to to estimate prey dynamics by combining bottom-trawl research survey data with predator stomach contents data while accounting for predator functional response. This model’s first application was to estimate sand lance abundance on the Grand Bank (NAFO divisions 3LNO), and the estimated abundance indices can be found as a csv file here as well. This package is associated with Robertson et al. (2022) in which we describe the model’s functionality, limitations, and application potential.
Install the devtools package and run:
devtools::install_github("MatthewRobertson2452/tastyfish")
library(tastyfish)
If you are having problems with installation, you can install the package locally as a ZIP file by clicking the Code menu and “download ZIP” from the github page. You can then extract the folder in a local directory while recording the directory name (which I will reference as download_dir). To install, then use
devtools::install_local(path=download_dir, dependencies=FALSE)
library(tastyfish)
This example will run through how an index was developed for Northern sand lance using stomach contents data from American plaice.
For this example, I will be analyzing trawl, full stomach contents, and called stomach contents data. They are all formatted similarly
library(tastyfish)
data("example_dat")
trawl<-example_dat$trawl
call_sto<-example_dat$call_sto
full_sto<-example_dat$full_sto
head(trawl)
## year pa
## 1 1996 1
## 2 1996 0
## 3 1996 0
## 4 1996 1
## 5 1996 1
## 6 1996 0
We are going to include all of this data into the model. To simplify processes later on, we will concatenate the data into vectors representing the number of samples, presence/absence, and year data. The order in which these are concatenated will matter for some functions later on and should remain consistent throughout.
n=c(length(trawl$year), length(call_sto$year), length(full_sto$year))
pa=c(trawl$pa, call_sto$pa, full_sto$pa)
year=c(trawl$year, call_sto$year, full_sto$year)
Since both called and full stomach contents data will be treated as
having the same functional response we will identify them with the same
ID for the model. This ID will be different from the trawl data. The
order of the numbers in the id
vector need to match the order that the
other data were concatenated.
id=c(0,1,1)
If we wanted the different types of stomach contents data to estimate separate functional response shapes, or if we had one type of stomach contents data for two predators, we would use:
id=c(0,1,2)
We also may want to remember what each data source represents, and will therefore create a vector of names that match the length of our id vector from the prior step. This will not impact the model but it does impact plotting functions.
names=c("Trawl", "Call", "Full")
We then need to organize these data into one dataframe to develop the
data format needed for analyses. This is done by inputting the different
objects into make_data()
.
fishy_dat<-make_data(pa=pa,
year=year,
n=n,
id=id,
names=names)
Now that we have that pre-processing done, we can input that information
with some starting values to create the tmb data and parameter lists
using model_data
and model_param
.
tmb_data<-model_data(fishy_dat=fishy_dat, n=n, type="nonlinear")
param_list<- model_param(tmb_data, lbeta=log(2), lchi=log(0.5), type="nonlinear")
If you are estimating multiple functional response forms (e.g. from multiple predators), lbeta and lchi can be treated as vectors where the order of starting values will match the order that the data was input as.
For example, a model with two predator’s would need to be written as:
tmb_data<-model_data(fishy_dat=fishy_dat, n=n, type="nonlinear")
param_list<- model_param(tmb_data, lbeta=rep(log(2),2), lchi=rep(log(0.5),2), type="nonlinear")
We can then run the TMB model to calculate the prey index of abundance while also estimating and accouting for the functional response of the predator by using the following code:
obj<- TMB::MakeADFun(data = c(model = "NLFPM", # which model to use
tmb_data),
parameters = param_list,
DLL = "tastyfish_TMBExports",
random=c("iye")) # package's DLL
opt<-nlminb(obj$par,obj$fn,obj$gr,control = list(trace=10,eval.max=2000,iter.max=1000),silent=TRUE)
rep<-obj$report()
sdrep<-TMB::sdreport(obj)
I have also included code to quickly visualize the outputs of the model, for the functional response:
plot_curve(rep, fishy_dat)
And time-series:
plot_ts(fishy_dat, rep, sdrep, year=unique(trawl$year),type="nonlinear")
We can compare these results to what we would obtain if assuming a
linear functional response by running the LFPM. We simply remove the
nonlinear functional response shape parameters from the param_list
function and change the type to “linear”:
tmb_data<-model_data(fishy_dat=fishy_dat, n=n, type="linear")
param_list<- model_param(tmb_data, type="linear")
We then change the TMB model to “LFPM”
obj<- TMB::MakeADFun(data = c(model = "LFPM", # which model to use
tmb_data),
parameters = param_list,
DLL = "tastyfish_TMBExports",
random=c("iye")) # package's DLL
opt<-nlminb(obj$par,obj$fn,obj$gr,control = list(trace=10,eval.max=2000,iter.max=1000),silent=TRUE)
rep<-obj$report()
sdrep<-TMB::sdreport(obj)
Finally, we can plot our output be specifying that we used a linear
functional response model in plot_ts
function.
plot_ts(fishy_dat, rep, sdrep, year=unique(trawl$year), type="linear")