Tutorial: Time series fitting
Learning Objectives
- Obtain experience with the procedures for fitting an ecosystem model to time series data
- Develop and understanding of the difference between drivers and reference time series
- Develop and understanding about the impact of missing drivers
The purpose of this exercise is to explore how to fit an ecosystem model to time series data. For this, we use the ecosystem model of Anchovy Bay that we constructed in a previous tutorial (download). When gathering time series for use in Ecosim, we construct a spreadsheet, and save it as a comma-separated value (CSV) file that can be read by the time-dynamic module of EwE, Ecosim.
The CSV file must have a specific format, see Table 1. You can copy this table to Excel and save it as a CSV file.
Table 1. Time series file for the Anchovy Bay tutorial. You can copy the content below, paste it (Paste Special > Text) into MS Excel, then save it as a csv file. Import the csv file to EwE as described in the tutorial.
Title | Sealers | Seal B | Trawlers | Cod B | Whiting B | Shrimp C | T bottom | dummy |
Weight | 0 | 1 | 0 | 1 | 1 | 1 | 1 | 1 |
Pool code | 1 | 2 | 2 | 3 | 4 | 7 | 5 | 4 |
Pool code 2 | ||||||||
Type | 3 | 0 | 3 | 0 | 0 | 6 | 2 | 2 |
1970 | 1 | 1 | 1 | 10 | 1 | 0.3 | 17.1 | 1 |
1971 | 1 | 1.05 | 17 | 1 | ||||
1972 | 1 | 1.103 | 16.9 | 1 | ||||
1973 | 0.75 | 1.158 | 16.7 | 1 | ||||
1974 | 0.5 | 1.216 | 17 | 1 | ||||
1975 | 0.25 | 0.8 | 1.276 | 16.6 | 1 | |||
1976 | 0 | 1.34 | 16.4 | 1 | ||||
1977 | 0 | 1.407 | 16.2 | 1 | ||||
1978 | 0 | 1.477 | 6 | 0.8 | 16.3 | 1 | ||
1979 | 0 | 1.551 | 16.1 | 1 | ||||
1980 | 0 | 1 | 1.629 | 16.3 | 1 | |||
1981 | 0 | 1.71 | 16.4 | 1 | ||||
1982 | 0 | 1.796 | 16.1 | 1 | ||||
1983 | 0 | 1.886 | 16 | 1 | ||||
1984 | 0 | 1.98 | 16.2 | 1 | ||||
1985 | 0 | 2.079 | 16.5 | 1 | ||||
1986 | 0 | 2.183 | 16.6 | 1 | ||||
1987 | 0 | 2.292 | 16.6 | 1 | ||||
1988 | 0 | 2.407 | 16.8 | 1 | ||||
1989 | 0 | 2.527 | 4 | 0.7 | 17 | 1 | ||
1990 | 0 | 2.577 | 17.2 | 1 | ||||
1991 | 0 | 2.629 | 17.6 | 1 | ||||
1992 | 0 | 2 | 2.682 | 17.8 | 1 | |||
1993 | 0 | 2.735 | 18 | 1 | ||||
1994 | 0 | 2.79 | 18.2 | 1 | ||||
1995 | 0 | 2.846 | 18.6 | 1 | ||||
1996 | 0 | 2.903 | 18.7 | 1 | ||||
1997 | 0 | 2.961 | 2 | 0.6 | 18.6 | 1 | ||
1998 | 0 | 3.02 | 18 | 1 | ||||
1999 | 0 | 3.08 | 18.9 | 1 | ||||
2000 | 0 | 3.142 | 19 | 1 | ||||
2001 | 0 | 3.205 | 19.2 | 1 | ||||
2002 | 0 | 3 | 3.269 | 19 | 1 | |||
2003 | 0 | 3.334 | 19.3 | 1 | ||||
2004 | 0 | 3.401 | 19.2 | 1 | ||||
2005 | 0 | 3.469 | 1 | 0.4 | 2 | 19 | 1 | |
2006 | 0 | 3.5 | 18.6 | 1 | ||||
2007 | 0 | 3.5 | 18.5 | 1 | ||||
2008 | 0 | 3.55 | 18.7 | 1 | ||||
2009 | 0 | 3.6 | 18.6 | 1 | ||||
2010 | 0 | 4 | 3.65 | 2.4 | 18 | 1 |
The first row gives titles of the time series, the second row is optional and gives weights that will be used as weighting factors for the estimation of summed squared residuals (SS). The pool code is the fleet number (for effort, here Sealers and Trawlers), the functional group number (for biomasses and mortalities, here Seal, Cod, Whiting, and Shrimp), or a forcing function number, (here “dummy”). The type is a code, explained in the time series chapter of the User Guide – it can be replaced by descriptive abbreviations as explained in the User Guide, e.g., BiomassRel instead of type 0. The following rows give the time series data by year. Notice that drivers (effort and fishing mortalities) should be given for all years, (or effort will be set to 0 for missing years), while reference time series can be for some years only.
For this tutorial, you will need a time series file, so copy the data in Table 1 to a blank Excel spreadsheet and save it as a CSV file.
Open EwE, then open (Menu > File > Open model) the Anchovy Bay database and model in Ecopath. Proceed to Ecosim > Input > Time series where you can import your time series file by selecting Import, and browse to the anchovybay.csv file. (If needed, you can download the CSV file from this link).
When the time series tab opens, check out each of the time series (see thumbnails at the bottom). We have effort, relative biomasses, and a catch series.
Before you start fitting the model, make a run of Ecosim (Ecosim > Output > Run Ecosim > Run. Examine the plots and notably work your way through Ecosim > Output > Ecosim group plots. To illustrate what to look for, examine the seals. You’ll see that the Ecosim run (line) shows a small increase in seal biomass, while the time series (dots) indicates that seals have quadrupled. This indicates that the vulnerability setting for seals is too low, i.e. that seals with the default vulnerability of 2 are assumed to be too close to their carrying capacity. [Read the density dependence chapter for more].
The default vulnerability multiplier setting (2.0) tells Ecosim that seals can at most double the predation mortality they are causing on their prey, while the time series says they have increased their biomass four times. Can you see the discrepancy?
You can try to increase vulnerability for seals as predators and see what happens, and you should go back to run Ecosim ever so often, while doing the fitting to find out what is happening.
In Ecosim now do the following,
- Ecosim > Input > Vulnerabilities form, click the top left cell to select all cells. Then enter 2 in the top right field to the left of where it says Apply, then click Apply, to reset all vulnerabilities to the default 2.
- To find the vulnerability multipliers that have the biggest impact on model fit (summed squared residuals, SS, see the estimating vulnerabilities chapter), go to Ecosim > Tools > Fit to time series, set the No of blocks to 1, then click Sensitivity of SS to V (V is the vulnerability multiplier, see the vulnerability and vulnerability multipliers chapter). Leave by predator checked, then click Search. The search routine will find the group for which the vulnerability parameter has the biggest impact on the SS. Click OK, and this information will be passed on to the search form.
- On the search form, click Search, wait till the search routine has converged, click no to further searching. Note the SS and AICc estimates.
Repeat the steps above with increasing number of search blocks and parameters until you are searching for all predators, then try searching predator-prey combinations, and carry on until you have established a pattern.
Up to now, we have evaluated density-dependent effects, including how ecological factors and fisheries impacts fitting, next is to add environmental effects by estimating ‘primary production anomalies‘ that let system productivity vary over time. We can do that by searching for spline points. For this we need a forcing function, and we can use the “dummy”, go Ecosim > Input > Forcing function > Apply FF (primary producer) , and click the phytoplankton row. Then select the “dummy” and use the green arrow to move it to applied functions.
Back on the Ecosim > Output > Tools > Fit to time series form, uncheck the vulnerability search and check the anomaly search. Now,
- Click the Forcing function tab. Then click Reset all to reset the forcing function.
- Click the Anomaly Search tab on the top. Then set the number of spline points to 2.
- Then click Search, and wait till the search routine has converged, click no to further searching when prompted. Note the SS and AICc estimates.
Repeat the steps above with increasing number of spline points, e.g., 2, 3, 5, 10, and 20. Also try to estimate a primary production anomaly for all years by setting the number of spline points to 0. If the initial search for parameters doesn’t get it to move away from the starting point (flat line), then repeat (click yes when asked).
Finally, try doing a combined vulnerability search and anomaly search. First reset as earlier, then set the number of spline points to 3, and try searching for 2-8 vulnerability parameters, using the search by predator option.
Again, collect the results, examine the SS and AICc, and consider.
When fitting many parameters, you likely found that the vulnerability multiplier setting for whales was set very high. This makes whales perhaps double in abundance, and the fitting procedure likely chose this outcome because of the secondary effect whales have on other species through predation. This may or may not be real, and it certainly shows the importance of examining the outcome, and the value of having information at hand.
If you examine the fitting for cod and whiting, bear in mind that they are caught by the same fleet, and subject to the same change in effort. Cod is estimated to have a low vulnerability, whiting a high. Why? Both have decreased considerably according to the time series, but if you examine the Ecosim > Output > Ecosim group plots, you’ll notice that predation mortality is relatively more important for whiting than for cod (and fishing mortality, vice versa). So, as fishing effort for the trawlers increased significantly this would have the biggest impact on cod. To sustain this, the search routine chose a low vulnerability setting for cod. It would have less impact on whiting, but whiting decreased according to the time series. Therefore, the vulnerability setting had to be increased.
Sensible? Well, you need to have knowledge or at least a gut feeling of what has happened in your ecosystem. Data is good.
Other parameters
There are other factors that impact the fittings, but we do not (yet) have a built-in way of fitting to those. Here are some examples of parameters of importance,
Ecosim, Group info
The Ecosim Group info parameters are explained in a tutorial.
Feeding time adjustment rate: does a consumer change the amount of time it spends feeding when feeding conditions change? Most fish tend not to, so we usually set this parameter to 0 for all but marine mammals and birds, top predators, and juvenile/larval fish.
Density-dependent catchability: For schooling fish and for species showing range-contraction, the CPUE from fisheries are often misleading with regards to abundance. A corollary of this is that the catchability will change with their abundance (density-dependence!) If you have reason to think that may be the case for a group, you can use this parameter.
Handling time: The parameter QBmax/QB0 (for handling time) can be used to ensure a Holling type-II functional response. Lower the parameter to 2 or 3 to evaluate the effect.
Switching power: If a predator shows a Holling type-III functional response (switching), you can set this here. So, if the predator stops eating things as they become rare, and switches to feed relatively more of them as they become abundant, try using this parameter.
About AIC and AICc
[latex]AIC = nlog(RSS/n) + 2k + constant*n[/latex] [1]
where k is the number of parameters estimated and n is the number of observations being fitted to (i.e. n is the number of time series values, this being number of series used multiplied by the number of years for each). The constant*n can be ignored if n is the same (i.e. the observation data to be fitted to is the same) and we are comparing between alternative hypotheses.
So, using AIC to compare among alternative hypotheses (model parameterizations) in Ecosim, we need to calculate:
[latex]AIC = nlog(minSS (from Ecosim)/n) + 2k[/latex]
AICc is AIC with a second order correction for small sample sizes, to start with:
[latex]AIC_c = AIC + 2k(k-1)/(n-k-1)[/latex] where n is the number of observations
Since AICc converges to AIC as n gets large, AICc should be employed regardless of sample size (Burnham and Anderson, 2004).
n is the number of observations being fitted to (i.e. n is the number of time series values, this being number of series used multiplied by the number of years for each)
SS (sums of squared residuals ([Observed – Estimated]2) from Ecosim)
k is the number of parameters estimated. This is the number of distinct vulnerabilities plus number of primary production anomaly points (last year-first year or the number of spline points if these are used.
Quiz
- Venables, W. N., & Ripley, B. D. (2002). Modern applied statistics with S, 4th ed. New York, NY: Springer. ↵