Using the tdsa package for time-dependent sensitivity analysis: a basic example

The tdsa package is designed to automate most of the steps needed to perform time-dependent sensitivity analysis (TDSA). Please refer to Ng et al. (in press), Ng et al. (in review), or the vignette “What is time-dependent sensitivity analysis (TDSA)?” for more details on TDSA. In this vignette, we demonstrate how to use the package with a simple example; a more complex example will be presented in a later vignette. Our focus will be on the input arguments that the user needs to provide, as well as the structure of the output.

Translocating individuals into a sink habitat that is being restored

Background

Conservation translocation refers to the human-mediated movement of a threatened species from one area to another, in order to restore or expand the species’ range, or to bolster a small or declining population (Armstrong & Seddon, 2007). While this is often done to reduce the threat of extinction, sometimes the ecosystem service rendered by the species may also provide a strong motivation. For example, although there is still significant debate, the recovery of woody plant species like aspens and cottonwoods at the Yellowstone National Park has been attributed to the biological control of the elk population by wolves reintroduced in the 1990s (Ripple & Beschta, 2012).

Timing is expected to be important for success of a translocation effort. In particular, attempts should only be made after the original threats that led to the species decline or extirpation have been sufficiently reduced (IUCN, 2013). Motivated by this, we created a simple hypothetical model to demonstrate how TDSA can be used to highlight the importance of translocation timing (Ng et al., in press). This is also the example in the help page of the function state_sens.

Consider an organism in a sink habitat, where the per-capita loss rate (mortality and emigration combined) exceeds the per-capita unregulated birth rate, so the population is only maintained through immigration. However, the mortality rate is expected to decrease over time due to ongoing habitat restoration or threat reduction efforts, so the population should eventually become self-sustaining. The population dynamics is hence given by \[\frac{dy}{dt} = by(t)(1 - ay(t)) - \mu(t)y(t) + \sigma,\] where \(y(t)\) is the population at time \(t\), \(b\) the unregulated per-capita birth rate, \(a\) the coefficient for reproductive competition, \(\mu(t)\) the time-varying per-capita loss rate, and \(\sigma\) the immigration rate. We assume that \(\mu(t)\) starts off above \(b\) (so it is a sink habitat), but decreases as a sigmoid and eventually falls below \(b\) (so the population becomes self-sustaining).

The organism provides an important ecosystem service. Over a management period from \(t_0\) to \(t_1\), we ascribe a value to the organism (the reward function) \[J=\int_{t_0}^{t_1} \underbrace{w y(t)}_{\substack{\text{reward}\\\text{integrand}}} dt + \underbrace{v y(t_1)}_{\substack{\text{terminal}\\ \text{payoff}}}.\] Here, it is assumed that each individual provides the service at a rate \(w\), so the integral gives the total amount of service accumulated over the period. However, we also want to ascribe value to maintaining a large population at the end of the management period, so the second term corresponds to a terminal payoff where \(v\) is the ascribed value per individual.

Say we want to translocate individuals to the habitat to speed up the population recovery and increase the reward \(J\). What is the best time to do so in order to maximise the increase in the reward? As early as possible? Or only when the loss rate has become low enough that the population can sustain itself? A one-off translocation causes a small, sudden increase in the population size, so it is useful to look at the time-dependent state sensitivity. Alternatively, we can interpret the translocation as a brief spike in the immigration rate \(\sigma\), so we can also look at the time-dependent parameter sensitivity of \(\sigma\).

Installing and loading the tdsa package

The tdsa package is already on CRAN, so it can be easily installed using the command install.packages("tdsa"). For the rest of this vignette, we assume that the package has already been installed.

Before proceeding further, we load the tdsa package.

library(tdsa)

Preparing the input arguments

To perform time-dependent sensitivity analysis using the function state_sens, we need to prepare the input arguments:

dynamic_fn, a function of the form function(t, y, parms, ...) that returns the RHS of the dynamic equations. For consistency with the syntax of the deSolve package (which many users might already be familiar with), the output must be a list, the first element a numeric vector of length equal to the number of state variables (which is just one in this example).
parms, an object used to specify input parameters for dynamic_fn. (We will discuss other ways to specify input parameters in the later vignette.)
reward_fn, a function of the form function(t, y, ...) that returns the integrand in the reward function.
terminal_fn, a function of the form function(y) that returns the terminal payoff.
y_0, a numeric vector used to specify the initial conditions \(y(t_0)\).
times, a numeric vector used to specify the discretised time steps between \(t_0\) and \(t_1\) (inclusive) used in the simulation.

# Parameter values for the dynamic equations.
parms = list(
  b = 1,                                          # Per-capita birth rate.
  a = 0.1,                                        # Competition coefficient.
  mu = function(t){0.5 + 1/(1 + exp((t-10)/2))},  # Per-capita loss rate.
  sigma = 0.2                                     # Immigration rate.
)

# Function that returns the dynamic equations.
dynamic_fn = function(t, y, parms){
  b = parms[["b"]]
  a = parms[["a"]]
  sigma = parms[["sigma"]]
  mu = parms[["mu"]](t)
  
  dy = b*y*(1- a*y) - mu*y + sigma
  return( list(dy) )
}

# Initial conditions.
y_0 = 0.37  # Approximate steady-state population before restoration efforts.

# Function that returns the reward integrand.
reward_fn = function(t, y){
  w = 1  # Per-capita rate at which the ecosystem service is provided.
  return( w * y )
}

# Function that returns the terminal payoff.
terminal_fn = function(y){
  v = 1.74  # Ascribed value per individual at the end of the period.
  return( v * y )
}

# Time steps over management period. We discretise it into 1001 time steps
# (so the step size is 0.03).
times = seq(0, 30, length.out=1001)

Calculating time-dependent state sensitivities

We first calculate time-dependent state sensitivities using the function state_sens. Since this is a continuous-time model, we choose model_type = "continuous". Since this is a simple model that should only takes seconds to run, there is no need to show progress indicators in the console, so we set verbose = FALSE.

state_sens_out = state_sens(
  model_type = "continuous",
  dynamic_fn = dynamic_fn,
  parms = parms,
  reward_fn = reward_fn,
  terminal_fn = terminal_fn,
  y_0 = y_0,
  times = times,
  verbose = FALSE
)

The output is a list that contains some of the input arguments such as dynamic_fn and times that will be needed to calculate parameter sensitivities later on, as well as two matrices. The first matrix, called state, gives the state variable \(y\) at each time step in times. Users of the deSolve package should find this familiar, except that we have removed the first column containing the time steps. The second matrix, called tdss, gives the state sensitivity \(\lambda\) also at each time step in times.

str(state_sens_out)

## List of 7
##  $ model_type        : chr "continuous"
##  $ dynamic_fn        :function (t, y, parms)  
##   ..- attr(*, "srcref")= 'srcref' int [1:8] 10 14 18 1 14 1 10 18
##   .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x0000000014914140> 
##  $ parms             :List of 4
##   ..$ b    : num 1
##   ..$ a    : num 0.1
##   ..$ mu   :function (t)  
##   .. ..- attr(*, "srcref")= 'srcref' int [1:8] 5 8 5 47 8 47 5 5
##   .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x0000000014914140> 
##   ..$ sigma: num 0.2
##  $ dynamic_fn_arglist: list()
##  $ times             : num [1:1001] 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 ...
##  $ state             : num [1:1001, 1] 0.37 0.37 0.37 0.37 0.37 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "1"
##  $ tdss              : num [1:1001, 1] 1.89 1.89 1.89 1.89 1.9 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr "1"

In the left panel below, we plot the per-capita unregulated birth and loss rates, from parms. They intersect at around \(t=10\). The middle panel shows the population size, from the abovementioned state matrix. The right panel shows the time-dependent state sensitivity, from the tdss matrix. We see that the sensitivity also peaks at around \(t=10\), so translocation is most effective when the population has just become self-sustaining, an intuitive result.

# Set graphical parameters.
par(mfrow=c(1,3), cex=1)
par(mar=c(3.2,3.2,2,2), mgp=c(2,0.7,0), cex.lab=1.2)

# Plot the per-capita unregulated birth and loss rates.
plot(times, parms[["mu"]](times), type="l", lwd=2,
     xlab="Time (year)", ylab="Demographic rate (/year)")
abline(h=parms[["b"]], col="red", lwd=2)
legend("topright", col=c("red", "black"), lwd=2, bty="n",
       legend=c("Birth rate", "Loss rate"))

# Plot the population size.
plot(times, state_sens_out[["state"]][,1], type="l", lwd=2,
     xlab="Time (year)", ylab="Population size y")

# Plot the time-dependent state sensitivity. Peaks at around t=10, which is
# roughly when mu and b intersects, so the population has just become
# self-sustaining.
plot(times, state_sens_out[["tdss"]][,1], type="l", lwd=2,
     xlab="Time (year)", ylab="State sensitivity of y")

Calculating time-dependent parameter sensitivities

Once we have calculated the state sensitivities, calculating the parameter sensitivities is easy—all we have to do is to use the output of the function state_sens as the input argument of the function parm_sens. This is because in addition to the state sensitivities, we have also included in the output of state_sens the original dynamic_fn that implicitly relates the model parameters to the state variables.

parm_sens_out = parm_sens(state_sens_out = state_sens_out,
                          verbose = FALSE)

The output of parm_sens is a list containing two elements. The first element is times, while the second element tdps is an object that gives the time-dependent parameter sensitivities. The structure of tdps depends on the structure of parms.

If parms is a numeric object (vector, matrix or array), or a function that returns a numeric object (i.e., a time-varying parameter), then tdps is an array with one more index than the object, with the first index indicating the time step. So a parameter vector (in parms) generates a matrix in the output, in which the first column is time and the other columns are the sensitivity to each parameter in the parameter vector at that time. Similarly a matrix becomes a 3-index array, etc.
If parms is a list containing any combination of the above, then tdps is a list with the above rule applied element-wise.

In this example, parms is a list containing single-number parameters b, a, sigma and mu (since mu is a function that returns a single number), so tdps is a list containing a one-column matrix for each parameter.

str(parm_sens_out)

## List of 2
##  $ times: num [1:1001] 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24 0.27 ...
##  $ tdps :List of 4
##   ..$ b    : num [1:1001, 1] 0.672 0.673 0.674 0.675 0.676 ...
##   ..$ a    : num [1:1001, 1] -0.258 -0.259 -0.259 -0.26 -0.26 ...
##   ..$ mu   : num [1:1001, 1] -0.698 -0.699 -0.7 -0.701 -0.702 ...
##   ..$ sigma: num [1:1001, 1] 1.89 1.89 1.89 1.89 1.9 ...

Below, we plot the parameter sensitivity of \(\sigma\). We find that it is identical to the state sensitivity, which is not surprising because from the dynamic equations, we know that a brief spike in \(\sigma\) should lead to a sudden increase in the population size and hence the same increase in \(J\).

# Set graphical parameters.
par(mar=c(3.2,3.2,2,2), mgp=c(2,0.7,0), cex.lab=1.2)

# Plot the parameter sensitivity of sigma.
plot(times, parm_sens_out[["tdps"]][["sigma"]][,1], type="l", lwd=2,
     xlab="Time (year)", ylab="Param. sensitivity of sigma")