,
Matthieu Bruneaux [aut, cre] | Title: | Isotopic Tracer Analysis Using MCMC |
|---|---|
| Description: | Implements Bayesian models to analyze data from tracer addition experiments. The implemented method was originally described in the article "A New Method to Reconstruct Quantitative Food Webs and Nutrient Flows from Isotope Tracer Addition Experiments" by López-Sepulcre et al. (2020) <doi:10.1086/708546>. |
| Authors: | Andrés López-Sepulcre [aut] ,
Matthieu Bruneaux [aut, cre] |
| Maintainer: | Matthieu Bruneaux <[email protected]> |
| License: | GPL-3 |
| Version: | 1.1.7 |
| Built: | 2024-11-05 06:22:10 UTC |
| Source: | https://gitlab.com/matthieu-bruneaux/isotracer |
The isotracer package allows modelling of fluxes across a network of compartments. Parameters are estimated using a Bayesian MCMC approach.
Maintainer: Matthieu Bruneaux [email protected] (ORCID)
Authors:
Andrés López-Sepulcre [email protected] (ORCID)
López-Sepulcre, A., M. Bruneaux, S. M. Collins, R. El-Sabaawi, A. S. Flecker, and S. A. Thomas. The American Naturalist (2020). "A New Method to Reconstruct Quantitative Food Webs and Nutrient Flows from Isotope Tracer Addition Experiments." https://doi.org/10.1086/708546.
Stan Development Team (2018). RStan: the R interface to Stan. R package version 2.18.2. https://mc-stan.org
Useful links:
Report bugs at https://gitlab.com/matthieu-bruneaux/isotracer/-/issues
networkModelStanfit objectsSubset method for networkModelStanfit objects
## S3 method for class 'networkModelStanfit' x[i, j, drop = TRUE]## S3 method for class 'networkModelStanfit' x[i, j, drop = TRUE]
x |
A |
i |
A vector of iteration indices. |
j |
A vector of parameter names or indices. |
drop |
Boolean. |
A networkModelStanfit object.
Note that new global parameters are not given any default prior.
add_covariates(nm, ..., use_regexpr = TRUE)add_covariates(nm, ..., use_regexpr = TRUE)
nm |
A |
... |
One or several formulas defining the covariates. |
use_regexpr |
Boolean, use regular expression to match the parameters affected by the formulas? |
A networkModel object.
# Using a subset of the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph") # Taking initial condtions from the 'lalaja' dataset at t=0 # Grouping by transect id inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m # Default model params(m, simplify = TRUE) # Adding an effect of the "transect" covariate on some parameters m <- add_covariates(m, upsilon_epi_to_pseph ~ transect) params(m, simplify = TRUE)# Using a subset of the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph") # Taking initial condtions from the 'lalaja' dataset at t=0 # Grouping by transect id inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m # Default model params(m, simplify = TRUE) # Adding an effect of the "transect" covariate on some parameters m <- add_covariates(m, upsilon_epi_to_pseph ~ transect) params(m, simplify = TRUE)
When applied to a steady-state compartment, this is equivalent to changing the steady state. Negative values are allowed, so one can add a "pulse" to a steady-state compartment and then later add a similar but negative "pulse" to simulate a drip in a stream for example.
add_pulse_event(nm, time, comp = NULL, unmarked, marked, which = NULL, pulses)add_pulse_event(nm, time, comp = NULL, unmarked, marked, which = NULL, pulses)
nm |
A |
time |
Numeric, time at which the pulse is happening. |
comp |
One compartment name only. |
unmarked |
Numeric, quantity of unmarked marker added. |
marked |
Numeric, quantity of marked marker added. |
which |
Vector of integers giving the nm rows to update. Default is to update all rows. |
pulses |
Optionally, a tibble containing the pulse information in columns. If provided, 'comp', 'time', 'unmarked' and 'marked' must be strings giving the corresponding column names. |
A networkModel object.
m <- trini_mod m$events <- NULL pulses <- tibble::tribble( ~ stream, ~ transect, ~ comp, ~ time, ~ qty_14N, ~ qty_15N, "UL", "transect.1", "NH4", 11, 0, -0.00569, "UL", "transect.2", "NH4", 11, 0, -0.00264, "UL", "transect.3", "NH4", 11, 0, -0.000726, "UL", "transect.1", "NO3", 11, 0, -0.00851, "UL", "transect.2", "NO3", 11, 0, -0.01118, "UL", "transect.3", "NO3", 11, 0, -0.01244, ) m <- add_pulse_event(m, pulses = pulses, comp = "comp", time = "time", unmarked = "qty_14N", marked = "qty_15N") mm <- trini_mod m$events <- NULL pulses <- tibble::tribble( ~ stream, ~ transect, ~ comp, ~ time, ~ qty_14N, ~ qty_15N, "UL", "transect.1", "NH4", 11, 0, -0.00569, "UL", "transect.2", "NH4", 11, 0, -0.00264, "UL", "transect.3", "NH4", 11, 0, -0.000726, "UL", "transect.1", "NO3", 11, 0, -0.00851, "UL", "transect.2", "NO3", 11, 0, -0.01118, "UL", "transect.3", "NO3", 11, 0, -0.01244, ) m <- add_pulse_event(m, pulses = pulses, comp = "comp", time = "time", unmarked = "qty_14N", marked = "qty_15N") m
This network model is the model used in the Quick Start tutorial
vignette. It is ready to be run at once with run_mcmc.
aquarium_modaquarium_mod
An object of class networkModel (inherits from tbl_df, tbl, data.frame) with 1 rows and 4 columns.
The code used to built the model is given in the example section below.
The aquarium_run dataset is a corresponding MCMC run.
library(tibble) library(dplyr) exp <- tibble::tribble( ~time.day, ~species, ~biomass, ~prop15N, 0, "algae", 1.02, 0.00384, 1, "algae", NA, 0.0534, 1.5, "algae", 0.951, NA, 2, "algae", 0.889, 0.0849, 2.5, "algae", NA, 0.0869, 3, "algae", 0.837, 0.0816, 0, "daphnia", 1.74, 0.00464, 1, "daphnia", NA, 0.00493, 1.5, "daphnia", 2.48, NA, 2, "daphnia", NA, 0.00831, 2.5, "daphnia", 2.25, NA, 3, "daphnia", 2.15, 0.0101, 0, "NH4", 0.208, 0.79, 1, "NH4", 0.227, NA, 1.5, "NH4", NA, 0.482, 2, "NH4", 0.256, 0.351, 2.5, "NH4", NA, 0.295, 3, "NH4", 0.27, NA ) inits <- exp %>% dplyr::filter(time.day == 0) obs <- exp %>% dplyr::filter(time.day > 0) aquarium_mod <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia -> NH4") %>% set_init(inits, comp = "species", size = "biomass", prop = "prop15N") %>% set_obs(obs, comp = "species", size = "biomass", prop = "prop15N", time = "time.day")library(tibble) library(dplyr) exp <- tibble::tribble( ~time.day, ~species, ~biomass, ~prop15N, 0, "algae", 1.02, 0.00384, 1, "algae", NA, 0.0534, 1.5, "algae", 0.951, NA, 2, "algae", 0.889, 0.0849, 2.5, "algae", NA, 0.0869, 3, "algae", 0.837, 0.0816, 0, "daphnia", 1.74, 0.00464, 1, "daphnia", NA, 0.00493, 1.5, "daphnia", 2.48, NA, 2, "daphnia", NA, 0.00831, 2.5, "daphnia", 2.25, NA, 3, "daphnia", 2.15, 0.0101, 0, "NH4", 0.208, 0.79, 1, "NH4", 0.227, NA, 1.5, "NH4", NA, 0.482, 2, "NH4", 0.256, 0.351, 2.5, "NH4", NA, 0.295, 3, "NH4", 0.27, NA ) inits <- exp %>% dplyr::filter(time.day == 0) obs <- exp %>% dplyr::filter(time.day > 0) aquarium_mod <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia -> NH4") %>% set_init(inits, comp = "species", size = "biomass", prop = "prop15N") %>% set_obs(obs, comp = "species", size = "biomass", prop = "prop15N", time = "time.day")
This is an MCMC run on aquarium_mod. The code used to run the
MCMC is: aquarium_run <- run_mcmc(aquarium_mod, thin = 4) (note that
thin = 4 was only used here to reduce the size of the data file
shipped with the package, but for a real-life analysis keeping the default
thin = 1 is usually recommended). The code used to build the model
itself is given in the help page for aquarium_mod.
aquarium_runaquarium_run
An object of class networkModelStanfit (inherits from mcmc.list) of length 4.
## Not run: plot(aquarium_run) summary(aquarium_run) ## End(Not run)## Not run: plot(aquarium_run) summary(aquarium_run) ## End(Not run)
Convert a compatible object to a tbl_graph object (from the tidygraph package)
as_tbl_graph(x, ...)as_tbl_graph(x, ...)
x |
Object to convert to a tbl_graph. |
... |
Passed to the appropriate method. |
A tbl_graph object.
Convert a network topology to a tbl_graph
## S3 method for class 'topology' as_tbl_graph(x, ...)## S3 method for class 'topology' as_tbl_graph(x, ...)
x |
A network topology. |
... |
Not used. |
A tbl_graph object.
tidy_flows object to an mcmc.list
Convert a tidy_flows object to an mcmc.list
## S3 method for class 'tidy_flows' as.mcmc.list(x, ...)## S3 method for class 'tidy_flows' as.mcmc.list(x, ...)
x |
A tidy flow object, as returned by |
... |
Not used for now. |
A mcmc.list object, with ordered iterations.
tidy_steady_states object to an mcmc.list
Convert a tidy_steady_states object to an mcmc.list
## S3 method for class 'tidy_steady_states' as.mcmc.list(x, ...)## S3 method for class 'tidy_steady_states' as.mcmc.list(x, ...)
x |
A tidy steady states object, as returned by
|
... |
Not used for now. |
A mcmc.list object, with ordered iterations.
List the available priors for model parameters
available_priors()available_priors()
A tibble containing information about the available priors.
available_priors()available_priors()
Combine mcmc.list objects
## S3 method for class 'mcmc.list' c(...)## S3 method for class 'mcmc.list' c(...)
... |
|
A mcmc.list object.
This is an experimental function. It attempts to calculate steady-state compartment sizes using the set parameter values and the initial compartment sizes. Use it with caution!
calculate_steady_state(nm)calculate_steady_state(nm)
nm |
A network model, with set parameter values. |
Note about how steady state sizes for split compartments are calculated: the steady size of the active portion is calculated divide it is divided by the active fraction (portion.act parameter) to get the total size including the refractory portion. In this case we get a "steady-state" refractory portion, consistent with steady state size of active fraction and with portion.act parameter.
A tibble containing steady-state compartment sizes.
m <- aquarium_mod m <- set_prior(m, constant_p(0), "lambda") m <- set_params(m, sample_params(m)) proj <- project(m, end = 40) plot(proj) z <- calculate_steady_state(m) z z$stable_sizesm <- aquarium_mod m <- set_prior(m, constant_p(0), "lambda") m <- set_params(m, sample_params(m)) proj <- project(m, end = 40) plot(proj) z <- calculate_steady_state(m) z z$stable_sizes
Return the compartments of a network model
comps(nm)comps(nm)
nm |
A |
A list of character vectors, with one list element per row of the input network model (list elements are in the same order as the input network model rows). Each list element containing the names of the compartments in the topology defined in the corresponding row of the input network model.
aquarium_mod comps(aquarium_mod) trini_mod comps(trini_mod)aquarium_mod comps(aquarium_mod) trini_mod comps(trini_mod)
This is equivalent to having a fixed parameter.
constant_p(value)constant_p(value)
value |
The constant value of the parameter. |
A list defining the prior.
constant_p(2)constant_p(2)
For details and references about quantities used in expressing isotopic ratios, see:
delta2prop(x = NULL, Rstandard = NULL)delta2prop(x = NULL, Rstandard = NULL)
x |
Vector of delta values. |
Rstandard |
String describing the isotopic measurement, e.g. "d15N", "d13C" and used to set automatically Rstandards (see the Section "Ratios for reference standards" for more details). Alternatively, a numeric value to use for Rstandard, e.g. 0.0036765. |
- Figure 1 in Coplen, Tyler B. “Guidelines and Recommended Terms for Expression of Stable-Isotope-Ratio and Gas-Ratio Measurement Results.” Rapid Communications in Mass Spectrometry 25, no. 17 (September 15, 2011): 2538–60. https://doi.org/10.1002/rcm.5129.
- Table 2.1 in Fry, Brian. Stable Isotope Ecology. New York: Springer-Verlag, 2006. //www.springer.com/gp/book/9780387305134.
A vector of same length of x, containing the proportion (numeric between 0 and 1) of heavy isotope based on the delta values and the Rstandard provided.
The ratios for reference standards are taken from the Table 2.1 from Fry 2006. Note that the values used for oxygen isotopes are from the standard mean ocean water (SMOW).
Standards recognized by this function are: c("d15N", "d2H", "d13C",
"d17O.SMOW", "d18O.SMOW", "d33S", "d34S", "d36S")
deltas <- c(78, 5180, 263, 1065, NA, 153, 345) # Rstandard can be specified with a string for some preset references prop15N <- delta2prop(deltas, "d15N") prop13C <- delta2prop(deltas, "d13C") # Rstandard can also be specified manually for non-preset references prop15N_manual <- delta2prop(deltas, 0.0036765) prop13C_manual <- delta2prop(deltas, 0.011180) # Call delta2prop() to get the detail of available references delta2prop()deltas <- c(78, 5180, 263, 1065, NA, 153, 345) # Rstandard can be specified with a string for some preset references prop15N <- delta2prop(deltas, "d15N") prop13C <- delta2prop(deltas, "d13C") # Rstandard can also be specified manually for non-preset references prop15N_manual <- delta2prop(deltas, 0.0036765) prop13C_manual <- delta2prop(deltas, 0.011180) # Call delta2prop() to get the detail of available references delta2prop()
Note that DIC might not be indicated for network models, as the posteriors are often not multinormal distributions.
dic(..., weight = TRUE)dic(..., weight = TRUE)
... |
One or several |
weight |
Boolean, if TRUE calculate DIC weights based on Link and Barker 2010 (Link, W. A., and R. J. Barker. 2010. Bayesian Inference With Ecological Applications. Amsterdam Boston Heidelberg London: Elsevier/Academic Press). |
LOO is probably not a good choice either since the data is akin to a time series (so data points are not independent). Maybe WAIC could be an option? (TODO: read about this.)
DIC is calculated as:
DIC = Dbar + pD
where D are deviance values calculated as -2 * loglik for each MCMC iteration, Dbar is the mean deviance value and pD is the effective number of parameters in the model and can be calculated as var(D)/2 (Gelman 2003).
A tibble with one row per mcmc.list object provided in
.... This tibble is sorted by DIC, so the row order might be
different from the mcmc.list objects order.
# Define two different models m1 <- aquarium_mod m2 <- set_topo(m1, c("NH4 -> algae -> daphnia -> NH4", "algae -> NH4")) m2 <- set_priors(m2, priors(m1)) m2 <- set_priors(m2, normal_p(0, 0.5), "upsilon_algae_to_NH4") # Run the models r1 <- run_mcmc(m1, chains = 2) r2 <- run_mcmc(m2, chains = 2) # Model comparison with DIC dic(r1, r2)# Define two different models m1 <- aquarium_mod m2 <- set_topo(m1, c("NH4 -> algae -> daphnia -> NH4", "algae -> NH4")) m2 <- set_priors(m2, priors(m1)) m2 <- set_priors(m2, normal_p(0, 0.5), "upsilon_algae_to_NH4") # Run the models r1 <- run_mcmc(m1, chains = 2) r2 <- run_mcmc(m2, chains = 2) # Model comparison with DIC dic(r1, r2)
Dataset built from the article "Phosphate absorption in eelgrass" by McRoy and Barsdate (1970)
eelgrasseelgrass
Tibble with columns
Light treatment: "light" or "dark".
The location where 32P phosphate was added: in the "upper" water compartment or in the "lower" water compartment.
Obsered compartment, one of "leaves_stem", "roots_rhizome", "upper_water", or "lower_water".
Elapsed time in minutes since the 32P addition.
Number of 32P atoms per mg (estimated from Figure 2 of the original paper).
Compartment mass in mg (taken from Table 1 of the original paper). Assumed constant during the experiment.
Number of 32P atoms in the compartment. Calculated from the two previous columns.
In brief, the experimental setup consists in individual eelgrass plants placed in 250 ml containers. Each container is partitioned by a layer of paraffin into an upper water compartment (containing the leaves and stems) and a lower water compartment (containing the roots and rhizomes).
Radioactive phosphorus (32P) is added as phosphate either in the upper or lower water compartment in each container. Containers were incubated either in light or dark conditions.
Tissue samples were collected and dried at various time points and 32P activity was measured (Figure 2 in the original paper). Biomass estimates in initial conditions were given in Table 1 of the original paper.
The data for 32P abundance per mg is extracted from Figure 2 of the original article. Atom counts per mg were derived from cpm per mg using a half-life value of 14.268 days for 32P.
For simplicity and in order to be able to match the 32P data with the biomass data (see below), only four compartments are considerd in the package dataset. Upper and lower water compartments match the compartments from the original article. "Leaf and stem" pools the original compartments "leaf tip", "leaf middle", "leaf base", and "stem". "Roots and rhizome" pools the original compartments "root" and "rhizome". Pooling is done by averaging the cpm per mg data, thereby making the rough approximation that each component of the pool contributes the same biomass as the other components.
The biomass data is taken from Table 1 in the original paper. Experimental containers had 160 cc of seawater in the upper compartment and 80 cc of seawater in the lower compartment. Based on comparison with data from Risgaard-Petersen 1998, I assumed that the biomasses for tissues were given in dry weight. I assumed that this was also the case for the cpm/mg data (i.e. cpm/mg of dry weight).
Data was taken from the figures and tables of the original paper. The original paper is: McRoy, C. Peter, and Robert J. Barsdate. “Phosphate Absorption in Eelgrass1.” Limnology and Oceanography 15, no. 1 (January 1, 1970): 6–13. https://doi.org/10.4319/lo.1970.15.1.0006.
Define an exponential prior
exponential_p(lambda)exponential_p(lambda)
lambda |
Lambda parameter (rate) of the exponential distribution. The mean of the exponential distribution is 1/lambda. |
A list defining the prior.
exponential_p(0.5)exponential_p(0.5)
Filter (alias for filter function from dplyr)
.data |
Data to filter. |
... |
Passed to dplyr::filter. |
preserve |
Ignored. |
See the returned value for dplyr::filter.
This function can be used to filter any tibble (e.g. network model object) that has a "group" column. See the Examples for more details and syntax.
filter_by_group(.data, ...)filter_by_group(.data, ...)
.data |
A tibble that has a 'group' column, such as a 'networkModel' object. |
... |
Conditional expressions for filtering (see the Examples). |
A tibble similar to the input object, but with rows filtered based
on ....
trini_mod trini_mod$group groups(trini_mod) filter_by_group(trini_mod, stream == "LL", transect == "transect.1") filter_by_group(trini_mod, transect == "transect.1") ## Not run: # The code below would raise an error because there is no "color" grouping variable. filter_by_group(trini_mod, color == "red") ## End(Not run)trini_mod trini_mod$group groups(trini_mod) filter_by_group(trini_mod, stream == "LL", transect == "transect.1") filter_by_group(trini_mod, transect == "transect.1") ## Not run: # The code below would raise an error because there is no "color" grouping variable. filter_by_group(trini_mod, color == "red") ## End(Not run)
Filter method for output of tidy_data_and_posterior_predict()
## S3 method for class 'ppcNetworkModel' filter(.data, ..., .preserve = FALSE)## S3 method for class 'ppcNetworkModel' filter(.data, ..., .preserve = FALSE)
.data |
A ppcNetworkModel object. |
... |
Passed to dplyr::filter. |
.preserve |
Ignored. |
A pccNetworkModel object filtered appropriately based on the [["vars"]] tibble.
prior objectPretty formatting of a prior object
## S3 method for class 'prior' format(x, ...)## S3 method for class 'prior' format(x, ...)
x |
An object of class |
... |
Not used. |
A character string for pretty printing of a prior.
prior_tibble objectPretty formatting of a prior_tibble object
## S3 method for class 'prior_tibble' format(x, ...)## S3 method for class 'prior_tibble' format(x, ...)
x |
An object of class |
... |
Not used. |
A character string for pretty printing of a prior tibble.
Note the name of the function to define a prior (gamma_p), in order
to avoid confusion with the R mathematical function gamma.
gamma_p(alpha, beta)gamma_p(alpha, beta)
alpha |
Shape parameter (equivalent to the |
beta |
Rate parameter (equivalent to the |
A list defining the prior.
gamma_p(9, 2) hist(sample_from_prior(gamma_p(9, 2), 1e3))gamma_p(9, 2) hist(sample_from_prior(gamma_p(9, 2), 1e3))
A quick-and-dirty way of visualizing relative flows in a network
ggflows(x, layout = "auto", edge = "fan", max_width, legend = TRUE, ...)ggflows(x, layout = "auto", edge = "fan", max_width, legend = TRUE, ...)
x |
A tibble with the flow estimates, with columns "from", "to", and "flow". |
layout |
Optional, layout to use (e.g. "sugiyama", "kk", "stress") |
edge |
"curve" (the default), "line" or "fan". |
max_width |
Optional, numeric giving the maximum edge width (minimum width is always 1). |
legend |
Boolean, display edge width legend? |
... |
Not used. |
A ggplot2 plot.
if (requireNamespace("ggraph")) { z <- tibble::tribble( ~from, ~to, ~flow, "leavesAndStem", "rootsAndRhizome", 333.929866077124, "lowerWater", "rootsAndRhizome", 4425.15780019304, "rootsAndRhizome", "leavesAndStem", 525.208837577916, "upperWater", "leavesAndStem", 11224.0814971855 ) ggflows(z) ggflows(z, max_width = 15) }if (requireNamespace("ggraph")) { z <- tibble::tribble( ~from, ~to, ~flow, "leavesAndStem", "rootsAndRhizome", 333.929866077124, "lowerWater", "rootsAndRhizome", 4425.15780019304, "rootsAndRhizome", "leavesAndStem", 525.208837577916, "upperWater", "leavesAndStem", 11224.0814971855 ) ggflows(z) ggflows(z, max_width = 15) }
A quick plot using ggraph
ggtopo(x, layout = "auto", edge = "fan", ...)ggtopo(x, layout = "auto", edge = "fan", ...)
x |
A network model or a topology matrix. |
layout |
Optional, layout to use (e.g. "sugiyama", "kk", "stress") |
edge |
"fan" (the default) or "line" or "curve". |
... |
Passed to the methods. |
A ggplot2 plot.
if (requireNamespace("ggraph")) { ggtopo(aquarium_mod, edge = "line") }if (requireNamespace("ggraph")) { ggtopo(aquarium_mod, edge = "line") }
A quick plot using ggraph
## S3 method for class 'networkModel' ggtopo(x, layout = "auto", edge = "fan", ...)## S3 method for class 'networkModel' ggtopo(x, layout = "auto", edge = "fan", ...)
x |
A topology matrix. |
layout |
Optional, layout to use (e.g. "sugiyama", "kk", "stress") |
edge |
"curve" (the default) or "line". |
... |
Not used for now. |
A ggplot2 plot.
if (requireNamespace("ggraph")) { ggtopo(aquarium_mod, edge = "line") ggtopo(trini_mod) }if (requireNamespace("ggraph")) { ggtopo(aquarium_mod, edge = "line") ggtopo(trini_mod) }
A quick plot using ggraph
## S3 method for class 'topology' ggtopo(x, layout = "auto", edge = "fan", ...)## S3 method for class 'topology' ggtopo(x, layout = "auto", edge = "fan", ...)
x |
A topology matrix. |
layout |
Optional, layout to use (e.g. "sugiyama", "kk", "stress") |
edge |
"curve" (the default), "line" or "fan". |
... |
Not used for now. |
A ggplot2 plot.
if (requireNamespace("ggraph")) { z <- topo(aquarium_mod) ggtopo(z) ggtopo(z, edge = "line") z <- topo(trini_mod) ggtopo(z) # For finer control, one can build a tbl_graph from the topology and # use ggraph directly x <- as_tbl_graph(z) library(ggraph) ggraph(x) + geom_edge_link() }if (requireNamespace("ggraph")) { z <- topo(aquarium_mod) ggtopo(z) ggtopo(z, edge = "line") z <- topo(trini_mod) ggtopo(z) # For finer control, one can build a tbl_graph from the topology and # use ggraph directly x <- as_tbl_graph(z) library(ggraph) ggraph(x) + geom_edge_link() }
networkModel objectGet the grouping for a networkModel object
## S3 method for class 'networkModel' groups(x)## S3 method for class 'networkModel' groups(x)
x |
A |
A tibble giving the grouping variable(s) for the input network
model. This tibble is in the same order as the rows of the input network
model. If the input network model did not have any grouping variable,
returns NULL.
groups(aquarium_mod) groups(trini_mod)groups(aquarium_mod) groups(trini_mod)
Define a half-Cauchy prior (on [0;+Inf])
hcauchy_p(scale)hcauchy_p(scale)
scale |
Median of the half-Cauchy distribution. |
A list defining the prior.
hcauchy_p(scale = 0.5)hcauchy_p(scale = 0.5)
Dataset built from the article "Fish introductions and light modulate food web fluxes in tropical streams: a whole-ecosystem experimental approach” by Collins et al. (2016).
lalajalalaja
Tibble with columns
Stream identity. It is always "UL" (for "Upper lalaja") in
this dataset. See the model trini_mod also shipped with the
package for the full dataset from the original Collins et al. study,
including data from the Lower Lajaja stream.
Transect identity. Three transects were sampled downstream
of the drip location: c("transect.1", "transect.2",
"transect.3").
Foodweb compartments. Eight compartments are included in this dataset: "NH4", dissolved ammonium; "NH3", dissolved nitrate; "epi", epilithon (primary producers growing on the surface of rocks on the stream bed); "FBOM", fine benthic organic material; "tricor", Tricorythodes (invertebrate); "pseph", Psephenus (invertebrate); "petro", Petrophila (invertebrate); "arg", Argia (invertebrate).
Size of compartment, in mg of nitrogen per m2.
Proportion of 15N nitrogen in a compartment nitrogen pool (i.e. 15N / (15N + 14N)).
Sampling time, in days.
In the original study, 15N-enriched ammonium was dripped into two mountain
streams in Trinidad (Upper Lalaja stream and Lower Lalaja stream) and
samples of the different foodweb compartments were taken during the drip and
after the drip in several transects in each stream. The transects were
located at different locations downstream of each drip. There were three
transects per stream. The drip phase lasted 10 days, and the post-drip phase
lasted 30 days. The complete dataset from the original study is available in the
trini_mod model shipped with the isotracer package.
The lalaja dataset is a subset of the full dataset and is used for
illustrative purpose in the "Trinidadian streams" case study, which is part
of the documentation of isotracer. It contains only the data for the
Upper Lalaja stream, and for some but not all of the foodweb compartments.
For more details about the dripping regime and how to use this dataset in a
network model, one should refer to the case study in the isotracer
package documentation.
This network model contains data from the original article: Collins, Sarah M., Steven A. Thomas, Thomas Heatherly, Keeley L. MacNeill, Antoine O.H.C. Leduc, Andrés López-Sepulcre, Bradley A. Lamphere, et al. 2016. “Fish Introductions and Light Modulate Food Web Fluxes in Tropical Streams: A Whole-Ecosystem Experimental Approach.” Ecology, <doi:10.1002/ecy.1530>.
This dataset was also used in the paper: López-Sepulcre, Andrés, Matthieu Bruneaux, Sarah M. Collins, Rana El-Sabaawi, Alexander S. Flecker, and Steven A. Thomas. 2020. “A New Method to Reconstruct Quantitative Food Webs and Nutrient Flows from Isotope Tracer Addition Experiments.” The American Naturalist 195 (6): 964–85. <doi:10.1086/708546>.
Dataset built from the Dryad depository entry associated with the article "Protein degradation rate in Arabidopsis thaliana leaf growth and development" by Li et al. (2017)
li2017li2017
li2017 is the main dataset and is a tibble with columns:
Protein identifier. Can be matched to a more explicit
protein description in li2017_prots.
Sample identity. Different samples were used for relative abundance measurements and labelled fraction measurements.
Relative abundance compared to a reference sample.
Proportion of 15N in the protein.
Time elapsed since growth medium switch to 15N, in days.
Leaf identity (3rd, 5th, or 7th leaf of individual plants).
li2017_prots maps protein identifiers to protein descriptions
and is a tibble with columns:
Protein identifier. Can be matched with the same column in
li2017.
Protein description
li2017_counts is a summary table counting the number of
available data points for relative abundance and labelled fraction for
each protein in li2017. It is a tibble with columns:
Protein identifier. Can be matched with the same column in
li2017.
Number of relative abundance data points for a given protein.
Number of labelled fraction data points for a given protein.
In this study, the authors used a growth medium containing 15N to grow 21-day old Arabidopsis plants which were grown on a natural 14N/15N medium until that day. The third, fifth and seventh leaves were sampled from individuals at different time points after the medium switch (0, 1, 3 and 5 days). Proteins were identified and labelled fractions were measured using mass spectrometry. Relative protein abundances were determined in comparison with a reference sample.
The aim of the authors was to quantify in vivo degradation rates for as many proteins as possible (1228 proteins in the original paper) and examine which determinants had an effect or not on protein degradation rates (e.g. protein domains, protein complex membership, ...).
Three datasets were extracted from the large dataset available on Dryad for
packaging inside isotracer: li2017, li2017_prots, and
li2017_counts.
Data was taken from the following Dryad repository: Li, Lei, Clark J. Nelson, Josua Troesch, Ian Castleden, Shaobai Huang, and A. Harvey Millar. “Data from: Protein Degradation Rate in Arabidopsis Thaliana Leaf Growth and Development.” Dryad, 2018. https://doi.org/10.5061/DRYAD.Q3H85.
The Dryad repository was associated with the following paper: Li, Lei, Clark J. Nelson, Josua Trösch, Ian Castleden, Shaobai Huang, and A. Harvey Millar. “Protein Degradation Rate in Arabidopsis Thaliana Leaf Growth and Development.” The Plant Cell 29, no. 2 (February 1, 2017): 207–28. https://doi.org/10.1105/tpc.16.00768.
Math generics for mcmc.list objects
## S3 method for class 'mcmc.list' Math(x, ...)## S3 method for class 'mcmc.list' Math(x, ...)
x |
|
... |
Other arguments passed to corresponding methods |
A mcmc.list object (with the added class
derived.mcmc.list).
Note that the colors represent the strength of the correlations (from 0 to 1), but do not inform about their sign. The method used to calculate correlation coefficients is Spearman's rho.
mcmc_heatmap(x, col = NULL, ...)mcmc_heatmap(x, col = NULL, ...)
x |
A |
col |
Optional, vectors of colors defining the color ramp. Default uses the divergent palette "Blue-Red 2" from the colorspace package. |
... |
Passed to |
Called for side effect (plotting).
Get a table with parameters which are missing priors
missing_priors(nm)missing_priors(nm)
nm |
A |
A tibble containing the parameters which are missing a prior. If no priors are missing, the tibble contains zero row.
# Using a subset of the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph") # No prior is set by default priors(m) # Set some priors m <- set_priors(m, normal_p(0, 10), "lambda") priors(m) # Which parameters are missing a prior? missing_priors(m)# Using a subset of the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph") # No prior is set by default priors(m) # Set some priors m <- set_priors(m, normal_p(0, 10), "lambda") priors(m) # Which parameters are missing a prior? missing_priors(m)
The first step in building a network model is to create a new, empty
networkModel object. This model can then be completed using functions
such as set_topo(), set_init(), etc...
new_networkModel(quiet = FALSE)new_networkModel(quiet = FALSE)
quiet |
Boolean, if |
An empty networkModel object. It is basically a zero-row
tibble with the appropriate columns.
m <- new_networkModel() m class(m)m <- new_networkModel() m class(m)
Define a truncated normal prior (on [0;+Inf])
normal_p(mean, sd)normal_p(mean, sd)
mean |
Mean of the untruncated normal. |
sd |
Standard deviation of the untruncated normal. |
A list defining the prior.
normal_p(mean = 0, sd = 4)normal_p(mean = 0, sd = 4)
prior object in tibblesFunction used for displaying prior object in tibbles
## S3 method for class 'prior' obj_sum(x)## S3 method for class 'prior' obj_sum(x)
x |
An object of class |
Input formatted with format(x).
mcmc.list objectsOps generics for mcmc.list objects
## S3 method for class 'mcmc.list' Ops(e1, e2)## S3 method for class 'mcmc.list' Ops(e1, e2)
e1 |
First operand |
e2 |
Second operand |
A mcmc.list object (with the added class
derived.mcmc.list).
## Not run: # aquarium_run is a coda::mcmc.list object shipped with the isotracer package a <- aquarium_run plot(a) # The calculations below are just given as examples of mathematical # operations performed on an mcmc.list object, and do not make any sense # from a modelling point of view. plot(a[, "upsilon_algae_to_daphnia"] - a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] + a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] / a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] * a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] - 10) plot(a[, "upsilon_algae_to_daphnia"] + 10) plot(a[, "upsilon_algae_to_daphnia"] * 10) plot(a[, "upsilon_algae_to_daphnia"] / 10) plot(10 - a[, "upsilon_algae_to_daphnia"]) plot(10 + a[, "upsilon_algae_to_daphnia"]) plot(10 * a[, "upsilon_algae_to_daphnia"]) plot(10 / a[, "upsilon_algae_to_daphnia"]) ## End(Not run)## Not run: # aquarium_run is a coda::mcmc.list object shipped with the isotracer package a <- aquarium_run plot(a) # The calculations below are just given as examples of mathematical # operations performed on an mcmc.list object, and do not make any sense # from a modelling point of view. plot(a[, "upsilon_algae_to_daphnia"] - a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] + a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] / a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] * a[, "lambda_algae"]) plot(a[, "upsilon_algae_to_daphnia"] - 10) plot(a[, "upsilon_algae_to_daphnia"] + 10) plot(a[, "upsilon_algae_to_daphnia"] * 10) plot(a[, "upsilon_algae_to_daphnia"] / 10) plot(10 - a[, "upsilon_algae_to_daphnia"]) plot(10 + a[, "upsilon_algae_to_daphnia"]) plot(10 * a[, "upsilon_algae_to_daphnia"]) plot(10 / a[, "upsilon_algae_to_daphnia"]) ## End(Not run)
Implementation of the '==' operator for priors
## S3 method for class 'prior' Ops(e1, e2)## S3 method for class 'prior' Ops(e1, e2)
e1, e2
|
Objects of class "prior". |
Boolean (or throws an error for unsupported operators).
p <- constant_p(0) q <- constant_p(4) p == q p <- hcauchy_p(2) q <- hcauchy_p(2) p == qp <- constant_p(0) q <- constant_p(4) p == q p <- hcauchy_p(2) q <- hcauchy_p(2) p == q
topology objectsOps generics for topology objects
## S3 method for class 'topology' Ops(e1, e2)## S3 method for class 'topology' Ops(e1, e2)
e1 |
First operand |
e2 |
Second operand |
Boolean (or throws an error for unsupported operators).
topo(aquarium_mod) == topo(trini_mod) topo(aquarium_mod) == topo(aquarium_mod)topo(aquarium_mod) == topo(trini_mod) topo(aquarium_mod) == topo(aquarium_mod)
Return the parameters of a network model
params(nm, simplify = FALSE)params(nm, simplify = FALSE)
nm |
A |
simplify |
If |
A tibble containing the parameter names and their current value (if
set). If simplify is TRUE, only return a sorted character
vector containing the parameters names.
params(aquarium_mod) params(trini_mod) params(trini_mod, simplify = TRUE)params(aquarium_mod) params(trini_mod) params(trini_mod, simplify = TRUE)
prior object in tibblesFunction used for displaying prior object in tibbles
## S3 method for class 'prior' pillar_shaft(x, ...)## S3 method for class 'prior' pillar_shaft(x, ...)
x |
An object of class |
... |
Not used. |
An object prepared with pillar::new_pillar_shaft_simple.
Plot observations/trajectories/predictions from a network model
## S3 method for class 'networkModel' plot(x, ...)## S3 method for class 'networkModel' plot(x, ...)
x |
A |
... |
Passed to |
Called for side effect (plotting).
split_to_unit_plot
Plot output from split_to_unit_plot
## S3 method for class 'ready_for_unit_plot' plot(x, ...)## S3 method for class 'ready_for_unit_plot' plot(x, ...)
x |
A |
... |
Passed to |
Called for side effect (plotting).
Draw from the posterior predictive distribution of the model outcome
posterior_predict(object, ...)posterior_predict(object, ...)
object |
Model from which posterior predictions can be made. |
... |
Passed to the appropriate method. |
Usually methods will implement a draw parameter, and the
returned object is a "draw" by N matrix where N is the number of data
points predicted per draw.
Draw from the posterior predictive distribution of the model outcome
## S3 method for class 'networkModelStanfit' posterior_predict(object, newdata, draw = NULL, cores = NULL, ...)## S3 method for class 'networkModelStanfit' posterior_predict(object, newdata, draw = NULL, cores = NULL, ...)
object |
A networkModelStanfit object. |
newdata |
Should be the model used to fit the networkStanfit object. |
draw |
Integer, number of draws to perform from the posterior. Default is 100. |
cores |
Number of cores to use for parallel calculations. Default is
|
... |
Not used for now. |
A "draw" by N matrix where N is the number of data points predicted per draw.
Add a column with predictions from a fit
## S3 method for class 'networkModel' predict( object, fit, draws = NULL, error.draws = 5, probs = 0.95, cores = NULL, dt = NULL, grid_size = NULL, at = NULL, end = NULL, ... )## S3 method for class 'networkModel' predict( object, fit, draws = NULL, error.draws = 5, probs = 0.95, cores = NULL, dt = NULL, grid_size = NULL, at = NULL, end = NULL, ... )
object |
Network model |
fit |
Model fit (mcmc.list object) |
draws |
Integer, number of draws from the posteriors |
error.draws |
Integer, number of draws from the error distribution, for a given posterior draw. |
probs |
Credible interval (default 0.95). |
cores |
Number of cores to use for parallel calculations. Default is
|
dt, grid_size
|
Time step size or grid points, respectively. |
at |
Timepoints at which the predictions should be returned. |
end |
Final timepoint used in the projections. |
... |
Not used. |
A network model object with an added column "prediction".
networkModel objectsPrint method for networkModel objects
## S3 method for class 'networkModel' print(x, ...)## S3 method for class 'networkModel' print(x, ...)
x |
A |
... |
Passsed to the next method. |
Called for the side effect of printing a network model object.
prior objectPretty printing of a prior object
## S3 method for class 'prior' print(x, ...)## S3 method for class 'prior' print(x, ...)
x |
An object of class |
... |
Not used. |
Mostly called for its side effect of printing, but also returns its input invisibly.
prior_tibble objectPretty printing of a prior_tibble object
## S3 method for class 'prior_tibble' print(x, ...)## S3 method for class 'prior_tibble' print(x, ...)
x |
An object of class |
... |
Not used. |
Mostly called for its side effect of printing, but also returns its input invisibly.
topology objectPretty printing of a topology object
## S3 method for class 'topology' print(x, help = TRUE, ...)## S3 method for class 'topology' print(x, help = TRUE, ...)
x |
An object of class |
help |
If TRUE, display a short help after the topology object explaining e.g. the steady state or the split compartment symbols. |
... |
Not used. |
Mostly called for its side effect (printing).
Return the tibble containing the priors of a networkModel
priors(nm, fix_set_params = FALSE, quiet = FALSE)priors(nm, fix_set_params = FALSE, quiet = FALSE)
nm |
A |
fix_set_params |
If TRUE, parameters for which a value is set are given a fixed value (i.e. their prior is equivalent to a point value). |
quiet |
Boolean to control verbosity. |
A tibble giving the current priors defined for the input network model.
priors(aquarium_mod) priors(trini_mod)priors(aquarium_mod) priors(trini_mod)
Calculate the trajectories of a network model
project( nm, dt = NULL, grid_size = NULL, at = NULL, end = NULL, flows = "no", cached_ts = NULL, cached_ee = NULL, ignore_pulses = FALSE )project( nm, dt = NULL, grid_size = NULL, at = NULL, end = NULL, flows = "no", cached_ts = NULL, cached_ee = NULL, ignore_pulses = FALSE )
nm |
A |
dt, grid_size
|
Either the time step size for trajectory calculations
( |
at |
Optional, vector of time values at which the trajectory must be evaluated. |
end |
Time value for end point. If not provided, the last observation or event is used. |
flows |
Return flow values? The default is "no" and no flows are calculated. Other values are "total" (total flows summed up from beginning to end timepoint), "average" (average flows per time unit, equal to total flows divided by the projection duration), and "per_dt" (detailled flow values are returned for each interval dt of the projection). |
cached_ts, cached_ee
|
Used for optimization by other functions, not for use by the package user. |
ignore_pulses |
Default to FALSE (i.e. apply pulses when projecting the network system). It is set to TRUE when calculating steady-state flows. |
A network model object with a "trajectory" column.
m <- aquarium_mod m <- set_params(m, sample_params(m)) z <- project(m) z <- project(m, flows = "per_dt") z <- project(m, flows = "total") z <- project(m, flows = "average")m <- aquarium_mod m <- set_params(m, sample_params(m)) z <- project(m) z <- project(m, flows = "per_dt") z <- project(m, flows = "total") z <- project(m, flows = "average")
Return the distribution family for observed proportions
prop_family(nm, quiet = FALSE)prop_family(nm, quiet = FALSE)
nm |
A |
quiet |
Boolean for being quiet about explaining the role of eta
(default is |
A character string describing the distribution family used to model observed proportions.
prop_family(aquarium_mod) prop_family(trini_mod)prop_family(aquarium_mod) prop_family(trini_mod)
This function performs the inverse of the operation performed by
delta2prop().
prop2delta(x = NULL, Rstandard = NULL)prop2delta(x = NULL, Rstandard = NULL)
x |
Vector of proportions values. |
Rstandard |
String describing the isotopic measurement, e.g. "d15N", "d13C" and used to set automatically Rstandards (see the Section "Ratios for reference standards" for more details). Alternatively, a numeric value to use for Rstandard, e.g. 0.0036765. |
A vector of same length of x, containing the delta values based on the proportions of heavy isotope provided as x and the Rstandard provided.
prop15N <- c(0.00395, 0.02222, 0.00462, 0.00753, NA, 0.00422, 0.00492) # Rstandard can be specified with a string for some preset references d15N <- prop2delta(prop15N, "d15N") d15N # Rstandard can also be specified manually for non-preset references d15N_manual <- prop2delta(prop15N, 0.0036765) d15N_manual # Call delta2prop() to get the detail of available references delta2prop()prop15N <- c(0.00395, 0.02222, 0.00462, 0.00753, NA, 0.00422, 0.00492) # Rstandard can be specified with a string for some preset references d15N <- prop2delta(prop15N, "d15N") d15N # Rstandard can also be specified manually for non-preset references d15N_manual <- prop2delta(prop15N, 0.0036765) d15N_manual # Call delta2prop() to get the detail of available references delta2prop()
Draw a Sankey plot with basic defaults
quick_sankey(flows, ...)quick_sankey(flows, ...)
flows |
A tibble containing flows (output from
|
... |
Passed to |
Mostly called for its side effect (plotting), but also returns invisible the scene object describing the Sankey plot. Note that the structure of this object is experimental and might change in the future!
Run a MCMC sampler on a network model using Stan
run_mcmc( model, iter = 2000, chains = 4, method = "matrix_exp", euler_control = list(), cores = NULL, stanfit = FALSE, vb = FALSE, ... )run_mcmc( model, iter = 2000, chains = 4, method = "matrix_exp", euler_control = list(), cores = NULL, stanfit = FALSE, vb = FALSE, ... )
model |
A |
iter |
A positive integer specifying the number of iterations for each chain (including warmup). The default is 2000. |
chains |
A positive integer specifying the number of Markov chains. The default is 4. |
method |
A character string indicating the method to use to solve ODE
in the Stan model; available methods are "matrix_exp" and "euler". The
default is "matrix_exp", which uses matrix exponential and is reasonably
fast for small networks. For large networks, the "euler" method can be
used. It implements a simple forward Euler method to solve the ODE and can
be faster than the matrix exponential approach, but extra caution must be
taken to check for numerical accuracy (e.g. testing different |
euler_control |
An optional list containing extra parameters when using
|
cores |
Number of cores to use for parallel use. Default is
|
stanfit |
If TRUE, returns a 'stanfit' object instead of the more classical 'mcmc.list' object. Note that when an 'mcmc.list' object is returned, the original 'stanfit' object is still accessible as an attribute of that object (see Examples). |
vb |
Boolean, if TRUE will use |
... |
Arguments passed to 'rstan::sampling' (e.g. iter, chains). |
An object of class 'stanfit' returned by 'rstan::sampling' if
stanfit = TRUE, otherwise the result of converting this
stanfit object with stanfit_to_named_mcmclist (i.e. an object
of class networkModelStanfit and mcmc.list, which still
carries the original 'stanfit' object stored as an attribute).
aquarium_mod ## Not run: # The 'aquarium_run' object is shipped with the package, so you don't # actually need to run the line below to obtain it aquarium_run <- run_mcmc(aquarium_mod) plot(aquarium_run) summary(aquarium_run) # The original stanfit object returned by Stan sfit <- attr(aquarium_run, "stanfit") sfit # The stanfit object can be used for diagnostics, LOO cross-validation, etc. rstan::loo(sfit) ## End(Not run)aquarium_mod ## Not run: # The 'aquarium_run' object is shipped with the package, so you don't # actually need to run the line below to obtain it aquarium_run <- run_mcmc(aquarium_mod) plot(aquarium_run) summary(aquarium_run) # The original stanfit object returned by Stan sfit <- attr(aquarium_run, "stanfit") sfit # The stanfit object can be used for diagnostics, LOO cross-validation, etc. rstan::loo(sfit) ## End(Not run)
Generate samples from a network model
sample_from( nm, at, dt = NULL, grid_size = NULL, end = NULL, error.draws = 1, cached_ts = NULL, cached_ee = NULL )sample_from( nm, at, dt = NULL, grid_size = NULL, end = NULL, error.draws = 1, cached_ts = NULL, cached_ee = NULL )
nm |
A |
at |
Vector of time values at which the samples should be taken. |
dt, grid_size
|
Time step size or grid points, respectively. |
end |
Final timepoint used in the projections. |
error.draws |
Integer, number of draws from the error distribution for each sample (default: 1). |
cached_ts, cached_ee
|
Used for optimization by other functions, not for use by the package user. |
A tibble containing the generated samples.
library(magrittr) mod <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia -> NH4") inits <- tibble::tribble( ~comps, ~sizes, ~props, ~treatment, "NH4", 0.2, 0.8, "light", "algae", 1, 0.004, "light", "daphnia", 2, 0.004, "light", "NH4", 0.5, 0.8, "dark", "algae", 1.2, 0.004, "dark", "daphnia", 1.3, 0.004, "dark") mod <- set_init(mod, inits, comp = "comps", size = "sizes", prop = "props", group_by = "treatment") mod <- add_covariates(mod, upsilon_NH4_to_algae ~ treatment) mod <- mod %>% set_params(c("eta" = 0.2, "lambda_algae" = 0, "lambda_daphnia" = 0, "lambda_NH4" = 0, "upsilon_NH4_to_algae|light" = 0.3, "upsilon_NH4_to_algae|dark" = 0.1, "upsilon_algae_to_daphnia" = 0.13, "upsilon_daphnia_to_NH4" = 0.045, "zeta" = 0.1)) spl <- mod %>% sample_from(at = 1:10) spllibrary(magrittr) mod <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia -> NH4") inits <- tibble::tribble( ~comps, ~sizes, ~props, ~treatment, "NH4", 0.2, 0.8, "light", "algae", 1, 0.004, "light", "daphnia", 2, 0.004, "light", "NH4", 0.5, 0.8, "dark", "algae", 1.2, 0.004, "dark", "daphnia", 1.3, 0.004, "dark") mod <- set_init(mod, inits, comp = "comps", size = "sizes", prop = "props", group_by = "treatment") mod <- add_covariates(mod, upsilon_NH4_to_algae ~ treatment) mod <- mod %>% set_params(c("eta" = 0.2, "lambda_algae" = 0, "lambda_daphnia" = 0, "lambda_NH4" = 0, "upsilon_NH4_to_algae|light" = 0.3, "upsilon_NH4_to_algae|dark" = 0.1, "upsilon_algae_to_daphnia" = 0.13, "upsilon_daphnia_to_NH4" = 0.045, "zeta" = 0.1)) spl <- mod %>% sample_from(at = 1:10) spl
Sample from a prior object
sample_from_prior(x, n = 1)sample_from_prior(x, n = 1)
x |
A |
n |
Integer, number of samples to draw. |
A numeric vector of length n.
sample_from_prior(constant_p(1)) sample_from_prior(constant_p(1), 10) sample_from_prior(hcauchy_p(0.5), 1) hist(sample_from_prior(hcauchy_p(0.5), 20)) hist(sample_from_prior(uniform_p(0, 3), 1000)) hist(sample_from_prior(scaled_beta_p(3, 7, 2), 1e4))sample_from_prior(constant_p(1)) sample_from_prior(constant_p(1), 10) sample_from_prior(hcauchy_p(0.5), 1) hist(sample_from_prior(hcauchy_p(0.5), 20)) hist(sample_from_prior(uniform_p(0, 3), 1000)) hist(sample_from_prior(scaled_beta_p(3, 7, 2), 1e4))
Sample parameter values from priors
sample_params(nm)sample_params(nm)
nm |
A |
A named vector containing parameter values.
library(magrittr) p <- sample_params(aquarium_mod) p proj <- aquarium_mod %>% set_params(p) %>% project(end = 10) plot(proj)library(magrittr) p <- sample_params(aquarium_mod) p proj <- aquarium_mod %>% set_params(p) %>% project(end = 10) plot(proj)
Draw a Sankey plot for a network and estimated flows
sankey( topo, nodes = NULL, flows = NULL, layout = NULL, new = TRUE, debug = FALSE, node_f = 1, edge_f = 1, node_s = "auto", edge_n = 32, cex_lab = NULL, cex.lab = NULL, fit = TRUE )sankey( topo, nodes = NULL, flows = NULL, layout = NULL, new = TRUE, debug = FALSE, node_f = 1, edge_f = 1, node_s = "auto", edge_n = 32, cex_lab = NULL, cex.lab = NULL, fit = TRUE )
topo |
A topology. |
nodes |
Optional, a tibble containing the properties of the nodes. It should have a 'comp' column with the same entries as the topology. It cannot have 'x' and 'y' entries. If it has a 'label' entry, it will replace the 'comp' values for node labels. |
flows |
A tibble containing the values of the flows in the topology. If NULL (the default), all flows have same width in the plot. |
layout |
String, node-placing algorithm to use from the ggraph package
(e.g. "stress"). The ggraph package itself uses some algoritms from the
igraph package. See the Details in the help of
|
new |
Boolean, create a new page for the plot? |
debug |
Boolean, if TRUE then draw a lot of shapes to help with debugging. |
node_f, edge_f
|
Multiplicative factor to adjust node and edge size. |
node_s |
String defining how node size is calculated. The effect of the string also depends on the chosen layout. |
edge_n |
Integer, number of interpolation points along each edge. |
cex_lab, cex.lab
|
Expansion factor for label size (both arguments are synonyms). |
fit |
Boolean, if TRUE try to fit all the graphical elements inside the canvas. |
Mostly called for its side effect (plotting), but also returns invisible the scene object describing the Sankey plot. Note that the structure of this object is experimental and might change in the future!
library(magrittr) topo <- topo(trini_mod) sankey(topo, debug = TRUE) sankey(topo, layout = "stress") sankey(topo(aquarium_mod), layout = "stress", edge_f = 0.5) m <- new_networkModel() %>% set_topo(c("subs -> NH3 -> subs", "NH3 -> Q, E", "E -> Q -> E", "E -> D, M")) %>% set_steady("subs") %>% set_prop_family("normal_sd") ggtopo(m) sankey(topo(m), layout = "stress") # Debug visualization ## Helper functions flows_from_topo <- function(x) { x <- unclass(x) # Remove the "topo" class to treat it as a matrix n_comps <- ncol(x) links <- which(x > 0) from <- links %/% n_comps + 1 to <- links %% n_comps links <- tibble::tibble(from = from, to = to) for (i in seq_len(nrow(links))) { if (links$to[i] == 0) { links$from[i] <- links$from[i] - 1 links$to[i] <- n_comps } stopifnot(x[links$to[i], links$from[i]] > 0) } flows <- tibble::tibble(from = colnames(x)[links$from], to = rownames(x)[links$to]) return(flows) } nodes_from_topo <- function(x) { nodes <- tibble::tibble(comp = colnames(x), label = colnames(x)) return(nodes) } t <- topo(trini_mod) nodes <- nodes_from_topo(t) nodes$label <- as.list(nodes$label) nodes$label[[2]] <- latex2exp::TeX("$\\beta$") nodes$size <- runif(nrow(nodes), 1, 2) flows <- flows_from_topo(t) flows$width <- runif(nrow(flows), 0.2, 2) z <- sankey(t, nodes = nodes, flows = flows, layout = "left2right", debug = TRUE, node_f = 1, edge_f = 0.9, edge_n = 32, cex_lab = 1.5) # Stress layout y <- new_networkModel() %>% set_topo(c("subs -> NH3 -> subs", "NH3 -> Q, E", "E -> Q -> E", "E -> D, M")) %>% set_steady("subs") %>% set_prop_family("normal_sd") y <- topo(y) nodes <- nodes_from_topo(y) nodes$size <- runif(nrow(nodes), 1, 10) ggtopo(y, edge = "fan") flows <- flows_from_topo(y) flows$width <- runif(nrow(flows), 0.2, 5) z <- sankey(y, nodes = nodes, flows = flows, debug = FALSE, edge_n = 32, edge_f = 0.4, node_s = "prop") # Another example r <- new_networkModel() %>% set_topo("infusion -> plasma -> body -> plasma") %>% set_steady(c("infusion", "body")) r <- topo(r) ggtopo(r, edge = "fan") sankey(r, debug = TRUE, edge_f = 0.2)library(magrittr) topo <- topo(trini_mod) sankey(topo, debug = TRUE) sankey(topo, layout = "stress") sankey(topo(aquarium_mod), layout = "stress", edge_f = 0.5) m <- new_networkModel() %>% set_topo(c("subs -> NH3 -> subs", "NH3 -> Q, E", "E -> Q -> E", "E -> D, M")) %>% set_steady("subs") %>% set_prop_family("normal_sd") ggtopo(m) sankey(topo(m), layout = "stress") # Debug visualization ## Helper functions flows_from_topo <- function(x) { x <- unclass(x) # Remove the "topo" class to treat it as a matrix n_comps <- ncol(x) links <- which(x > 0) from <- links %/% n_comps + 1 to <- links %% n_comps links <- tibble::tibble(from = from, to = to) for (i in seq_len(nrow(links))) { if (links$to[i] == 0) { links$from[i] <- links$from[i] - 1 links$to[i] <- n_comps } stopifnot(x[links$to[i], links$from[i]] > 0) } flows <- tibble::tibble(from = colnames(x)[links$from], to = rownames(x)[links$to]) return(flows) } nodes_from_topo <- function(x) { nodes <- tibble::tibble(comp = colnames(x), label = colnames(x)) return(nodes) } t <- topo(trini_mod) nodes <- nodes_from_topo(t) nodes$label <- as.list(nodes$label) nodes$label[[2]] <- latex2exp::TeX("$\\beta$") nodes$size <- runif(nrow(nodes), 1, 2) flows <- flows_from_topo(t) flows$width <- runif(nrow(flows), 0.2, 2) z <- sankey(t, nodes = nodes, flows = flows, layout = "left2right", debug = TRUE, node_f = 1, edge_f = 0.9, edge_n = 32, cex_lab = 1.5) # Stress layout y <- new_networkModel() %>% set_topo(c("subs -> NH3 -> subs", "NH3 -> Q, E", "E -> Q -> E", "E -> D, M")) %>% set_steady("subs") %>% set_prop_family("normal_sd") y <- topo(y) nodes <- nodes_from_topo(y) nodes$size <- runif(nrow(nodes), 1, 10) ggtopo(y, edge = "fan") flows <- flows_from_topo(y) flows$width <- runif(nrow(flows), 0.2, 5) z <- sankey(y, nodes = nodes, flows = flows, debug = FALSE, edge_n = 32, edge_f = 0.4, node_s = "prop") # Another example r <- new_networkModel() %>% set_topo("infusion -> plasma -> body -> plasma") %>% set_steady(c("infusion", "body")) r <- topo(r) ggtopo(r, edge = "fan") sankey(r, debug = TRUE, edge_f = 0.2)
If a random variable X follows a scaled beta distribution with parameters (alpha, beta, scale), then X/scale follows a beta distribution with parameters (alpha, beta).
scaled_beta_p(alpha, beta, scale = 1)scaled_beta_p(alpha, beta, scale = 1)
alpha |
Alpha parameter of the unscaled beta distribution. |
beta |
Beta parameter of the unscaled beta distribution. |
scale |
The upper boundary of the prior. |
A list defining the prior.
scaled_beta_p(0.8, 20, scale = 10)scaled_beta_p(0.8, 20, scale = 10)
Select parameters based on their names
## S3 method for class 'mcmc.list' select(.data, ...)## S3 method for class 'mcmc.list' select(.data, ...)
.data |
A |
... |
Strings used to select variables using pattern matching with
|
An mcmc.list object, with the same extra class(es) as
.data (if any).
Indicating a non-zero value for half-life will add a decay to the marked portion of the tracer element. The decay constant is calculated from the half-life value as:
set_half_life(nm, hl, quiet = FALSE)set_half_life(nm, hl, quiet = FALSE)
nm |
A |
hl |
Half-life value, in the same time unit as the observations are (or will be) given. Setting half-life to zero is equivalent to using a stable isotope (no decay used in the model). |
quiet |
Boolean for verbosity. |
lambda_decay = log(2) / half_life
Note that for correct calculations the half-life value should be given in the same time unit (e.g. hour, day) that the time unit used for observations.
A networkModel object.
library(magrittr) x <- new_networkModel() %>% set_topo("32P -> root -> leaf") %>% set_half_life(hl = 14.268) xlibrary(magrittr) x <- new_networkModel() %>% set_topo("32P -> root -> leaf") %>% set_half_life(hl = 14.268) x
Set initial conditions in a network model
set_init(nm, data, comp, size, prop, group_by = NULL)set_init(nm, data, comp, size, prop, group_by = NULL)
nm |
A |
data |
A tibble containing the initial conditions |
comp |
String, name of the |
size |
String, name of the |
prop |
String, name of the |
group_by |
Optional vector of string giving the names of the columns to use for grouping the data into replicates |
A networkModel object.
# Using the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph", "FBOM -> tricor", "petro, tricor -> arg") # Taking initial condtions from the 'lalaja' dataset at t=0 inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m# Using the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph", "FBOM -> tricor", "petro, tricor -> arg") # Taking initial condtions from the 'lalaja' dataset at t=0 inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m
Set observations in a network model
set_obs(nm, data, comp, size, prop, time, group_by)set_obs(nm, data, comp, size, prop, time, group_by)
nm |
A |
data |
A tibble containing the observations. If NULL, remove observations from the model. |
comp |
String, name of the |
size |
String, name of the |
prop |
String, name of the |
time |
String, name of the |
group_by |
Optional vector of string giving the names of the columns to use for grouping the data into replicates |
A networkModel object.
# Using the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph", "FBOM -> tricor", "petro, tricor -> arg") # Taking initial condtions from the 'lalaja' dataset at t=0 inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m # Taking observations from 'lalaja' m <- set_obs(m, lalaja[lalaja[["time.days"]] > 0, ], time = "time.days") m plot(m)# Using the topology from the Trinidad case study m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi, FBOM", "epi -> petro, pseph", "FBOM -> tricor", "petro, tricor -> arg") # Taking initial condtions from the 'lalaja' dataset at t=0 inits <- lalaja[lalaja[["time.days"]] == 0, ] inits m <- set_init(m, inits, comp = "compartment", size = "mgN.per.m2", prop = "prop15N", group_by = "transect") m # Taking observations from 'lalaja' m <- set_obs(m, lalaja[lalaja[["time.days"]] > 0, ], time = "time.days") m plot(m)
Set the parameters in a network model
set_params(nm, params, force = TRUE, quick = FALSE)set_params(nm, params, force = TRUE, quick = FALSE)
nm |
A |
params |
A named vector or a tibble with columns c("parameter", "value") containing the (global) parameter values. |
force |
Boolean, if FALSE will not overwrite already set parameters. |
quick |
Boolean, if TRUE take some shortcuts for faster parameter settings when called by another function. This should usually be left to the default (FALSE) by a regular package user. |
A networkModel object.
m <- aquarium_mod p <- sample_params(m) m2 <- set_params(m, p) m2$parametersm <- aquarium_mod p <- sample_params(m) m2 <- set_params(m, p) m2$parameters
Set prior(s) for a network model
set_prior(x, prior, param = "", use_regexp = TRUE, quiet = FALSE) set_priors(x, prior, param = "", use_regexp = TRUE, quiet = FALSE)set_prior(x, prior, param = "", use_regexp = TRUE, quiet = FALSE) set_priors(x, prior, param = "", use_regexp = TRUE, quiet = FALSE)
x |
A |
prior |
A prior built with e.g. uniform_p() or hcauchy_p(). Call
|
param |
String, target parameter or regexp to target several
parameters. Default is the empty string |
use_regexp |
Boolean, if |
quiet |
Boolean, if |
A networkModel object.
# Copy `aquarium_mod` m <- aquarium_mod priors(m) # Modify the priors of `m` m <- set_priors(m, exponential_p(0.5), "lambda") priors(m) # Re-apply priors from the original `aquarium_mod` prev_priors <- priors(aquarium_mod) prev_priors m <- set_priors(m, prev_priors) priors(m)# Copy `aquarium_mod` m <- aquarium_mod priors(m) # Modify the priors of `m` m <- set_priors(m, exponential_p(0.5), "lambda") priors(m) # Re-apply priors from the original `aquarium_mod` prev_priors <- priors(aquarium_mod) prev_priors m <- set_priors(m, prev_priors) priors(m)
Set the distribution family for observed proportions
set_prop_family(nm, family, quiet = FALSE)set_prop_family(nm, family, quiet = FALSE)
nm |
A |
family |
Allowed values are "gamma_cv", "beta_phi", "normal_cv", and "normal_sd". |
quiet |
Boolean, if |
A networkModel object.
library(magrittr) m <- new_networkModel() %>% set_topo(links = "NH4, NO3 -> epi -> pseph, tricor") m <- m %>% set_prop_family("beta_phi") m attr(m, "prop_family")library(magrittr) m <- new_networkModel() %>% set_topo(links = "NH4, NO3 -> epi -> pseph, tricor") m <- m %>% set_prop_family("beta_phi") m attr(m, "prop_family")
Set the distribution family for observed sizes
set_size_family(nm, family, by_compartment, quiet = FALSE, quiet_reset = FALSE)set_size_family(nm, family, by_compartment, quiet = FALSE, quiet_reset = FALSE)
nm |
A |
family |
Allowed values are "normal_cv" and "normal_sd". |
by_compartment |
Boolean, if TRUE then zeta is compartment-specific. |
quiet |
Boolean, if |
quiet_reset |
Boolean, write a message when model parameters (and covariates and priors) are reset? |
A networkModel object.
library(magrittr) m <- new_networkModel() %>% set_topo(links = "NH4, NO3 -> epi -> pseph, tricor") m <- m %>% set_size_family("normal_sd") m attr(m, "size_family") m <- m %>% set_size_family(by_compartment = TRUE) attr(m, "size_zeta_per_compartment")library(magrittr) m <- new_networkModel() %>% set_topo(links = "NH4, NO3 -> epi -> pseph, tricor") m <- m %>% set_size_family("normal_sd") m attr(m, "size_family") m <- m %>% set_size_family(by_compartment = TRUE) attr(m, "size_zeta_per_compartment")
This function automatically adds a default prior (uniform on [0,1]) for the active portion of split compartments.
set_split(nm, comps = NULL, which = NULL)set_split(nm, comps = NULL, which = NULL)
nm |
A |
comps |
Vector of strings, the names of the compartments to set split. |
which |
Vector of integers giving the nm rows to update. Default is to update all rows. |
A networkModel object.
library(magrittr) x <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia") %>% set_split("algae") topo(x)library(magrittr) x <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia") %>% set_split("algae") topo(x)
Flag some network compartments as being in a steady state
set_steady(nm, comps = NULL, which = NULL)set_steady(nm, comps = NULL, which = NULL)
nm |
A |
comps |
Vector of strings, names of the compartments to set steady. |
which |
Vector of integers giving the nm rows to update. Default is to update all rows. |
A networkModel object.
library(magrittr) x <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia") %>% set_steady("NH4") topo(x)library(magrittr) x <- new_networkModel() %>% set_topo("NH4 -> algae -> daphnia") %>% set_steady("NH4") topo(x)
Set the topology in a network model.
set_topo(nm, ..., from = NULL, to = NULL)set_topo(nm, ..., from = NULL, to = NULL)
nm |
A |
... |
One or more strings describing the links defining the network topology. Optionally, links can be given as a data frame. See the examples for more details about acceptable input formats. |
from |
Optional, string containing the column name for sources if links are provided as a data frame. |
to |
Optional, string containing the column name for destinations if links are provided as a data frame. |
A networkModel object.
# A single string can describe several links in one go. m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi -> pseph, tricor") m topo(m) # Several strings can be given as distinct arguments. m2 <- new_networkModel() %>% set_topo("NH4, NO3 -> epi -> pseph, tricor", "NH4 -> FBOM, CBOM", "CBOM <- NO3") m2 topo(m2) # Multiple strings can be also be combined into a single argument with `c()`. links <- c("NH4, NO3 -> epi -> pseph, tricor", "NH4 -> FBOM, CBOM", "CBOM <- NO3") m3 <- new_networkModel() %>% set_topo(links) m3 topo(m3) # A data frame can be used to specify the links. links <- data.frame(source = c("NH4", "NO3", "epi"), consumer = c("epi", "epi", "petro")) links m4 <- new_networkModel() %>% set_topo(links, from = "source", to = "consumer") m4 m4$topology[[1]]# A single string can describe several links in one go. m <- new_networkModel() %>% set_topo("NH4, NO3 -> epi -> pseph, tricor") m topo(m) # Several strings can be given as distinct arguments. m2 <- new_networkModel() %>% set_topo("NH4, NO3 -> epi -> pseph, tricor", "NH4 -> FBOM, CBOM", "CBOM <- NO3") m2 topo(m2) # Multiple strings can be also be combined into a single argument with `c()`. links <- c("NH4, NO3 -> epi -> pseph, tricor", "NH4 -> FBOM, CBOM", "CBOM <- NO3") m3 <- new_networkModel() %>% set_topo(links) m3 topo(m3) # A data frame can be used to specify the links. links <- data.frame(source = c("NH4", "NO3", "epi"), consumer = c("epi", "epi", "petro")) links m4 <- new_networkModel() %>% set_topo(links, from = "source", to = "consumer") m4 m4$topology[[1]]
Return the distribution family for observed sizes
size_family(nm, quiet = FALSE)size_family(nm, quiet = FALSE)
nm |
A |
quiet |
Boolean for being quiet about explaining the role of zeta
(default is |
A character string describing the distribution family used to model observed sizes.
size_family(aquarium_mod) size_family(trini_mod)size_family(aquarium_mod) size_family(trini_mod)
When running run_mcmc with stanfit = FALSE (typically for
debugging purposes), the parameters in the returned stanfit object
are named using a base label and an indexing system. This function provides
a way to convert this stanfit object into a more conventional
mcmc.list object where parameters are named according to their role
in the original network model used when running run_mcmc.
stanfit_to_named_mcmclist(stanfit)stanfit_to_named_mcmclist(stanfit)
stanfit |
A stanfit object returned by |
An mcmc.list object. It also has the original stanfit object
stored as an attribute "stanfit".
Extract data from a networkModel object into a tidy tibble.
tidy_data(x)tidy_data(x)
x |
A networkModel object. |
A tibble (note: row ordering is not the same as in the input).
tidy_data(aquarium_mod) tidy_data(trini_mod)tidy_data(aquarium_mod) tidy_data(trini_mod)
This function prepares both tidy data from a model and tidy posterior predictions from a model fit. Having those two tibbles prepared at the same time allows to merge them to ensure that observed data, predicted data and original variables other than observations are all in sync when using y and y_rep objects for bayesplot functions.
tidy_dpp(model, fit, draw = NULL, cores = NULL)tidy_dpp(model, fit, draw = NULL, cores = NULL)
model |
A networkModel object. |
fit |
A networkModelStanfit object. |
draw |
Integer, number of draws to sample from the posterior. |
cores |
Number of cores to use for parallel calculations. Default is
|
A list with y, y_rep and vars.
If neither n_per_chain and n are provided, all iterations are
used.
tidy_flows( nm, mcmc, n_per_chain = NULL, n = NULL, n_grid = 64, steady_state = FALSE, dt = NULL, grid_size = NULL, at = NULL, end = NULL, use_cache = TRUE, cores = NULL )tidy_flows( nm, mcmc, n_per_chain = NULL, n = NULL, n_grid = 64, steady_state = FALSE, dt = NULL, grid_size = NULL, at = NULL, end = NULL, use_cache = TRUE, cores = NULL )
nm |
A |
mcmc |
The corresponding output from |
n_per_chain |
Integer, number of iterations randomly drawn per chain. Note that iterations are in sync across chains (in practice, random iterations are chosen, and then parameter values extracted for those same iterations from all chains). |
n |
Integer, number of iterations randomly drawn from |
n_grid |
Size of the time grid used to calculate trajectories |
steady_state |
Boolean (default: FALSE). If TRUE, then steady state
compartment sizes are calculated for each iteration and steady state flows
are calculated from those compartment sizes. Note that any pulse that
might be specified in the input model |
dt, grid_size
|
Time step size or grid points, respectively. |
at |
Timepoints at which the predictions should be returned. |
end |
Final timepoint used in the projections. |
use_cache |
Boolean, use cache for faster calculations? |
cores |
Number of cores to use for parallel calculations. Default is
|
Warning: This function is still maturing and its interface and output might change in the future.
Note about how steady state sizes for split compartments are calculated: the steady size of the active portion is calculated divide it is divided by the active fraction (portion.act parameter) to get the total size including the refractory portion. In this case we get a "steady-state" refractory portion, consistent with steady state size of active fraction and with portion.act parameter.
A tidy table containing the mcmc iterations (chain, iteration,
parameters), the grouping variables from the network model and the
flows. The returned flow values are the average flow per unit of time
over the trajectory calculations (or steady state flows if
steady_state is TRUE).
tf <- tidy_flows(aquarium_mod, aquarium_run, n_per_chain = 25, cores = 2) tf tfmcmc <- as.mcmc.list(tf) plot(tfmcmc)tf <- tidy_flows(aquarium_mod, aquarium_run, n_per_chain = 25, cores = 2) tf tfmcmc <- as.mcmc.list(tf) plot(tfmcmc)
Extract a tidy output from an mcmc.list
tidy_mcmc(x, spread = FALSE, include_constant = TRUE)tidy_mcmc(x, spread = FALSE, include_constant = TRUE)
x |
An mcmc.list object |
spread |
Boolean, spread the parameters into separate columns? |
include_constant |
Boolean, include constant parameters as proper parameter traces? |
A tidy table containing one iteration per row
fit <- lapply(1:4, function(i) { z <- matrix(rnorm(200), ncol = 2) colnames(z) <- c("alpha", "beta") coda::as.mcmc(z) }) fit <- coda::as.mcmc.list(fit) tidy_mcmc(fit) tidy_mcmc(fit, spread = TRUE)fit <- lapply(1:4, function(i) { z <- matrix(rnorm(200), ncol = 2) colnames(z) <- c("alpha", "beta") coda::as.mcmc(z) }) fit <- coda::as.mcmc.list(fit) tidy_mcmc(fit) tidy_mcmc(fit, spread = TRUE)
Draw from the posterior predictive distribution of the model outcome
tidy_posterior_predict(object, newdata, draw = NULL, cores = NULL, ...)tidy_posterior_predict(object, newdata, draw = NULL, cores = NULL, ...)
object |
A networkModelStanfit object. |
newdata |
The original model used to fit the networkStanfit object. |
draw |
Integer, number of draws to sample from the posterior. Default is 100. |
cores |
Number of cores to use for parallel calculations. Default is
|
... |
Not used for now. |
A tidy table.
If neither n_per_chain and n are provided, all iterations are
used.
tidy_steady_states(nm, mcmc, n_per_chain = NULL, n = NULL)tidy_steady_states(nm, mcmc, n_per_chain = NULL, n = NULL)
nm |
A |
mcmc |
The corresponding output from |
n_per_chain |
Integer, number of iterations randomly drawn per chain. Note that iterations are in sync across chains (in practice, random iterations are chosen, and then parameter values extracted for those same iterations from all chains). |
n |
Integer, number of iterations randomly drawn from |
Note about how steady state sizes for split compartments are calculated: the steady size of the active portion is calculated divide it is divided by the active fraction (portion.act parameter) to get the total size including the refractory portion. In this case we get a "steady-state" refractory portion, consistent with steady state size of active fraction and with portion.act parameter.
A tidy table containing the mcmc iterations (chain, iteration, parameters), the grouping variables from the network model and the steady state sizes.
If neither n_per_chain and n are provided, all iterations are
used.
tidy_trajectories( nm, mcmc, n_per_chain = NULL, n = NULL, n_grid = 64, dt = NULL, grid_size = NULL, at = NULL, end = NULL, use_cache = TRUE, cores = NULL )tidy_trajectories( nm, mcmc, n_per_chain = NULL, n = NULL, n_grid = 64, dt = NULL, grid_size = NULL, at = NULL, end = NULL, use_cache = TRUE, cores = NULL )
nm |
A |
mcmc |
The corresponding output from |
n_per_chain |
Integer, number of iterations randomly drawn per chain. Note that iterations are in sync across chains (in practice, random iterations are chosen, and then parameter values extracted for those same iterations from all chains). |
n |
Integer, number of iterations randomly drawn from |
n_grid |
Size of the time grid used to calculate trajectories |
dt, grid_size
|
Time step size or grid points, respectively. |
at |
Timepoints at which the predictions should be returned. |
end |
Final timepoint used in the projections. |
use_cache |
Boolean, use cache for faster calculations? |
cores |
Number of cores to use for parallel calculations. Default is
|
Warning: This function is still maturing and its interface and output might change in the future.
A tidy table containing the mcmc iterations (chain, iteration, parameters), the grouping variables from the network model and the trajectories.
tt <- tidy_trajectories(aquarium_mod, aquarium_run, n = 10, cores = 2) tttt <- tidy_trajectories(aquarium_mod, aquarium_run, n = 10, cores = 2) tt
Return the list of topologies, or a unique topology if all identical
topo(nm, simplify = TRUE)topo(nm, simplify = TRUE)
nm |
A |
simplify |
Boolean, return only a unique topology if all topologies are identical or if there is only one? Default is TRUE. |
A list of the networkModel topologies or, if all topologies
are identical (or if there is only one) and simplify is TRUE, a
single topology (not wrapped into a single-element list).
aquarium_mod topo(aquarium_mod) trini_mod topo(trini_mod)aquarium_mod topo(aquarium_mod) trini_mod topo(trini_mod)
Plot mcmc.list objects
traceplot(x, ...)traceplot(x, ...)
x |
A |
... |
Passed to |
Called for side effect (plotting).
This model is used in the package case study about Trinidadian streams and is based on an original dataset taken from Collins et al. (2016).
trini_modtrini_mod
An object of class networkModel (inherits from tbl_df, tbl, data.frame) with 6 rows and 6 columns.
The model is complete, with topology, initial conditions, observations, covariates and priors.
It is ready for an MCMC run as shown in the example. Note that it might be a good idea to relax the priors for uptake rates from seston to Leptonema (e.g. using hcauchy_p(10)), seston being a compartment that is flowing with the stream water and that can be replenished from upstream.
This network model contains data from the original article: Collins, Sarah M., Steven A. Thomas, Thomas Heatherly, Keeley L. MacNeill, Antoine O.H.C. Leduc, Andrés López-Sepulcre, Bradley A. Lamphere, et al. 2016. “Fish Introductions and Light Modulate Food Web Fluxes in Tropical Streams: A Whole-Ecosystem Experimental Approach.” Ecology, <doi:10.1002/ecy.1530>.
This dataset was also used in the paper: López-Sepulcre, Andrés, Matthieu Bruneaux, Sarah M. Collins, Rana El-Sabaawi, Alexander S. Flecker, and Steven A. Thomas. 2020. “A New Method to Reconstruct Quantitative Food Webs and Nutrient Flows from Isotope Tracer Addition Experiments.” The American Naturalist 195 (6): 964–85. <doi:10.1086/708546>.
trini_mod ggtopo(trini_mod) ## Not run: # Warning: the run below can take quite a long time! # (about 15 min with 4 cores at 3.3 Ghz). run <- run_mcmc(trini_mod, iter = 500, chains = 4, cores = 4) ## End(Not run)trini_mod ggtopo(trini_mod) ## Not run: # Warning: the run below can take quite a long time! # (about 15 min with 4 cores at 3.3 Ghz). run <- run_mcmc(trini_mod, iter = 500, chains = 4, cores = 4) ## End(Not run)
prior object in tibblesFunction used for displaying prior object in tibbles
## S3 method for class 'prior' type_sum(x)## S3 method for class 'prior' type_sum(x)
x |
An object of class |
Input formatted with format(x).
Define a uniform prior
uniform_p(min, max)uniform_p(min, max)
min, max
|
Minimum and maximum boundaries for the uniform prior. |
A list defining the prior.
uniform_p(min = 0, max= 1)uniform_p(min = 0, max= 1)