bristlecone
Predator-Prey Dynamics: Snowshoe Hare and Lynx
Here we use the classic example of snowshoe hare and lynx predator-prey dynamics, to demonstrate the basic functions of Bristlecone. The dataset is a 90-year time-series of snowshoe hare and lynx pelts purchased by the Hudson's Bay Company of Canada. Data is in 1000s.
To get started, we first load and open the Bristlecone library in an F# script file (.fsx):
open Bristlecone // Opens Bristlecone core library and estimation engine
open Bristlecone.Language // Open the language for writing Bristlecone models
open Bristlecone.Time
Defining the ecological model
In Bristlecone, a single ecological model (representing a single hypothesis)
is defined through the ModelSystem type. A ModelSystem needs to include three key
components:
- Model equations. When working in continuous time, these are a system of Ordinary Differential Equations (ODEs).
- Parameters to be estimated. You must specify the starting bounds and constraints for every parameter included in the model equations.
- Likelihood function. The (negative log) likelihood function -logL represents the probability of observing the data given the parameter set. We use a negative log likelihood function, which is then minimised during optimisation.
In this example, we demonstrate using the Lotka–Volterra predator–prey model as the model system. For the -logL function we assume log-normal observation error for each of the two states, as we treat hare and lynx as dimensionless counts and assume multiplicative observation error; we model this with a lognormal likelihood. This is appropriate here because the fur‑return counts are strictly positive and vary over orders of magnitude. Any parameters passed into -logL functions are automatically estimated alongside model parameters; in this case these are sigma-hare and sigma-lynx.
Note. If we instead modelled density (e.g. individuals per km²), we would typically use a Gaussian or Student‑t likelihood on the raw scale, because taking logs of dimensional quantities is not physically meaningful.
[<Measure>] type prey = 1 // dimensionless counts
[<Measure>] type predator = 1 // dimensionless counts
// States
let H = state<prey> "hare"
let L = state<predator> "lynx"
// Likelihood parameters
let σH = parameter "σ[H]" Positive 0.001 0.100
let σL = parameter "σ[L]" Positive 0.001 0.100
let ``predator-prey`` =
// Parameters
let α = parameter "α" Positive 0.5<1/year> 1.5<1/year> // Maximum prey per‑capita growth rate
let β = parameter "β" Positive 0.01<1/(predator * year)> 0.05<1/(predator * year)> // Predation rate coefficient
let γ = parameter "γ" Positive 0.1<1/year> 1.0<1/year> // Predator natural mortality rate
let δ = parameter "δ" Positive 0.01<1/(prey * year)> 0.1<1/(prey * year)> // Predator growth efficiency
let ``dH/dt``: ModelExpression<prey / year> =
(P α - P β * State L) * This<prey>
let ``dL/dt``: ModelExpression<predator / year> =
(-P γ + P δ * State H) * This<predator>
// Likelihood
let NLL =
ModelLibrary.NegLogLikelihood.LogNormal (Require.state H) σH +
ModelLibrary.NegLogLikelihood.LogNormal (Require.state L) σL
Model.empty
|> Model.addRateEquation H ``dH/dt``
|> Model.addRateEquation L ``dL/dt``
|> Model.estimateParameter α
|> Model.estimateParameter β
|> Model.estimateParameter δ
|> Model.estimateParameter γ
|> Model.useLikelihoodFunction NLL
|> Model.compile
Setting up the Bristlecone Engine
A bristlecone engine provides a fixed setup for estimating parameters from data. Use the same engine for all model fits within a single study. This engine uses a gradident descent method (Nelder Mead simplex), and a basic Runge-Kutta 4 integration method provided by MathNet Numerics.
let engine: EstimationEngine.EstimationEngine<DatingMethods.Annual,int<year>,year,1> =
Bristlecone.mkContinuous ()
|> Bristlecone.withBristleconeOptimiser
|> Bristlecone.withConditioning Conditioning.NoConditioning
|> Bristlecone.withSeed 1500 // We are setting a seed for this example - see below
|> Bristlecone.withTimeConversion DateMode.Conversion.Annual.toYears
Note. We have set a seed for random number generation for this worked example. This ensures that the results are the same each time this documentation is generated.
Does it all work? Testing the engine and model
Before being confident in the ability of our estimation engine to be able to arrive at the correct solution, we must run a full test of the model and estimation engine.
Bristlecone includes the Bristlecone.testModel function, which
we use here. Given a model system and estimation engine, the function
generates a random parameter set (θ) many times; for each θ, the
'true' time-series are generated. The test result indicates the effectiveness
of the configuration at estimating θ given the auto-generated
time-series. If there is divergence, there is likely an
issue with your model or the Bristlecone configuration.
Bristlecone includes many settings to configure the test
procedure. A simple test configuration is set as Test.defaultSettings,
but here we will configure some additional settings:
let testSettings =
Test.createWithTimeMode DateMode.annualDateMode (Yearly 1<year>) (DatingMethods.Annual 1970<year>)
|> Test.seriesLength 15
|> Test.t1 (Require.state H) 50.
|> Test.t1 (Require.state L) 75.
|> Test.withObservationError (Require.state H) (Test.Error.logNormal σH)
|> Test.withObservationError (Require.state L) (Test.Error.logNormal σL)
|> Test.rule (Require.state H) (Test.GenerationRules.between 10. 10000.)
|> Test.rule (Require.state L) (Test.GenerationRules.between 10. 10000.)
In our TestSettings, we have specified the initial time point (t = 0)
for both modelled time-series. We have also added noise around
each generated time-series, and specified that each time-series
should be 30 years in length.
With these test settings, we can now run the test.
let endCondition = Optimisation.EndConditions.Profiles.mcmc 5000<iteration> engine.LogTo
let testResult = ``predator-prey`` |> Bristlecone.tryTestModel engine endCondition testSettings
(forceOk testResult).Parameters
In the example using this seed, the ecological parameter set under test is:
1.1070 (α), 0.03299 (β), 0.7930 (γ), 0.01487 (δ) [ error: 0.0706 (σH), 0.0868 (σL) ] Here, the sigma values represent the standard deviation of the noise around hare and lynx respectively.
From the test results, the key questions is whether the ecological dynamics were recovered. To identify this we can look at these factors:
- Was the estimated -logL near the true -logL?
- Were all of the true parameters recovered (within a tolerance)?
- Do the estimated and observed trajectories look similar?
- Are there trends in the error structure over time?
First, the true -logL was 65.8821, while the estimated minimum was 66.5655, which are close given that observation error was added.
Second, the parameters were recovered to reasonable levels, again given that observation error was added to a short 15-point time-series:
val it: Test.ParameterTestResult list =
[{ Identifier = "α"
RealValue = 1.149500728
EstimatedValue = 1.107033385 }; { Identifier = "β"
RealValue = 0.0329926312
EstimatedValue = 0.03244952328 };
{ Identifier = "γ"
RealValue = 0.7930286527
EstimatedValue = 0.8282329879 }; { Identifier = "δ"
RealValue = 0.0148735689
EstimatedValue = 0.0160377167 };
{ Identifier = "σ[H]"
RealValue = 0.07067647576
EstimatedValue = 0.05908173058 }; { Identifier = "σ[L]"
RealValue = 0.08681011945
EstimatedValue = 0.09076999054 }]
Third, we can look at the predicted versus observed time-series:
The fit is relatively tight. The test result type also includes
.ErrorStructure, which shows the per-observation root squared
error in the units state^2, where state here is the count of hare
/ lynx.
We can also examine raw parameter traces to visualise the contours
of the parameter space. Note that the traces are
organised from newest to oldest. The .Trace is organised
into a list, where each item is a stage or sub-stage of the
optimisation method applied. Here, this includes first a pre-tuning phase,
then a tuning phase, each heating level, each annealing level,
and finally a 'clean' MCMC chain.