R/tfa_lambda.R

In popdemo: Demographic Modelling Using Projection Matrices

################################################################################
#' 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)
}