Applied Population Ecology in R

Applied Population Ecology in R

Olaf P. Jensen

Email:

Models

Process

1.  Exponential (r is the intrinsic rate of increase, C is catch)

or , where births (B) =

and deaths (D) =

2.  Logistic (Schaefer), K is the carrying capacity

Observation

1.

Concepts

1.  Observation and state/process models

2.  Simulation and assessment models

3.  Non-linear least squares parameter estimation

4.  Adaptive Management experiments

Questions

1.  Given a reasonably long (25 year) time series of catch and fishing effort data, how accurately can you estimate unfished population size (K), the intrinsic rate of increase (r), the harvest rate that results in maximum sustainable yield (UMSY) and the ratio of current harvest rate to UMSY if there is no process or observation error?

2.  How does your answer change if we add lognormal error to the catch observations?

3.  Can a temporary reduction in fishing effort (i.e., an adaptive management experiment) improve parameter estimates?

Instructions

Numbers correspond to question numbers above.

1. Open the R code that you developed for the logistic model described in the model section above. Recall that there are leading parameters (at the beginning). Add derived parameters at the end. The derived parameters are functions of the leading parameters:

MSY = r*K/4

UMSY = r/2

NMSY = K/2

The simulation model comes first in the code because it generates the data.

Now write the assessment model (after the simulation model). The model structure should be the same as that of the simulation model. When the values of the parameters of the assessment model are the same as those of the simulation model, the predicted biomass and predicted catch should exactly match the true biomass and true catch. Verify this by comparing the numbers and looking at the graph.

Before you attempt to answer the questions, get a feel for how the model behaves. It’s easier if you first remove all fishing mortality by setting the catchability parameter (trueq) to zero. Start the model at low population size by changing the biomass at year one to a low number. Now change the parameters true_r and true_K and see how the model responds. What types of dynamics are possible?

Now let’s create the objective function, i.e., the function that we use to describe how well the model output matches our observations. In this case the “observations” are output from a simulation model. Your objective function in this case should be the sum of squared deviations (SS). The deviations are simply the difference between the “observed” and predicted catches. We can change the assessment model parameters (est_r, est_K, est_q) until we find a combination that minimizes the SS. Estimating parameters in this way is called nonlinear least squares fitting. See if you can find different parameter combinations that result in a good fit between the predicted and observed catches (use the graph). When you find a combination that appears to fit, write down the parameter values. Try to come up with at least four different combinations that all fit pretty well. Is there any pattern to the parameter combinations that fit? What do the different parameter combinations say about MSY and UMSY?

What if instead of trying this brute force approach, we use the optim() function in R to find the best values for the three leading parameters? Type ?optim in the console learn about the syntax for this function. To use optim(), we will first need to create a function that applies the logistic model and returns the sum of squared deviations.

Can optim find a unique solution, i.e., does it find the same parameter combination no matter what values for the three parameters you start with? Does it find the correct values for the three leading parameters?

Now let’s see what happens when we add a small amount of lognormal error to the catch observations using a random number generator function in R. The function is rnorm. Type ?rnorm in the console learn about the syntax for this function. We will implement lognormal error in the catch as follows:

Catch (with error) = Catch (without error) * exp(sigma)

where sigma is a random normal variable with mean zero and standard deviation of 1.

Now can optim find a unique and/or correct solution?

Note: creating random lognormal variables by using EXP(Y), where Y is a random normal variable results in a biased distribution, i.e. if the mean of Y is zero, the mean of EXP(Y) is not exactly one. (If you’re curious, test to see if this is true). Various bias corrections are available.

In the previous exercise, we used a “one-way trip” catch history, where effort increased monotonically while catch peaked and then declined. Such catch histories are typical, and unfortunately not very informative. Let’s see how an adaptive management experiment, where we essentially close the fishery for a few years, effects the stock assessment. Open the worksheet “AM_Experiment.” The layout of this worksheet is identical to “ML_fit_to_CPUE.” Notice that instead of a monotonically increasing effort time series, there is a buildup of effort, during which catch peaks and then begins to decline. Soon after catches decline, effort is reduced to very low levels for three years, after which it resumes at a fixed intermediate level. Compare your ability to estimate the two leading parameters and the derived parameters using this catch history vs. the monotonically increasing catch history when there is observation error. As before, generate one vector of observation errors (using the random number generator) and use it in both models.