Nothing
#' Checks that event times are correct
#'
#' @description Checks \code{event_times} and \code{event_times2} are of the
#' appropriate class and have expected characteristics for correct calculation
#' of NLTT in \code{\link{nltt_diff_exact_extinct}}.
#' @author Pedro Neves and Richèl Bilderbeek and Thijs Janzen
#'
#' @param event_times event times of the first phylogeny
#' @param event_times2 event times of the second phylogeny
#' @param time_unit the time unit of the branching times
#' \itemize{
#' \item{"ago: "}{the branching times are postive,
#' as these are in time units ago}
#' \item{"since: "}{the branching times are negative,
#' as these are in time units since present}
#' }
#'
#' @return Nothing. Throws error with helpful error message if
#' \code{event_times} and \code{event_times2} are not correct.
#' @export
check_input_event_times <- function(event_times, event_times2, time_unit) {
# This function doesn't handle phylo objects
if (!is.numeric(event_times) || !is.numeric(event_times2)) {
stop("event times should be a numeric vector of event times, not a
phylogeny. Use nltt_diff_exact_norm_brts for phylo objects instead.")
}
if (time_unit == "since") {
if (!all(event_times <= 0.0) || !all(event_times2 <= 0.0)) {
stop("event times must be negative, ",
"for example -3 time units since the present")
}
} else if (!all(event_times >= 0.0) || !all(event_times2 >= 0.0)) {
stop("event times must be positive, ",
"for example 3 time units ago")
}
}
#' Calculates the exact difference between the nLTT
#' curves of the event times. This includes extinction events.
#' @description Takes branching times such as (for example) as returned by the
#' DDD package.
#' @author Pedro Neves and Richèl Bilderbeek and Thijs Janzen
#' @inheritParams default_params_doc
#' @param event_times event times of the first phylogeny
#' @param species_number the number of species at each event time of the first
#' phylogeny
#' @param event_times2 event times of the second phylogeny
#' @param species_number2 the number of species at each event time of the second
#' phylogeny
#' @param distance_method how the difference between the two nLTTs is summed
#' \itemize{
#' \item{"abs: "}{the absolute distance between the two nLTTs is summed}
#' \item{"squ: "}{the squared distance between the two nLTTs is summed}
#' }
#' @param normalize should the output be normalized? Default is TRUE.
#'
#' @examples
#'
#' # Generate data
#' n <- 10
#' b_times_n <- (seq(1, n) / n)
#' lineages_n <- b_times_n
#' b_times2_n <- b_times_n * b_times_n
#' lineages2_n <- b_times2_n
#'
#' # Calculate nLTT
#' out <- nLTT::nltt_diff_exact_extinct(
#' event_times = b_times_n,
#' species_number = lineages_n,
#' event_times2 = b_times2_n,
#' species_number2 = lineages2_n,
#' time_unit = "ago",
#' distance_method = "abs"
#' )
#'
#' @export
nltt_diff_exact_extinct <- function(
event_times,
species_number,
event_times2,
species_number2,
distance_method = "abs",
time_unit = "since",
normalize = TRUE) {
nLTT::check_time_unit(time_unit)
if (!is.numeric(event_times)) {
stop(
"event times should be a numeric vector of event times, not a
phylogeny. Use nltt_diff_exact_norm_brts for phylo objects instead."
)
}
testit::assert(is.numeric(event_times2))
nLTT::check_input_event_times(
event_times = event_times,
event_times2 = event_times2,
time_unit = time_unit
)
# Each branching time must have a number of lineages to accompany it
testit::assert(length(event_times) == length(species_number))
testit::assert(length(event_times2) == length(species_number2))
if (time_unit == "ago") {
# Conformize te b_times for the classic calculation
event_times <- c(-1.0 * rev(sort(event_times)), 0)
event_times2 <- c(-1.0 * rev(sort(event_times2)), 0)
species_number <- c(species_number, utils::tail(species_number, n = 1))
species_number2 <- c(species_number2, utils::tail(species_number2, n = 1))
}
# Each branching time must have a number of species to accompany it
testit::assert(length(event_times) == length(species_number))
testit::assert(length(event_times2) == length(species_number2))
# We calculate with 'since the present' as a time unit
testit::assert(all(event_times <= 0.0))
testit::assert(all(event_times2 <= 0.0))
testit::assert(all(species_number >= 0.0))
testit::assert(all(species_number2 >= 0.0))
if (normalize) {
#normalize event times
event_times_n <- 1.0 - event_times / min(event_times)
#normalize lineages
species_number_n <- species_number / max(species_number)
# Normalizations must have worked
testit::assert(all(event_times_n >= 0.0 & event_times_n <= 1.0))
testit::assert(all(species_number_n >= 0.0 & species_number_n <= 1.0))
#normalize branching times
event_times2_n <- 1 - event_times2 / min(event_times2)
#normalize lineages
species_number2_n <- species_number2 / max(species_number2)
nltt_diff_exact_calc_extinct(
event_times = event_times_n,
species_number = species_number_n,
event_times2 = event_times2_n,
species_number2 = species_number2_n,
distance_method = distance_method
)
} else {
nltt_diff_exact_calc_extinct(
event_times = event_times,
species_number = species_number,
event_times2 = event_times2,
species_number2 = species_number2,
distance_method = distance_method
)
}
}
#' Calculates the exact difference between the nLTT
#' curves of the event times. This includes extinction events.
#' @author Thijs Janzen and Richèl Bilderbeek and Pedro Neves
#' @param event_times event times of the first phylogeny
#' @param species_number the number of species at each event time of the first
#' phylogeny
#' @param event_times2 event times of the second phylogeny
#' @param species_number2 the number of species at each event time of the second
#' phylogeny
#' @param distance_method (string) absolute, or squared distance?
#'
#' @examples
#'
#' # Generate data
#' n <- 10
#' b_times_n <- (seq(1, n) / n)
#' lineages_n <- b_times_n
#' b_times2_n <- b_times_n * b_times_n
#' lineages2_n <- b_times2_n
#'
#' # Calculate nLTT
#' out <- nLTT::nltt_diff_exact_calc_extinct(
#' event_times = b_times_n,
#' species_number = lineages_n,
#' event_times2 = b_times2_n,
#' species_number2 = lineages2_n,
#' distance_method = "abs"
#' )
#' #'
#' @export
nltt_diff_exact_calc_extinct <- function(
event_times,
species_number,
event_times2,
species_number2,
distance_method) {
# Are equally sized?{
if (requireNamespace("testit", quietly = TRUE)) {
testit::assert(length(event_times) == length(species_number))
testit::assert(length(event_times2) == length(species_number2))
# Are non-trivial?
testit::assert(length(event_times) > 0)
testit::assert(length(species_number) > 0)
testit::assert(length(event_times2) > 0)
testit::assert(length(species_number) > 0)
# Are sorted?
testit::assert(all(sort(event_times) == event_times))
testit::assert(all(sort(event_times2) == event_times2))
}
# Make a list of all branching times, and remove duplicates
# TODO: Use more effecient merge, as both ranges are sorted
all_event_times <- unique(sort(c(event_times, event_times2)))
diff <- 0.0
#iterate through all branching times
for (k in 2:length(all_event_times)) {
tim <- all_event_times[k]
#find the index of the first branching time
#that is up to the focal branching time
suppressWarnings(
index1 <- max(which(event_times < tim))
)
suppressWarnings(
index2 <- max(which(event_times2 < tim)) #same for the other tree
)
#find the number of lineages at time "tim" for tree 1
lins1 <- species_number[max(index1, 1)]
#find the number of lineages at time "tim" for tree 2
lins2 <- species_number2[max(index2, 1)]
#the amount of time that this difference persisted
dt <- all_event_times[k] - all_event_times[k - 1]
if (distance_method == "abs") {
diff <- diff + dt * abs(lins1 - lins2) #update the difference
}
if (distance_method == "squ") {
diff <- diff + dt * (lins1 - lins2) * (lins1 - lins2)
}
}
diff
}
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.