Nothing
R/stability.R
In fluxweb: Estimate Energy Fluxes in Food Webs
#' Estimates network stability
#'
#' Computes resilience of the system through Jacobian matrix eigenvalues.
#'
#'
#' @param val.mat A matrix describing fluxes between species (usually a result of \code{\link[fluxweb]{fluxing}} function).
#' @param biomasses A vector of species biomasses.
#' @param losses A vector or an array of species energy losses (excluding predation).
#' @param efficiencies A vector or an array of conversion efficiencies of species in the adjacency matrix. These values describe the proportion of consumed energy that is converted to biomass of the consumer.
#' @param growth.rate A vector defining growth rate of basal species.
#' @param bioms.prefs Logical, if \code{TRUE} (default) preferences are scaled according to species biomasses.
#' @param bioms.losses Logical, if \code{TRUE} (default) losses are scaled with biomass.
#' @param ef.level Set to \code{"prey"} if efficiencies are defined by prey, \code{"pred"} if they are a property of the predator.
#' @param full.output Logical, if \code{TRUE} function return supplementary informations.
#'
#' @return Maximum eigenvalue of the Jacobian matrix of a Lotka Voltera like system of equations. If full.output, Jacobian eigenvalues and eigenvectors are returned.
#'
#'
#'@details
#'
#'
#For more information about the set of the underlying system of differential equations and mathematical derivation, please read Gauzens et al. 2017 SI2. <https://www.biorxiv.org/content/early/2017/12/06/229450>
#'
#'
#'\itemize{
#'
#'\item{\code{losses}:} Express species energetic losses not related to consumption. Usually metabolic or death rates.
#'When an array is provided, losses associated to each species correspond to line sums.
#'
#'\item{\code{efficiencies}:} Determines how efficient species are to convert energy (see \code{ef.level} for more details).
#'Providing an array will assume values depending on both prey and predator identity.
#'
#'\item{\code{growth.rate}:} Growth rates of basal species defined. Length of the vector should be equal to the number of species.
#'expects positive numeric values for index corresponding to basal species, NA otherwise
#'
#'\item{\code{bioms.pref}:} If \code{TRUE}, preferences \eqn{w_{ij}} of predator j on prey i are scaled according to species biomass using the following formula:
#'\deqn{
#'w_{i,j} = \frac{mat[i,j] * biomasses[i]}{\sum_k mat[i,k]* biomasses[k]}
#'}
#'\item{\code{bioms.losses}:} If \code{TRUE}, function will assume that losses are defined per biomass unit.
#'Thus, total losses will be thereafter multiplied by biomass values for each species.
#'
#'\item{\code{ef.level}:} If \code{"prey"} (resp \code{"pred"}), the total amount of energy that can be metabolized from a trophic link
#'will be determined by prey (resp pred) identity. \code{"link.specific"} assumes that efficiencies are defined for each trophic interaction
#'and implies \code{efficiencies} parameter to be a matrix
#'
#'\item{\code{full.output}:} If \code{TRUE}, function result is a list of eigenvalues and eigenvectors of the Jacobian matrix.
#'}
#'
#' @examples
# # first compute species per unit biomass metabolic rates using the metabolic theory:
#' losses = 0.15 * groups.level$bodymasses^(-0.25)
#'
#' # growth rates of basal sppecies
#' growth.rates = rep(NA, dim(groups.level$mat)[1])
#' growth.rates[colSums(groups.level$mat) == 0] = 0.5
#'
#' val.mat = fluxing(groups.level$mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' bioms.pref = TRUE,
#' ef.level = "pred")
#'
#' stability.value(val.mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' growth.rates,
#' ef.level = "pred")
#'
#' @author Benoit Gauzens, \email{benoit.gauzens@gmail.com}
#'
#' @export
stability.value = function(val.mat,
biomasses,
losses,
efficiencies,
growth.rate,
bioms.prefs = TRUE,
bioms.losses = TRUE,
ef.level = "prey",
full.output = FALSE){
# mat
if (! is.numeric(val.mat)){
stop("'val.mat' must be numeric")
}
if (dim(val.mat)[1] != dim(val.mat)[2]){
stop("val.mat should be a square matrix")
}
# biomasses
if (! is.vector(biomasses)){
stop("biomasses should be a vector")
} else{
if (length(biomasses) != dim(val.mat)[1]){
stop("length of biomasses vector should equal to dimensions of val.mat")
}
}
if (! is.numeric(biomasses)){
stop("'biomasses' must be numeric")
} else if (any(biomasses <= 0)){
stop("'biomasses' must be all >0")
}
# losses
if (! is.numeric(losses)){
stop("'losses' should be numeric")
} else if (any(losses < 0)){
stop("'losses' contain negative value(s)")
}
# efficiencies
if (! is.numeric(efficiencies)){
stop("'efficiencies' must be numeric")
if (!(any(efficiencies < 0) || any(efficiencies > 1))) {
stop("'efficiencies' must be all in interval [0,1]")
}
if (is.vector(efficiencies)){
if (length(efficiencies) != dim(val.mat)[1]){
stop("'efficiencies' vector length sould be equal to number of species (dimension of val.mat)")
}
} else if (dim(efficiencies != dim(val.mat))){
stop("'efficiencies' matrix dimension different from 'val.mat'")
}
}
if (!(ef.level %in% c('prey', 'pred', 'link.specific'))){
stop("ef.level should be set to 'pred', 'prey' or 'link.specific'")
}
if (ef.level == 'prey' && is.matrix(efficiencies)){ # if user did not change the ef.level = "prey" optional argument but provide a matrix for efficiencies:
warning("'ef.level' is set to 'prey' and expect a vector of efficiencies but get a matrix instead.\n ef.level was then set to 'link.specific'")
ef.level = 'link.specific'
}
# growth.rate
if (! is.numeric(growth.rate)){
stop("growth.rate should be numeric")
} else if (any(growth.rate[colSums(val.mat) == 0] < 0)){
stop("all growth.rate values are expected to be positive (or 0)")
}
if (!is.vector(growth.rate)){
stop("growth.rate should be a vector")
} else if (length(growth.rate) != dim(val.mat)[1]){
stop("length of growth.rate vector should be equal to the number of species")
}
if (sum(!is.na(growth.rate[colSums(val.mat)!=0])) > 0){
warning('some growth rates are defined for non basal species. Values will be ignored')
}
if (sum(is.na(growth.rate[which(colSums(val.mat)==0)])) > 0){
stop('growth rates of some basal species are not defined')
}
### define loss vector as the sum of species losses:
if (! is.vector(losses)){
losses = rowSums(losses)
}
if (bioms.losses == T){
losses = losses*biomasses
}
nb_s = dim(val.mat)[1]
nb_b = sum(colSums(val.mat) == 0)
jacob = matrix(0, nb_s, nb_s)
jacob = create.jacob(val.mat, biomasses, losses, efficiencies, growth.rate, bioms.prefs = TRUE, bioms.losses = TRUE, ef.level = "prey")
# return(eigen(jacob))
if (full.output){
# should also return the Jacobian
return(eigen(jacob))
} else{
return(max(Re(eigen(jacob)$values)))
}
}
#' making network stability
#'
#' Find the smallest scalar multiplying a variable from losses insuring system stability
#'
#' @param val.mat A matrix describing fluxes between species (usually a result of \code{\link[fluxweb]{fluxing}} function).
#' @param biomasses A vector of species biomasses.
#' @param losses A vector or an array of species energy losses (excluding predation).
#' @param efficiencies A vector or an array of conversion efficiencies of species in the adjacency matrix. These values describe the proportion of consumed energy that is converted to biomass of the consumer.
#' @param growth.rate A vector defining growth rate of basal species.
#' @param losses.scale Defines a Column from \code{losses} on which scalar multiplication will be tested. (default \code{NULL} if the value is independent of losses).
#' @param bioms.prefs Logical, if \code{TRUE} (default) preferences are scaled accordingly to species biomasses.
#' @param bioms.losses Logical, if \code{TRUE} (default) losses are scaled with biomass.
#' @param ef.level Set to \code{"prey"} if efficiencies are defined by prey, \code{"pred"} if they are a property of the predator.
#' @param interval Search interval for returned value.
#' @param ... Optional parameters for function \code{\link[stats]{uniroot}}
#'
#'
#' @return
#' A list from \code{\link[stats]{uniroot}} function.
#'
#' @details
#' The function assumes a monotonous increase of stability with multiplication by a scalar value. Solution is estimated from the \code{\link[stats]{uniroot}} function, and stability using the \code{\link[fluxweb]{fluxing}} function
#' Thus, accordingly to \code{\link[stats]{uniroot}} solving criteria, if stability values at the two extremum parts of the interval are of same sign, an error is raised.
#'
#' Behavior of the multiplicative term depends on the type of losses:
#' \itemize{
#' \item{\code{losses.scale = NULL} and \code{is.vector(losses)}:} multiplication will be applied to the \code{losses} vector.
#' \item{\code{losses.scale = NULL} and \code{is.matrix(losses)}:} multiplication will be independent of any columns from \code{losses}.
#' \item{\code{losses.scale = FALSE} :} value used for multiplication always independent of losses.
#' \item{other values:} should refer to an element of losses.
#' }
#'
#' @seealso \code{\link[stats]{uniroot}} for root estimate and \code{\link[fluxweb]{stability.value}} for assessing system stability.
#'
#'
#' @examples
#'
# # first compute species per unit biomass metabolic rates using the metabolic theory:
#' losses = 0.15 * groups.level$bodymasses^(-0.25)
#'
#' # growth rates of basal sppecies
#' growth.rates = rep(NA, dim(groups.level$mat)[1])
#' growth.rates[colSums(groups.level$mat) == 0] = 0.5
#'
#' val.mat = fluxing(groups.level$mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' bioms.pref = TRUE,
#' ef.level = "pred")
#' make.stability(val.mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' growth.rates,
#' ef.level = "pred")
#'
#' @export
#'
#'
make.stability = function(val.mat,
biomasses,
losses,
efficiencies,
growth.rate,
losses.scale = NULL,
bioms.prefs = TRUE,
bioms.losses = TRUE,
ef.level = "prey",
interval = c(1e-12,1),
...){
stability.value.wrapper = function(x, unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level){
tot.losses = x*col + unvariant
return(stability.value(val.mat, biomasses, tot.losses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level))
}
# mat
if (! is.numeric(val.mat)){
stop("'val.mat' must be numeric")
}
if (dim(val.mat)[1] != dim(val.mat)[2]){
stop("val.mat should be a square matrix")
}
nb_s = dim(val.mat)[1]
# biomasses
if (! is.vector(biomasses)){
stop("biomasses should be a vector")
} else{
if (length(biomasses) != dim(val.mat)[1]){
stop("length of biomasses vector should equal to dimensions of val.mat")
}
}
if (! is.numeric(biomasses)){
stop("'biomasses' must be numeric")
} else if (any(biomasses <= 0)){
stop("'biomasses' must be all >0")
}
# losses
if (! is.numeric(losses)){
stop("'losses' should be numeric")
} else if (any(losses < 0)){
stop("'losses' contain negative value(s)")
}
# efficiencies
if (! is.numeric(efficiencies)){
stop("'efficiencies' must be numeric")
if (!(any(efficiencies < 0) || any(efficiencies > 1))) {
stop("'efficiencies' must be all in interval [0,1]")
}
if (is.vector(efficiencies)){
if (length(efficiencies) != dim(val.mat)[1]){
stop("'efficiencies' vector length sould be equal to number of species (dimension of val.mat)")
}
} else if (dim(efficiencies != dim(val.mat))){
stop("'efficiencies' matrix dimension different from 'val.mat'")
}
}
#ef.level
if ((ef.level != "prey") && (ef.level != "pred")){
stop("ef.level should be either 'pred' or 'prey'")
}
# growth.rate
if (! is.numeric(growth.rate)){
stop("growth.rate should be numeric")
} else if (any(growth.rate[colSums(val.mat) == 0] < 0)){
stop("all growth.rate values are expected to be positive (or 0)")
}
if (!is.vector(growth.rate)){
stop("growth.rate should be a vector")
} else if (length(growth.rate) != dim(val.mat)[1]){
stop("length of growth.rate vector should be equal to the number of species")
}
if (sum(!is.na(growth.rate[colSums(val.mat)!=0])) > 0){
warning('some growth rates are defined for non basal species. Values will be ignored')
}
if (sum(is.na(growth.rate[which(colSums(val.mat)==0)])) > 0){
stop('growth rates of some basal species are not defined')
}
# compute losses
# here, operation is:
# final.losses = colSums(invariant + value*col), where value is the value to estimate
if (missing(losses.scale)){
if (is.vector(losses)){
unvariant = rep(0, nb_s)
col = losses
}
if (is.matrix(losses)){
#generate a warning as behaviour will be the same as losses.scale = false?
unvariant = colSums(losses)
col = rep[1, nb_s]
}
} else{
if (is.numeric(losses.scale) && !is.integer(losses.scale)){
stop("losses.scale should be a either character or non integer numeric")
}
if (is.vector(losses) && losses.scale == FALSE){
unvariant = losses
col = rep(1, nb_s)
}
if (is.matrix(losses) && losses.scale == FALSE){
unvariant = colSums(losses)
col = rep(1, nb_s)
}
if (is.matrix(losses) && is.integer(losses.scale)){
if ((losses.scale <1 ) || (losses.scale > dim(val.mat[1]))){
stop("losses.scale should be between 1 and species number")
}
col = losses[,losses.scale]
unvariant = colSums(losses[, -losses.scale])
}
if (is.matrix(losses) && is.character(losses.scale)){
if (!(losses.scale %in% names(losses))){
stop("losses.scale not a column name of losses")
}
col = losses[,names(losses.scale) == losses.scale]
unvariant = colSums(losses[, colnames(losses) != losses.scale])
}
}
# get stability values to at the interval boundaries
stab.upper = stability.value.wrapper(max(interval), unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level)
stab.lower = stability.value.wrapper(min(interval), unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level)
if (stab.lower < 0 && stab.upper <0){
stop("unable to assess minimal stability value: system always stable in specified interval.")
}
if (stab.lower > 0 && stab.upper > 0){
stop("unable to assess minimal stability value: system never stable in specified interval.")
}
min.stab.value = stats::uniroot(stability.value.wrapper, interval, ..., unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level,
lower = min(interval), upper = max(interval), f.lower = stab.lower, f.upper = stab.upper)
return(min.stab.value)
}
#' Internal function to compute Jacobian
#'
#'@keywords internal
#'
#'
#'
create.jacob = function(val.mat, biomasses, losses, efficiencies, growth.rate, bioms.prefs = TRUE, bioms.losses = TRUE, ef.level = "prey"){
nb_s = dim(val.mat)[1]
nb_b = sum(colSums(val.mat) == 0)
jacob = matrix(0, nb_s, nb_s)
basal.index = which(colSums(val.mat) == 0) # store the position of basal species in arguments
# double loop here... maybe possible to vectorise
for (i in 1:nb_s){
for (j in 1: nb_s){
if (i == j){
if (i %in% basal.index){ # then species i is basal
jacob[i,j] = growth.rate[i] - sum(val.mat[i,])/biomasses[i]
} else{
if (ef.level == "pred") {
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i] - 1) / biomasses[i] +
efficiencies[i] * sum(val.mat[,i])/biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
if (ef.level == "prey" && is.vector(efficiencies)){
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i] - 1) / biomasses[i] +
sum(val.mat[,i] * efficiencies) / biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
if (ef.level == "link.specific"){
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i,i] - 1) / biomasses[i] +
sum(val.mat[,i] * efficiencies[,i]) / biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
}
} else{
if (ef.level == "pred") {jacob[i,j] = (efficiencies[i]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
if (ef.level == "prey" && is.vector(efficiencies)) {jacob[i,j] = (efficiencies[j]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
# should change the following line to consider a flag for efficiencies = "link.level" instead of checking if it's a matri or not
if (ef.level == "link.specific") {jacob[i,j] = (efficiencies[j,i]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
}
}
}
return(jacob)
}
Try the fluxweb package in your browser
Any scripts or data that you put into this service are public.
fluxweb documentation built on May 2, 2019, 6:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' Estimates network stability
#'
#' Computes resilience of the system through Jacobian matrix eigenvalues.
#'
#'
#' @param val.mat A matrix describing fluxes between species (usually a result of \code{\link[fluxweb]{fluxing}} function).
#' @param biomasses A vector of species biomasses.
#' @param losses A vector or an array of species energy losses (excluding predation).
#' @param efficiencies A vector or an array of conversion efficiencies of species in the adjacency matrix. These values describe the proportion of consumed energy that is converted to biomass of the consumer.
#' @param growth.rate A vector defining growth rate of basal species.
#' @param bioms.prefs Logical, if \code{TRUE} (default) preferences are scaled according to species biomasses.
#' @param bioms.losses Logical, if \code{TRUE} (default) losses are scaled with biomass.
#' @param ef.level Set to \code{"prey"} if efficiencies are defined by prey, \code{"pred"} if they are a property of the predator.
#' @param full.output Logical, if \code{TRUE} function return supplementary informations.
#'
#' @return Maximum eigenvalue of the Jacobian matrix of a Lotka Voltera like system of equations. If full.output, Jacobian eigenvalues and eigenvectors are returned.
#'
#'
#'@details
#'
#'
#For more information about the set of the underlying system of differential equations and mathematical derivation, please read Gauzens et al. 2017 SI2. <https://www.biorxiv.org/content/early/2017/12/06/229450>
#'
#'
#'\itemize{
#'
#'\item{\code{losses}:} Express species energetic losses not related to consumption. Usually metabolic or death rates.
#'When an array is provided, losses associated to each species correspond to line sums.
#'
#'\item{\code{efficiencies}:} Determines how efficient species are to convert energy (see \code{ef.level} for more details).
#'Providing an array will assume values depending on both prey and predator identity.
#'
#'\item{\code{growth.rate}:} Growth rates of basal species defined. Length of the vector should be equal to the number of species.
#'expects positive numeric values for index corresponding to basal species, NA otherwise
#'
#'\item{\code{bioms.pref}:} If \code{TRUE}, preferences \eqn{w_{ij}} of predator j on prey i are scaled according to species biomass using the following formula:
#'\deqn{
#'w_{i,j} = \frac{mat[i,j] * biomasses[i]}{\sum_k mat[i,k]* biomasses[k]}
#'}
#'\item{\code{bioms.losses}:} If \code{TRUE}, function will assume that losses are defined per biomass unit.
#'Thus, total losses will be thereafter multiplied by biomass values for each species.
#'
#'\item{\code{ef.level}:} If \code{"prey"} (resp \code{"pred"}), the total amount of energy that can be metabolized from a trophic link
#'will be determined by prey (resp pred) identity. \code{"link.specific"} assumes that efficiencies are defined for each trophic interaction
#'and implies \code{efficiencies} parameter to be a matrix
#'
#'\item{\code{full.output}:} If \code{TRUE}, function result is a list of eigenvalues and eigenvectors of the Jacobian matrix.
#'}
#'
#' @examples
# # first compute species per unit biomass metabolic rates using the metabolic theory:
#' losses = 0.15 * groups.level$bodymasses^(-0.25)
#'
#' # growth rates of basal sppecies
#' growth.rates = rep(NA, dim(groups.level$mat)[1])
#' growth.rates[colSums(groups.level$mat) == 0] = 0.5
#'
#' val.mat = fluxing(groups.level$mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' bioms.pref = TRUE,
#' ef.level = "pred")
#'
#' stability.value(val.mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' growth.rates,
#' ef.level = "pred")
#'
#' @author Benoit Gauzens, \email{benoit.gauzens@gmail.com}
#'
#' @export
stability.value = function(val.mat,
biomasses,
losses,
efficiencies,
growth.rate,
bioms.prefs = TRUE,
bioms.losses = TRUE,
ef.level = "prey",
full.output = FALSE){
# mat
if (! is.numeric(val.mat)){
stop("'val.mat' must be numeric")
}
if (dim(val.mat)[1] != dim(val.mat)[2]){
stop("val.mat should be a square matrix")
}
# biomasses
if (! is.vector(biomasses)){
stop("biomasses should be a vector")
} else{
if (length(biomasses) != dim(val.mat)[1]){
stop("length of biomasses vector should equal to dimensions of val.mat")
}
}
if (! is.numeric(biomasses)){
stop("'biomasses' must be numeric")
} else if (any(biomasses <= 0)){
stop("'biomasses' must be all >0")
}
# losses
if (! is.numeric(losses)){
stop("'losses' should be numeric")
} else if (any(losses < 0)){
stop("'losses' contain negative value(s)")
}
# efficiencies
if (! is.numeric(efficiencies)){
stop("'efficiencies' must be numeric")
if (!(any(efficiencies < 0) || any(efficiencies > 1))) {
stop("'efficiencies' must be all in interval [0,1]")
}
if (is.vector(efficiencies)){
if (length(efficiencies) != dim(val.mat)[1]){
stop("'efficiencies' vector length sould be equal to number of species (dimension of val.mat)")
}
} else if (dim(efficiencies != dim(val.mat))){
stop("'efficiencies' matrix dimension different from 'val.mat'")
}
}
if (!(ef.level %in% c('prey', 'pred', 'link.specific'))){
stop("ef.level should be set to 'pred', 'prey' or 'link.specific'")
}
if (ef.level == 'prey' && is.matrix(efficiencies)){ # if user did not change the ef.level = "prey" optional argument but provide a matrix for efficiencies:
warning("'ef.level' is set to 'prey' and expect a vector of efficiencies but get a matrix instead.\n ef.level was then set to 'link.specific'")
ef.level = 'link.specific'
}
# growth.rate
if (! is.numeric(growth.rate)){
stop("growth.rate should be numeric")
} else if (any(growth.rate[colSums(val.mat) == 0] < 0)){
stop("all growth.rate values are expected to be positive (or 0)")
}
if (!is.vector(growth.rate)){
stop("growth.rate should be a vector")
} else if (length(growth.rate) != dim(val.mat)[1]){
stop("length of growth.rate vector should be equal to the number of species")
}
if (sum(!is.na(growth.rate[colSums(val.mat)!=0])) > 0){
warning('some growth rates are defined for non basal species. Values will be ignored')
}
if (sum(is.na(growth.rate[which(colSums(val.mat)==0)])) > 0){
stop('growth rates of some basal species are not defined')
}
### define loss vector as the sum of species losses:
if (! is.vector(losses)){
losses = rowSums(losses)
}
if (bioms.losses == T){
losses = losses*biomasses
}
nb_s = dim(val.mat)[1]
nb_b = sum(colSums(val.mat) == 0)
jacob = matrix(0, nb_s, nb_s)
jacob = create.jacob(val.mat, biomasses, losses, efficiencies, growth.rate, bioms.prefs = TRUE, bioms.losses = TRUE, ef.level = "prey")
# return(eigen(jacob))
if (full.output){
# should also return the Jacobian
return(eigen(jacob))
} else{
return(max(Re(eigen(jacob)$values)))
}
}
#' making network stability
#'
#' Find the smallest scalar multiplying a variable from losses insuring system stability
#'
#' @param val.mat A matrix describing fluxes between species (usually a result of \code{\link[fluxweb]{fluxing}} function).
#' @param biomasses A vector of species biomasses.
#' @param losses A vector or an array of species energy losses (excluding predation).
#' @param efficiencies A vector or an array of conversion efficiencies of species in the adjacency matrix. These values describe the proportion of consumed energy that is converted to biomass of the consumer.
#' @param growth.rate A vector defining growth rate of basal species.
#' @param losses.scale Defines a Column from \code{losses} on which scalar multiplication will be tested. (default \code{NULL} if the value is independent of losses).
#' @param bioms.prefs Logical, if \code{TRUE} (default) preferences are scaled accordingly to species biomasses.
#' @param bioms.losses Logical, if \code{TRUE} (default) losses are scaled with biomass.
#' @param ef.level Set to \code{"prey"} if efficiencies are defined by prey, \code{"pred"} if they are a property of the predator.
#' @param interval Search interval for returned value.
#' @param ... Optional parameters for function \code{\link[stats]{uniroot}}
#'
#'
#' @return
#' A list from \code{\link[stats]{uniroot}} function.
#'
#' @details
#' The function assumes a monotonous increase of stability with multiplication by a scalar value. Solution is estimated from the \code{\link[stats]{uniroot}} function, and stability using the \code{\link[fluxweb]{fluxing}} function
#' Thus, accordingly to \code{\link[stats]{uniroot}} solving criteria, if stability values at the two extremum parts of the interval are of same sign, an error is raised.
#'
#' Behavior of the multiplicative term depends on the type of losses:
#' \itemize{
#' \item{\code{losses.scale = NULL} and \code{is.vector(losses)}:} multiplication will be applied to the \code{losses} vector.
#' \item{\code{losses.scale = NULL} and \code{is.matrix(losses)}:} multiplication will be independent of any columns from \code{losses}.
#' \item{\code{losses.scale = FALSE} :} value used for multiplication always independent of losses.
#' \item{other values:} should refer to an element of losses.
#' }
#'
#' @seealso \code{\link[stats]{uniroot}} for root estimate and \code{\link[fluxweb]{stability.value}} for assessing system stability.
#'
#'
#' @examples
#'
# # first compute species per unit biomass metabolic rates using the metabolic theory:
#' losses = 0.15 * groups.level$bodymasses^(-0.25)
#'
#' # growth rates of basal sppecies
#' growth.rates = rep(NA, dim(groups.level$mat)[1])
#' growth.rates[colSums(groups.level$mat) == 0] = 0.5
#'
#' val.mat = fluxing(groups.level$mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' bioms.pref = TRUE,
#' ef.level = "pred")
#' make.stability(val.mat,
#' groups.level$biomasses,
#' losses,
#' groups.level$efficiencies,
#' growth.rates,
#' ef.level = "pred")
#'
#' @export
#'
#'
make.stability = function(val.mat,
biomasses,
losses,
efficiencies,
growth.rate,
losses.scale = NULL,
bioms.prefs = TRUE,
bioms.losses = TRUE,
ef.level = "prey",
interval = c(1e-12,1),
...){
stability.value.wrapper = function(x, unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level){
tot.losses = x*col + unvariant
return(stability.value(val.mat, biomasses, tot.losses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level))
}
# mat
if (! is.numeric(val.mat)){
stop("'val.mat' must be numeric")
}
if (dim(val.mat)[1] != dim(val.mat)[2]){
stop("val.mat should be a square matrix")
}
nb_s = dim(val.mat)[1]
# biomasses
if (! is.vector(biomasses)){
stop("biomasses should be a vector")
} else{
if (length(biomasses) != dim(val.mat)[1]){
stop("length of biomasses vector should equal to dimensions of val.mat")
}
}
if (! is.numeric(biomasses)){
stop("'biomasses' must be numeric")
} else if (any(biomasses <= 0)){
stop("'biomasses' must be all >0")
}
# losses
if (! is.numeric(losses)){
stop("'losses' should be numeric")
} else if (any(losses < 0)){
stop("'losses' contain negative value(s)")
}
# efficiencies
if (! is.numeric(efficiencies)){
stop("'efficiencies' must be numeric")
if (!(any(efficiencies < 0) || any(efficiencies > 1))) {
stop("'efficiencies' must be all in interval [0,1]")
}
if (is.vector(efficiencies)){
if (length(efficiencies) != dim(val.mat)[1]){
stop("'efficiencies' vector length sould be equal to number of species (dimension of val.mat)")
}
} else if (dim(efficiencies != dim(val.mat))){
stop("'efficiencies' matrix dimension different from 'val.mat'")
}
}
#ef.level
if ((ef.level != "prey") && (ef.level != "pred")){
stop("ef.level should be either 'pred' or 'prey'")
}
# growth.rate
if (! is.numeric(growth.rate)){
stop("growth.rate should be numeric")
} else if (any(growth.rate[colSums(val.mat) == 0] < 0)){
stop("all growth.rate values are expected to be positive (or 0)")
}
if (!is.vector(growth.rate)){
stop("growth.rate should be a vector")
} else if (length(growth.rate) != dim(val.mat)[1]){
stop("length of growth.rate vector should be equal to the number of species")
}
if (sum(!is.na(growth.rate[colSums(val.mat)!=0])) > 0){
warning('some growth rates are defined for non basal species. Values will be ignored')
}
if (sum(is.na(growth.rate[which(colSums(val.mat)==0)])) > 0){
stop('growth rates of some basal species are not defined')
}
# compute losses
# here, operation is:
# final.losses = colSums(invariant + value*col), where value is the value to estimate
if (missing(losses.scale)){
if (is.vector(losses)){
unvariant = rep(0, nb_s)
col = losses
}
if (is.matrix(losses)){
#generate a warning as behaviour will be the same as losses.scale = false?
unvariant = colSums(losses)
col = rep[1, nb_s]
}
} else{
if (is.numeric(losses.scale) && !is.integer(losses.scale)){
stop("losses.scale should be a either character or non integer numeric")
}
if (is.vector(losses) && losses.scale == FALSE){
unvariant = losses
col = rep(1, nb_s)
}
if (is.matrix(losses) && losses.scale == FALSE){
unvariant = colSums(losses)
col = rep(1, nb_s)
}
if (is.matrix(losses) && is.integer(losses.scale)){
if ((losses.scale <1 ) || (losses.scale > dim(val.mat[1]))){
stop("losses.scale should be between 1 and species number")
}
col = losses[,losses.scale]
unvariant = colSums(losses[, -losses.scale])
}
if (is.matrix(losses) && is.character(losses.scale)){
if (!(losses.scale %in% names(losses))){
stop("losses.scale not a column name of losses")
}
col = losses[,names(losses.scale) == losses.scale]
unvariant = colSums(losses[, colnames(losses) != losses.scale])
}
}
# get stability values to at the interval boundaries
stab.upper = stability.value.wrapper(max(interval), unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level)
stab.lower = stability.value.wrapper(min(interval), unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level)
if (stab.lower < 0 && stab.upper <0){
stop("unable to assess minimal stability value: system always stable in specified interval.")
}
if (stab.lower > 0 && stab.upper > 0){
stop("unable to assess minimal stability value: system never stable in specified interval.")
}
min.stab.value = stats::uniroot(stability.value.wrapper, interval, ..., unvariant, col, val.mat, biomasses, efficiencies, growth.rate, bioms.prefs, bioms.losses, ef.level,
lower = min(interval), upper = max(interval), f.lower = stab.lower, f.upper = stab.upper)
return(min.stab.value)
}
#' Internal function to compute Jacobian
#'
#'@keywords internal
#'
#'
#'
create.jacob = function(val.mat, biomasses, losses, efficiencies, growth.rate, bioms.prefs = TRUE, bioms.losses = TRUE, ef.level = "prey"){
nb_s = dim(val.mat)[1]
nb_b = sum(colSums(val.mat) == 0)
jacob = matrix(0, nb_s, nb_s)
basal.index = which(colSums(val.mat) == 0) # store the position of basal species in arguments
# double loop here... maybe possible to vectorise
for (i in 1:nb_s){
for (j in 1: nb_s){
if (i == j){
if (i %in% basal.index){ # then species i is basal
jacob[i,j] = growth.rate[i] - sum(val.mat[i,])/biomasses[i]
} else{
if (ef.level == "pred") {
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i] - 1) / biomasses[i] +
efficiencies[i] * sum(val.mat[,i])/biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
if (ef.level == "prey" && is.vector(efficiencies)){
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i] - 1) / biomasses[i] +
sum(val.mat[,i] * efficiencies) / biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
if (ef.level == "link.specific"){
jacob[i,j] = -losses[i] + val.mat[i,j] * (efficiencies[i,i] - 1) / biomasses[i] +
sum(val.mat[,i] * efficiencies[,i]) / biomasses[i] - sum(val.mat[i,]) / biomasses[i]
}
}
} else{
if (ef.level == "pred") {jacob[i,j] = (efficiencies[i]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
if (ef.level == "prey" && is.vector(efficiencies)) {jacob[i,j] = (efficiencies[j]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
# should change the following line to consider a flag for efficiencies = "link.level" instead of checking if it's a matri or not
if (ef.level == "link.specific") {jacob[i,j] = (efficiencies[j,i]*val.mat[j,i] - val.mat[i,j])/biomasses[j]}
}
}
}
return(jacob)
}
Try the fluxweb package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.