The trophic package has several different uses. These are some examples of possible use cases. Please feel free to contribute your own additions (see here for details).

example: predicting biomass of a high-order predator with known food web

This is a simple example where we have a known (or assumed) food web and some estimates of trophic transfer efficiencies. The food web additionally contains information on the relative dominance of consumers over potential competitors. This information is captured in the objects food_web (a lower triangular matrix of zeros and ones), efficiency_mean (a series of proportions for trophic transfer efficiencies), and dominance_matrix (a matrix of relative values for each consumer). In food webs in trophic, we assume that “rows eat columns”, that is, a value in a given cell indicates that the node in that row eats the node in that column.

Example fixed food web
decapod.crustaceans small.fish carp reptiles mammals birds large.fish
decapod.crustaceans 0 0 0 0 0 0 0
small.fish 0 0 0 0 0 0 0
carp 0 0 0 0 0 0 0
reptiles 0 1 0 0 0 0 0
mammals 1 0 0 0 0 0 0
birds 1 1 0 0 0 0 0
large.bodied.native.fish 1 1 0 0 1 1 0
Example efficiency matrix
decapod.crust small.fish carp reptiles mammals birds large.bodied.fish
decapod.crust 0.1 0.15 0.12 0.09 0.05 0.07 0.04
small.fish 0.1 0.15 0.12 0.09 0.05 0.07 0.04
carp 0.1 0.15 0.12 0.09 0.05 0.07 0.04
reptiles 0.1 0.15 0.12 0.09 0.05 0.07 0.04
mammals 0.1 0.15 0.12 0.09 0.05 0.07 0.04
birds 0.1 0.15 0.12 0.09 0.05 0.07 0.04
large.bodied.fish 0.1 0.15 0.12 0.09 0.05 0.07 0.04
Example dominance matrix
decapod.crustaceans small.fish carp reptiles mammals birds large.fish
decapod.crustaceans 0 0 0 0 0 0 0
small.fish 0 0 0 0 0 0 0
carp 0 0 0 0 0 0 0
reptiles 0 2 0 0 0 0 0
mammals 1 0 0 0 0 0 0
birds 1 4 0 0 0 0 0
large.bodied.native.fish 8 4 0 0 12 4 0

If transfer efficiencies are unknown, they can be assumed to be 10% (0.1) by calling build_trophic_dynamics() with efficiency_matrix = NULL. Similarly, if dominance values are unknown, they can be assumed to be equal (all consumers will get an equal share of a shared resource) by setting dominance_matrix = NULL in the build_trophic_dynamics() function.

We additionally have estimates of primary production (production_mean and production_sd) that we use to initialise the food web. These values are required for all models in trophic. Estimates of primary production can be collated into a trophic object with the build_primary_producers() function, which can take named vectors of means and standard deviations. Alternatively, if production_mean and production_sd are not named, the model will assume that production values correspond to the first columns in the food web (one column for each value in production_mean and production_sd).

# load trophic package
library(trophic)

# prepare food web, transfer efficiencies, and dominance matrix
example_fw <- build_food_web(interaction_matrix = food_web)
example_eps <- build_efficiency_matrix(efficiency_mean = efficiency_mean,
                                       efficiency_sd = 0.05)
example_doms <- build_dominance_matrix(dominance = dominance_matrix)

# build primary producers object
example_pp <- build_primary_producers(production_mean = production_mean,
                                      production_sd = production_sd)

# collate trophic_dynamics object
example_trophic <- build_trophic_dynamics(food_web = example_fw,
                                          efficiency_matrix = example_eps,
                                          dominance_matrix = example_doms)

The collated trophic_dynamics and primary_producers objects can be used to predict production at any node in the food web with the estimate_production() function. This function takes two additional arguments: stochastic and nsim. The stochastic argument specifies which components of the projection should be estimated with uncertainty. These can be any combination of food_web, efficiency, and primary_production. If stochastic = NULL or stochastic = c("random_text") then the production estimates will be simulated without noise (based on mean values in trophic_dynamics and primary_producers).

# predict production of all nodes
predicted_production <- estimate_production(trophic_dynamics = example_trophic,
                                            primary_producers = example_pp,
                                            stochastic = c("efficiency",
                                                           "primary_production"))

Estimated production values can be converted to biomass using the estimate_biomass() function, with an appropriate pb_ratio object. A pb_ratio object specifies the ratio of production to biomass, i.e., how much production can a node sustain given its biomass. We use the inverse of this value to predict biomass from estimated production, scaled by a conversion factor of 10 to link carbon production (\(g C day^{-1}\)) to total biomass production (\(g biomass day^{-1}\)). This value was estimated from published data in Waters (1977; details_to_come).

The pb_ratio object can take three forms: fixed, gradient, or stochastic. A fixed pb_ratio will be set to the midpoint of the range argument, and conversions will not include any uncertainty. A gradient pb_ratio will take a set number of values in range, with the number of values determined by the length argument to build_pb_ratio(). The estimate_biomass() function will generate one set of nsim biomass values for each value in the pb_ratio gradient. A stochastic pb_ratio will sample values in range according to specified probabilities (the probs argument in build_pb_ratio()). The estimate_biomass() function will generate one set of nsim biomass values, where each of the nsim production values is assigned a random P:B value.

# create a P:B ratio object to link production to biomass
example_pb <- build_pb_ratio(range = c(0.25, 5.75), type = "fixed")

# use P:B ratio to estimate biomass from production
predicted_biomass <- estimate_biomass(production_estimates = predicted_production,
                                      pb_ratio = example_pb)

There are several helper functions to summarise the information generated by estimate_production or estimate_biomass. For example, it is possible to plot() the values for all nodes or a subset of nodes, where subsets can be based on row indexing or on node names. It is possible to extract values for any or all nodes, summarised by a function of your choice, using the extract_nodes() function.

# plot biomass for all nodes
plot(predicted_biomass)

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c(16:19))

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c("small.fish", "large.bodied.native.fish"))

# extract mean of biomass for a single node of interest
fish_biomass_mean <- extract_nodes(predicted_biomass, nodes = "large.bodied.native.fish",
                                   FUN = mean)

# extract quantiles of biomass for a single node of interest
fish_biomass_quantiles <- extract_nodes(predicted_biomass, nodes = "large.bodied.native.fish",
                                        FUN = quantile, p = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975))

example: predicting biomass of a high-order predator with uncertain food web

We can repeat the previous example but account for uncertainty in the stucture of the food web. In this case, the interaction matrix takes values in [0, 1] rather than being equal exactly to zero or one. For example, the example food web shown above can be altered so that links have a probability rather than definitely being present or absent.

Example stochastic food web
decapod.crustaceans small.fish carp reptiles mammals birds large.fish
decapod.crustaceans 0.00 0.00 0 0 0.00 0.00 0
small.fish 0.00 0.00 0 0 0.00 0.00 0
carp 0.00 0.00 0 0 0.00 0.00 0
reptiles 0.00 0.38 0 0 0.00 0.00 0
mammals 0.94 0.00 0 0 0.00 0.00 0
birds 0.80 0.96 0 0 0.00 0.00 0
large.bodied.native.fish 0.92 0.80 0 0 0.41 0.86 0

The production estimates are based on nsim replicates, with the food web drawn at random in each replicate based on the values (probabilities) in the interaction matrix.

# load trophic package
library(trophic)

# prepare food web, transfer efficiencies, and dominance matrix
stoch_food_web <- food_web * runif(length(food_web), min = 0, max = 1)
example_fw <- build_food_web(interaction_matrix = stoch_food_web)
example_eps <- build_efficiency_matrix(efficiency_mean = efficiency_mean,
                                       efficiency_sd = 0.05)
example_doms <- build_dominance_matrix(dominance = dominance_matrix)

# build primary producers object
example_pp <- build_primary_producers(production_mean = production_mean,
                                      production_sd = production_sd)

# collate trophic_dynamics object
example_trophic <- build_trophic_dynamics(food_web = example_fw,
                                          efficiency_matrix = example_eps,
                                          dominance_matrix = example_doms)

# predict production of all nodes
predicted_production <- estimate_production(trophic_dynamics = example_trophic,
                                            primary_producers = example_pp,
                                            stochastic = c("efficiency",
                                                           "primary_production"))

# create a P:B ratio object to link production to biomass
example_pb <- build_pb_ratio(range = c(0.25, 5.75), type = "fixed")

# use P:B ratio to estimate biomass from production
predicted_biomass <- estimate_biomass(production_estimates = predicted_production,
                                      pb_ratio = example_pb)

# plot biomass for all nodes
plot(predicted_biomass)

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c(16:19))

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c("small.fish", "large.bodied.native.fish"))

# extract mean of biomass for a single node of interest
fish_biomass_mean <- extract_nodes(predicted_biomass, nodes = "large.bodied.native.fish",
                                   FUN = mean)

# extract quantiles of biomass for a single node of interest
fish_biomass_quantiles <- extract_nodes(predicted_biomass, nodes = "large.bodied.native.fish",
                                        FUN = quantile, p = c(0.025, 0.1, 0.25, 0.5, 0.75, 0.9, 0.975))

example: calculating sensitivity of predicted biomass to transfer efficiencies

# load trophic package
library(trophic)

# prepare food web, transfer efficiencies, and dominance matrix
example_fw <- build_food_web(interaction_matrix = food_web)
example_eps <- build_efficiency_matrix(efficiency_mean = efficiency_mean,
                                       efficiency_sd = 0.05)
example_doms <- build_dominance_matrix(dominance = dominance_matrix)

# build primary producers object
example_pp <- build_primary_producers(production_mean = production_mean,
                                      production_sd = production_sd)

# collate trophic_dynamics object
example_trophic <- build_trophic_dynamics(food_web = example_fw,
                                          efficiency_matrix = example_eps,
                                          dominance_matrix = example_doms)

# predict production of all nodes
predicted_production <- estimate_production(trophic_dynamics = example_trophic,
                                            primary_producers = example_pp,
                                            stochastic = c("efficiency",
                                                           "primary_production"))

# create a P:B ratio object to link production to biomass
example_pb <- build_pb_ratio(range = c(0.25, 5.75), length = 4)

# use P:B ratio to estimate biomass from production
predicted_biomass <- estimate_biomass(production_estimates = predicted_production,
                                      pb_ratio = example_pb)

# plot biomass for all nodes
par(mfrow = c(2, 2))
plot(predicted_biomass)

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c(16:19))

# plot biomass for nodes of interest
plot(predicted_biomass, nodes = c("small.fish", "large.bodied.native.fish"))

# reset plot settings
par(mfrow = c(1, 1))

example: calculating sensitivity of predicted biomass to species’ dominances

# do something