Header menu logo bristlecone

Composing model parts together

Bristlecone has been designed to enable swift model composition, whereby smaller models may be composed into larger systems of models.

To get started, open the Bristlecone libraries.

open Bristlecone // Opens Bristlecone core library and estimation engine
open Bristlecone.Language // Open the language for writing Bristlecone models

Nested model systems

One key use of compositional modelling is to understand which (ecological) theories best represent one or many particular time-series data. A nested model is one where one or more processes may be substituted for different theoretical forms and competed to determine which combination of model components best represents the given data.

Using Bristlecone, we can simply define a nested model as a Bristlecone model but where the definition takes some additional model components as function parameters. For example, take the following model:

let someModel someFn =
    Parameter "K" + someFn This - Environment "T" * Parameter "J"

In the above expression, someFn is a function of signature ModelExpression -> ModelExpression. Looking at the equation, you can see that someFn takes one argument, the model expression This, which represents the current state of the model equation.

We can take this concept further to define a nestable model, in which a number of components may be swapped for alternative mathematical forms based on (ecological) theory. Here, we will apply a simple model of non-linear plant growth. To do this, we first define a base model as a ModelSystem that takes a number of interchangable components:

let baseModel growthLimit lossRate =
    Model.empty
    |> Model.addEquation "m" (Parameter "r" * (growthLimit This) - lossRate This)
    |> Model.estimateParameter "r" noConstraints 0.01 1.00
    |> Model.useLikelihoodFunction (ModelLibrary.Likelihood.sumOfSquares [ "x" ])

Importantly, unlike when normally building models, we do not call Model.compile at the end. This step will be done later.

We can then design the growth limit and loss rate components that will fit into this model. These are made using the modelComponent and subComponent functions as follows:

let growthLimit =
    modelComponent
        "Limitation to growth"
        [

          subComponent "Linear" (fun _ -> Constant 1.)

          subComponent "Monomolecular" (fun m -> Parameter "K" - m)
          |> estimateParameter "K" notNegative 10.0 50.0

          subComponent "Gompertz" (fun m -> m * log (Parameter "K" / m))
          |> estimateParameter "K" notNegative 10.0 50.0

          subComponent "Logisitc (three parameter)" (fun m -> m * (Constant 1. - (m / Parameter "K")))
          |> estimateParameter "K" notNegative 10.0 50.0 ]

If a component requires additional parameters over the base model, these may be added by piping into the estimateParameter function as above.

We can similarly define the other required component for the base model - the loss rate:

let lossRate =
    modelComponent
        "Biomass loss rate"
        [

          subComponent "None" (fun _ -> Constant 0.)

          subComponent "Density-dependent" (fun m -> m * Parameter "alpha")
          |> estimateParameter "alpha" notNegative 0.001 0.10 ]

Once the components are set up, we can compile the nested model by making all possible combinations of the model components. We do this using the Hypotheses module, which will compile a list of model hypotheses for us to test.

let hypotheses =
    baseModel
    |> Hypotheses.createFromComponent growthLimit
    |> Hypotheses.useAnother lossRate
    |> Hypotheses.compile

The resultant constructed hypotheses are shown in the following printout:

[Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Linear" }; { Component = "Biomass loss rate"
                                      Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Monomolecular" }; { Component = "Biomass loss rate"
                                             Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Gompertz" }; { Component = "Biomass loss rate"
                                        Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Logisitc (three parameter)" };
     { Component = "Biomass loss rate"
       Implementation = "None" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "alpha",
              Parameter (PositiveOnly, Transform, NotEstimated (0.001, 0.1)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Linear" }; { Component = "Biomass loss rate"
                                      Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "alpha",
              Parameter (PositiveOnly, Transform, NotEstimated (0.001, 0.1)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Monomolecular" };
     { Component = "Biomass loss rate"
       Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "alpha",
              Parameter (PositiveOnly, Transform, NotEstimated (0.001, 0.1)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Gompertz" }; { Component = "Biomass loss rate"
                                        Implementation = "Density-dependent" }]);
 Hypothesis
   ({ Parameters =
       Pool
         (map
            [(ShortCode "K",
              Parameter (PositiveOnly, Transform, NotEstimated (10.0, 50.0)));
             (ShortCode "alpha",
              Parameter (PositiveOnly, Transform, NotEstimated (0.001, 0.1)));
             (ShortCode "r",
              Parameter (Unconstrained, Transform, NotEstimated (0.01, 1.0)))])
      Equations = map [(ShortCode "m", <fun:Invoke@3682-2>)]
      Measures = map []
      NegLogLikelihood = <fun:baseModel@5> },
    [{ Component = "Limitation to growth"
       Implementation = "Logisitc (three parameter)" };
     { Component = "Biomass loss rate"
       Implementation = "Density-dependent" }])]

For further analysis, you may choose to use the orchestration functionality to run many hypotheses and replicates (multiple fits per hypothesis) in parallel.

Multiple items
module Bristlecone from Bristlecone
<namespacedoc><summary>The core library of Bristlecone, containing model-fitting functions.</summary></namespacedoc>
Main functionality of Bristlecone, including functions to scaffold `ModelSystem`s and for model-fitting (tests and real fits).


--------------------
namespace Bristlecone
module Language from Bristlecone
<summary> An F# Domain Specific Language (DSL) for scripting with Bristlecone. </summary>
val someModel: someFn: (ModelExpression -> ModelExpression) -> ModelExpression
val someFn: (ModelExpression -> ModelExpression)
Multiple items
union case ModelExpression.Parameter: string -> ModelExpression

--------------------
module Parameter from Bristlecone
union case ModelExpression.This: ModelExpression
union case ModelExpression.Environment: string -> ModelExpression
val baseModel: growthLimit: (ModelExpression -> ModelExpression) -> lossRate: (ModelExpression -> ModelExpression) -> ModelBuilder.ModelBuilder
val growthLimit: (ModelExpression -> ModelExpression)
val lossRate: (ModelExpression -> ModelExpression)
module Model from Bristlecone.Language
<summary> Terms for scaffolding a model system for use with Bristlecone. </summary>
val empty: ModelBuilder.ModelBuilder
val addEquation: name: string -> eq: ModelExpression -> builder: ModelBuilder.ModelBuilder -> ModelBuilder.ModelBuilder
val estimateParameter: name: string -> constraintMode: Parameter.Constraint -> lower: float -> upper: float -> builder: ModelBuilder.ModelBuilder -> ModelBuilder.ModelBuilder
val noConstraints: Parameter.Constraint
val useLikelihoodFunction: likelihoodFn: ModelSystem.LikelihoodFn -> builder: ModelBuilder.ModelBuilder -> ModelBuilder.ModelBuilder
namespace Bristlecone.ModelLibrary
module Likelihood from Bristlecone.ModelLibrary
<summary>Likelihood functions to represent a variety of distributions and data types.</summary>
<namespacedoc><summary>Pre-built model parts for use in Bristlecone</summary></namespacedoc>
val sumOfSquares: keys: string list -> ModelSystem.ParameterValueAccessor -> data: CodedMap<ModelSystem.PredictedSeries> -> float
<summary> Residual sum of squares. Provides a simple metric of distance between observed data and model predictions. </summary>
val growthLimit: string * (string * PluggableComponent<ModelExpression>) list
val modelComponent: name: 'a -> list: 'b -> 'a * 'b
val subComponent: name: 'a -> expression: ('b -> ModelExpression) -> 'a * PluggableComponent<'b>
<summary> Creates a nested component that can be inserted into a base model. </summary>
union case ModelExpression.Constant: float -> ModelExpression
val m: ModelExpression
val estimateParameter: name: string -> constraintMode: Parameter.Constraint -> lower: float -> upper: float -> 'a * PluggableComponent<'b> -> 'a * PluggableComponent<'b>
val notNegative: Parameter.Constraint
val log: value: 'T -> 'T (requires member Log)
val lossRate: string * (string * PluggableComponent<ModelExpression>) list
val hypotheses: Hypotheses.Hypothesis list
module Hypotheses from Bristlecone.Language
<summary>Types to represent a hypothesis, given that a hypothesis is a model system that contains some alternate formulations of certain components.</summary>
val createFromComponent: string * (string * PluggableComponent<'a>) list -> builder: (('a -> ModelExpression) -> 'b) -> Writer.Writer<'b,(Hypotheses.ComponentName * CodedMap<Parameter.Parameter>)> list
<summary>Implement a component where a model system requires one. A component is a part of a model that may be varied, for example between competing hypotheses.</summary>
<param name="comp">A tuple representing a component scaffolded with the `modelComponent` and `subComponent` functions.</param>
<param name="builder">A builder started with `createFromComponent`</param>
<typeparam name="'a"></typeparam>
<typeparam name="'b"></typeparam>
<returns>A builder to add further components or compile with `Hypothesis.compile`</returns>
val useAnother: string * (string * PluggableComponent<'a>) list -> builder: Writer.Writer<(('a -> ModelExpression) -> 'b),(Hypotheses.ComponentName * CodedMap<Parameter.Parameter>)> list -> Writer.Writer<'b,(Hypotheses.ComponentName * CodedMap<Parameter.Parameter>)> list
<summary>Add a second or further component to a model system where one is still required.</summary>
<param name="comp">A tuple representing a component scaffolded with the `modelComponent` and `subComponent` functions.</param>
<param name="builder">A builder started with `createFromComponent`</param>
<typeparam name="'a"></typeparam>
<typeparam name="'b"></typeparam>
<returns>A builder to add further components or compile with `Hypothesis.compile`</returns>
val compile: builder: Writer.Writer<ModelBuilder.ModelBuilder,(Hypotheses.ComponentName * CodedMap<Parameter.Parameter>)> list -> Hypotheses.Hypothesis list
<summary>Compiles a suite of competing model hypotheses based on the given components. The compilation includes only the required parameters in each model hypothesis, and combines all labels into a single model identifier.</summary>
<param name="builder">A builder started with `createFromComponent`</param>
<returns>A list of compiled hypotheses for this model system and specified components.</returns>

Type something to start searching.