Nothing
################################################################################
#' Transfer function analysis
#'
#' @description
#' Transfer function analysis of the dominant eigenvalue of a population matrix
#' projection model using a specified perturbation structure.
#'
#' @param A a square, irreducible, nonnegative numeric matrix of any dimension
#' @param d,e numeric vectors that determine the perturbation structure
#' (see details).
#' @param lambdarange a numeric vector giving the range of lambda values
#' (asymptotic growth rates) to be achieved (see details).
#' @param prange a numeric vector giving the range of perturbation magnitude
#' (see details)
#' @param digits specifies which values of lambda should be excluded from
#' analysis to avoid a computationally singular system (see details).
#'
#' @details
#' \code{tfa_lambda} calculates the transfer function of the dominant eigenvalue
#' of a matrix (\code{A}), given a perturbation structure (specified using
#' \code{d} and \code{e}), and either a range of target values for asymptotic
#' population growth (\code{lambdavalues}) or a range of desired perturbation
#' magnitude (\code{prange}). Currently \code{tfa_lambda} can only work with rank-
#' one, single-parameter perturbations (see Hodgson & Townley 2004).
#'
#' The perturbation structure is determined by \code{d\%*\%t(e)}. Therefore,
#' the rows to be perturbed are determined by \code{d} and the columns to be
#' perturbed are determined by \code{e}. The specific values in d and e will
#' determine the relative perturbation magnitude. So for example, if only entry
#' [3,2] of a 3 by 3 matrix is to be perturbed, then \code{d = c(0,0,1)} and
#' \code{e = c(0,1,0)}. If entries [3,2] and [3,3] are to be perturbed with the
#' magnitude of perturbation to [3,2] half that of [3,3] then \code{d = c(0,0,1)}
#' and \code{e = c(0,0.5,1)}. \code{d} and \code{e} may also be expressed as
#' numeric one-column matrices, e.g. \code{d = matrix(c(0,0,1), ncol=1)},
#' \code{e = matrix(c(0,0.5,1), ncol=1)}. See Hodgson et al. (2006) for more
#' information on perturbation structures.
#'
#' The perturbation magnitude is determined by \code{prange}, a numeric vector
#' that gives the continuous range of perterbation magnitude to evaluate over.
#' This is usually a sequence (e.g. \code{prange=seq(-0.1,0.1,0.001)}), but
#' single transfer functions can be calculated using a single perturbation
#' magnitude (e.g. \code{prange=-0.1}). Because of the nature of the equation,
#' \code{prange} is used to find a range of lambda from which the perturbation
#' magnitudes are back-calculated, so the output perturbation magnitude
#' \code{p} will match \code{prange} in length and range but not in numerical
#' value (see value). Alternatively, a vector \code{lambdarange} can be
#' specified, representing a range of target lambda values from which the
#' corresponding perturbation values will be calculated. Only one of either
#' \code{prange} or \code{lambdarange} may be specified.
#'
#' \code{tfa_lambda} uses the resolvent matrix in its calculation, which cannot be
#' computed if any lambda are equal to the dominant eigenvalue of
#' \code{A}. \code{digits} specifies the values of lambda that should be
#' excluded in order to avoid a computationally singular system. Any values
#' equal to the dominant eigenvalue of \code{A} rounded to an accuracy of
#' \code{digits} are excluded. \code{digits} should only need to be changed
#' when the system is found to be computationally singular, in which case
#' increasing \code{digits} should help to solve the problem.
#'
#' \code{tfa_lambda} will not work for reducible matrices.
#'
#' There is an S3 plotting method available (see \code{\link{plot.tfa}}
#' and examples below)
#'
#' @return
#' A list containing numerical vectors:
#' \describe{
#' \item{p}{the perturbation magnitude(s).}
#' \item{lambda}{the dominant eigenvalue(s) of the perturbed matrix(matrices)}.
#' }
#' (Note that \code{p} will not be exactly the same as \code{prange} when
#' \code{prange} is specified, as the code calculates p for a given lambda
#' rather than the other way around, with \code{prange} only used to determine
#' max, min and number of lambda values to evaluate.)
#'
#' @references
#' \itemize{
#' \item Townley & Hodgson (2004) J. Appl. Ecol., 41, 1155-1161.
#' \item Hodgson et al. (2006) J. Theor. Biol., 70, 214-224.
#' }
#'
#' @family TransferFunctionAnalyses
#' @family PerturbationAnalyses
#'
#' @seealso
#' S3 plotting method: \code{\link{plot.tfa}}
#'
#' @examples
#' # Create a 3x3 matrix
#' ( A <- matrix(c(0,1,2,0.5,0.1,0,0,0.6,0.6), byrow=TRUE, ncol=3) )
#'
#' # Calculate the transfer function of A[3,2] given a range of lambda
#' evals <- eigen(A)$values
#' lmax <- which.max(Re(evals))
#' lambda <- Re(evals[lmax])
#' lambdarange <- seq(lambda-0.1, lambda+0.1, 0.01)
#' ( transfer <- tfa_lambda(A, d=c(0,0,1), e=c(0,1,0), lambdarange=lambdarange) )
#'
#' # Plot the transfer function using the S3 method
#' plot(transfer)
#'
#' # Calculate the transfer function of perturbation to A[3,2] and A[3,3]
#' # with perturbation to A[3,2] half that of A[3,3], given a range of
#' # perturbation values
#' p<-seq(-0.6,0.4,0.01)
#' ( transfer2 <- tfa_lambda(A, d=c(0,0,1), e=c(0,0.5,1), prange=p) )
#'
#' # Plot p and lambda without using the S3 method
#' plot(transfer$lambda~transfer$p, type="l", xlab="p", ylab=expression(lambda))
#'
#' @concept transfer function
#' @concept systems control
#' @concept nonlinear
#' @concept perturbation
#' @concept population viability
#' @concept PVA
#' @concept ecology
#' @concept demography
#' @concept population
#'
#' @export tfa_lambda
#'
tfa_lambda <-
function(A,d,e,prange=NULL,lambdarange=NULL,digits=1e-10){
if(any(length(dim(A))!=2,dim(A)[1]!=dim(A)[2])) stop("A must be a square matrix")
order<-dim(A)[1]
if(!isIrreducible(A)) stop("Matrix A is reducible")
if(!isPrimitive(A)) warning("Matrix A is imprimitive")
eigvals<-eigen(A)$values
lmax<-which.max(Re(eigvals))
lambda<-Re(eigvals[lmax])
if(any(is.null(d),is.null(e))) stop("please specify a perturbation structure using d and e")
d<-as.matrix(d)
e<-as.matrix(e)
if(is.null(lambdarange)&is.null(prange)) stop("Please specify one of either lambdarange or prange")
if(!is.null(lambdarange)&!is.null(prange)) stop("Please specify only one of either lambdarange or prange")
if(!is.null(lambdarange)&is.null(prange)){
lambdarange<-lambdarange[round(lambdarange,-log10(digits))!=round(lambda,-log10(digits))]
pvals<-1/.tf(A,lambdarange,d,e)
}
if(is.null(lambdarange)&!is.null(prange)){
mineigvals<-eigen(A+(min(prange)*d%*%t(e)))$values
maxeigvals<-eigen(A+(max(prange)*d%*%t(e)))$values
minlmax<-which.max(Re(mineigvals))
minlambda<-Re(mineigvals[minlmax])
maxlmax<-which.max(Re(maxeigvals))
maxlambda<-Re(maxeigvals[maxlmax])
lambdarange<-seq(minlambda,maxlambda,(maxlambda-minlambda)/(length(prange)-1))
lambdarange<-lambdarange[round(lambdarange,-log10(digits))!=round(lambda,-log10(digits))]
pvals<-1/.tf(A,lambdarange,d,e)
}
final<-list(p=pvals,lambda=lambdarange)
class(final)<-c("tfa","list")
return(final)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.