Units and priors

2021-09-24

In this vignette, we will discuss two important points which are related to each other: the parameters physical units and how to choose parameter priors. Those two are related since the biological meaning of the priors necessarily depends on the units used for parameters.

Parameter units

Some parameters of the network models are unitless, but those which aren't have units which depend on the units used for the data.

This is an overview of the different types of parameters that are used in a network model:

Parameter name Description Unit Domain
upsilon_<src>_to_<dst> Uptake rate from src to dst Depends on the data \(\left[0;+\infty\right[\)
lambda_<src> Loss rate from src Depends on the data \(\left[0;+\infty\right[\)
portion.act_<src> Active portion of src Unitless \(\left[0;1\right]\)
eta Default: coefficient of variation (tracer proportion) Unitless \(\left[0;+\infty\right[\)
zeta Default: coefficient of variation (biomass) Unitless \(\left[0;+\infty\right[\)

Uptake rates

Uptake rates (upsilon_<src>_to_<dst>) are expressed in proportion of compartment material transferred per unit of time. For example: 0.05 \(\textrm{day}^{-1}\) means that 5% of src is transferred to dst per day.

The actual time unit corresponds to the time unit used in the observation data. For example, if the observation data give sampling times as hours since the beginning of the experiment, then the uptake rates will be in \(\textrm{hour}^{-1}\).

Loss rates

Loss rates (lambda_<src>) are also expressed in proportion of compartment material transferred (lost) per unit of time. For example: 0.1 \(\textrm{week}^{-1}\) means that 10% of src is lost and exits the network system per week.

Again, the actual time unit corresponds to the time unit used in the observation data. It is the same as for uptake rates.

Active portions

An active portion is the proportion of a compartment initial standing stock that is being involved in flows in the network, and must be comprised between 0 and 1. For example, an active portion of 0.25 means that only 25% of the initial compartment standing stock is used when calculating uptakes by other compartments and losses, while 75% of the initial compartment standing stock behaves as a refractory portion which is not involved in the network.

Active portions are useful to model compartments for which only some of the content is involved in exchanges at the time scale of the experiment (e.g. algal and bacterial layers growing on dead leaves in a stream). They can explain why consumers feeding selectively on this active portion might exhibit higher isotopic enrichments than the apparent enrichment of the consumed compartment.

Coefficients of variations

By default, \(\eta\) and \(\zeta\) are the coefficients of variation for the distribution of observed isotopic proportions and compartment sizes, respectively, around their predicted mean values. As coefficients of variation, they are the ratio between the standard deviation and the mean of the distributions they describe, and are thus unitless.

For \(\eta\), several parameterizations are actually implemented in isotracer and can be set with set_prop_family():

Choosing priors

Once we understand what the physical units of each parameter are, we can make sensible choices for priors.

For uptake and loss rates, half-Cauchy priors are reasonable choices since they allow for any positive value but gives more weight to smaller values and can be scaled by their median. For example, if we consider a stream ecosystem for which the time data was given in days, allowing for 25% of a compartment biomass to be renewed every day is already pretty generous, so a half-Cauchy with scale 0.25 is a slightly informative but reasonable prior.

For active portions, beta priors can be used if the user has any preliminary knowledge about expected active portions, but since the interval is bounded on \(\left[0;1\right]\) a uniform prior is proper and acceptable.

For coefficients of variation, half-Cauchy priors are again a reasonable choice. For example, a half-Cauchy with median 0.5 is not very informative but gives more weight to coefficients of variation below 50%.

Setting priors

For now, a few different types of priors can be specified by the user when building a network model. Let's illustrate this with our usual example of an aquarium with NH4, algae and Daphnia:

library(isotracer)
library(tidyverse)
mod <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia -> NH4")

What are the parameters of the model?

params(mod)
## [1] "eta"                      "lambda_algae"             "lambda_daphnia"          
## [4] "lambda_NH4"               "upsilon_algae_to_daphnia" "upsilon_daphnia_to_NH4"  
## [7] "upsilon_NH4_to_algae"     "zeta"

What are the default priors used?

priors(mod)
## # A tibble: 8 × 2
##   in_model                 prior                
##   <chr>                    <list>               
## 1 eta                      <hcauchy (scale=0.1)>
## 2 lambda_algae             <hcauchy (scale=0.1)>
## 3 lambda_daphnia           <hcauchy (scale=0.1)>
## 4 lambda_NH4               <hcauchy (scale=0.1)>
## 5 upsilon_algae_to_daphnia <hcauchy (scale=0.1)>
## 6 upsilon_daphnia_to_NH4   <hcauchy (scale=0.1)>
## 7 upsilon_NH4_to_algae     <hcauchy (scale=0.1)>
## 8 zeta                     <hcauchy (scale=0.1)>

By default, all parameters are given a half-Cauchy prior defined on \(\left[0;+\infty\right[\), with a scale of 0.1 (for a positive half-Cauchy prior, the scale is the median).

The priors available in isotracer are:

Priors of a network model are set using the set_prior() function:

mod <- mod %>% set_prior(uniform(0, 10), "zeta")
## Prior modified for parameter(s): 
##   - zeta
priors(mod)
## # A tibble: 8 × 2
##   in_model                 prior                   
##   <chr>                    <list>                  
## 1 eta                      <hcauchy (scale=0.1)>   
## 2 lambda_algae             <hcauchy (scale=0.1)>   
## 3 lambda_daphnia           <hcauchy (scale=0.1)>   
## 4 lambda_NH4               <hcauchy (scale=0.1)>   
## 5 upsilon_algae_to_daphnia <hcauchy (scale=0.1)>   
## 6 upsilon_daphnia_to_NH4   <hcauchy (scale=0.1)>   
## 7 upsilon_NH4_to_algae     <hcauchy (scale=0.1)>   
## 8 zeta                     <uniform (min=0,max=10)>

Note that by default set_prior() assigns priors to all parameters that match the name provided:

mod <- mod %>% set_prior(scaled_beta(4, 8, 5), "lambda")
## Prior modified for parameter(s): 
##   - lambda_algae
##   - lambda_daphnia
##   - lambda_NH4
priors(mod)
## # A tibble: 8 × 2
##   in_model                 prior                                 
##   <chr>                    <list>                                
## 1 eta                      <hcauchy (scale=0.1)>                 
## 2 lambda_algae             <scaled_beta (alpha=4,beta=8,scale=5)>
## 3 lambda_daphnia           <scaled_beta (alpha=4,beta=8,scale=5)>
## 4 lambda_NH4               <scaled_beta (alpha=4,beta=8,scale=5)>
## 5 upsilon_algae_to_daphnia <hcauchy (scale=0.1)>                 
## 6 upsilon_daphnia_to_NH4   <hcauchy (scale=0.1)>                 
## 7 upsilon_NH4_to_algae     <hcauchy (scale=0.1)>                 
## 8 zeta                     <uniform (min=0,max=10)>

If you want to avoid using partial matching, you can use the argument use_regexp = FALSE:

mod <- mod %>% set_prior(constant(0.2), "eta", use_regexp = FALSE)
## Prior modified for parameter(s): 
##   - eta
priors(mod)
## # A tibble: 8 × 2
##   in_model                 prior                                 
##   <chr>                    <list>                                
## 1 eta                      <constant (value=0.2)>                
## 2 lambda_algae             <scaled_beta (alpha=4,beta=8,scale=5)>
## 3 lambda_daphnia           <scaled_beta (alpha=4,beta=8,scale=5)>
## 4 lambda_NH4               <scaled_beta (alpha=4,beta=8,scale=5)>
## 5 upsilon_algae_to_daphnia <hcauchy (scale=0.1)>                 
## 6 upsilon_daphnia_to_NH4   <hcauchy (scale=0.1)>                 
## 7 upsilon_NH4_to_algae     <hcauchy (scale=0.1)>                 
## 8 zeta                     <uniform (min=0,max=10)>

Notice in the last example how eta was affected while zeta was not.