R/qage.R

In evolucionario/cladeage: CladeAge: Estimating Clade Age Using The Fossil Record

##' @title Quantile Function of Clade
##' @description Quantile function that computes the estimated age corresponding to a
##' particular probability for the upper bound of a distribution of ages.

##' @param p vector of probabilities.
##' @param ages vector of fossil ages.
##' @param method a character string specifying the method: \code{"StraussSadler"}
##' method (default),\code{"Solow"}, \code{"Beta"},\code{"NorrisPenGap}, \code{"NorrisGhostLin}, or \code{""RobertsSolow"}.
##' @param k the number of fossil ages to use in the “RobertsSolow” method, which only used the k oldest ages. 

##' @details
##' The "StraussSadler" (Strauss & Sadler 1989) method (default) assumes that
##' fossil ages are uniformly distributed, and a warning is returned if a Kolmogorov-Smirnov test rejects the uniformity hypothesis. The “Beta” method is a different parameterization of the Strauss & Sadler's method (Wang et al. 2009) that uses the qbeta function, since the ratio between the observed maximum fossil age and the clade age for the Strauss & Sadler model is distributed according to a Beta distribution with parameters N and 1 (Wang et al. 2009); this should give the same result as "StraussSadler". The "Solow" method (Solow 2003) does not assume uniformity but is based on the two oldest ages only. The “RobertsSolow” method does not assume uniformity and is based on the fact that the joint distribution of the k oldest fossil ages can be modeled with the same Weibull distribution (Roberts & Solow 2003). "NorrisPenGap" and "NorrisGhostLin" implement the method of Norris et al.(2015) based on the two oldest fossils only and the log-logistic distribution. "NorrisPenGap" is used when the precise phylogenetic placement of fossils is not known, whereas "NorrisGhostLin" is used when one fossil from each daughter lineage is used.

##' @return A numeric value (or vector of numeric values, if multiple p values
##' are provided) representing the age estimate of the clade origin given the
##' method a p value provided

##' @importFrom stats qbeta

##' @examples
##'   qage(p=c(0.1, 0.5, 0.9), ages=c(54, 30, 25, 14, 5));
##'   qage(p=c(0.1, 0.5, 0.9), ages=c(54, 30, 25, 14, 5), method="Beta");
##'   qage(p=c(0.1, 0.5, 0.9), ages=c(54, 30, 25, 14, 5), method="Solow");

##' @importFrom Rdpack reprompt
##' @references
##' \insertRef{Claramunt2015}{cladeage}
##' 
##' \insertRef{Gingerich1998}{cladeage}
##' 
##' \insertRef{Norris2015}{cladeage}
##' 
##' \insertRef{Solow2003}{cladeage}
##' 
##' \insertRef{Strauss1989}{cladeage}
##' 
##' \insertRef{Wang2007}{cladeage}
##' 
##' \insertRef{Wang2009}{cladeage}
##' 
##' \insertRef{Wang2010}{cladeage}
##' @export
qage <- function(p=0.5, ages,  method="StraussSadler", k=min(length(ages),5)) {
 	
 	if(length(ages) < 2) stop("More than one fossil age is needed for inference")
 	 	
 	ages <- sort(ages, decreasing=TRUE)

	if(method=="StraussSadler") {

		n <- length(ages)
		
		range <- max(ages)-min(ages)
		
		AGE <- range*(1-p)^-(1/n) + min(ages)
		}
		
	if(method=="Beta") {
		
		n <- length(ages)

		range <- max(ages) - min(ages)

		X <- range/qbeta(p=p, shape1=n, shape2=1, lower.tail=FALSE)

		AGE <- X + min(ages)
		
		 }
		
	if(method=="Solow") {
				
		AGE <- ages[1] + (ages[1]-ages[2])*p/(1-p) 

		}

	if(method=="NorrisPenGap") {	

		PenG <- ages[1] - ages[2]

		UltG <- exp(qlogis(p=p, location=log(PenG)))
		
		AGE <- UltG + ages[1] }
	
	if(method=="NorrisGhostLin") {
		
		GLin <- ages[1] - ages[2] 

		UltG <- exp(qlogis(p=p, location=log(GLin/2)))
		
		AGE <- UltG + ages[1]}
	
	if(method=="RobertsSolow") {
		
	# Select the k oldest observations (5 by default)
	
		ages <- ages[1:k]

	# Estimate the shape paremter of the joint Weibull distribution
			
		t <- numeric()
	
		for(i in 1:(k-2)) { t[i] <- (ages[1] - ages[k]) / (ages[1] - ages[i+1]) }

		v <- 1/(k-1) * sum(log(t))

	# Calculate Su (Robert & Solow 2003, = c(alpha) Solow 2005 eq. 18)
		Su <- (-log(1-p)/k)^-v

	# But Su must be >= 1 (otherwise the estimate may be younger than older fossil) so:
		Su <- pmax(Su, 1) # use pmax instead of max to allow for vectors of p and Su
	
	# Calculate the quantile (Solow 2005 eq. 17) (Something seems wrong)	
	#AGE <- ( ages[1] - Su*ages[k]) / ( 1 - Su)
	
	# Calculate the quantile (Roberts & Solow 2003)
	
		AGE <- ages[1] + ( ages[1] - ages[k]) / ( Su - 1 )

	}
	
	return(AGE)
}