From e44194e973ff3e976ef84741a1f56b373463ae23 Mon Sep 17 00:00:00 2001
From: Thomas Guillerme <guillert@tcd.ie>
Date: Tue, 5 Nov 2024 17:12:17 +0000
Subject: [PATCH] still trying to understand dtt

---
 R/dtt.dispRity.tmp.R | 552 ++++++++++++++++++++++++++-----------------
 R/dtt.dispRity_fun.R |   2 +-
 2 files changed, 335 insertions(+), 219 deletions(-)

diff --git a/R/dtt.dispRity.tmp.R b/R/dtt.dispRity.tmp.R
index e1068a28..741ebc16 100644
--- a/R/dtt.dispRity.tmp.R
+++ b/R/dtt.dispRity.tmp.R
@@ -1,198 +1,314 @@
-#' @title dtt dispRity (from \code{geiger::dtt})
-#'
-#' @description A wrapper for the \code{geiger::dtt} function working with any disparity metric.
-#'
-#' @param data A \code{dispRity} object or a \code{matrix}
-#' @param metric The disparity metric to be passed to \code{\link{dispRity}}.
-#' @param tree A \code{phylo} object matching the data and with a \code{root.time} element. Can be missing if \code{data} has a \code{tree} component.
-#' @param nsim The number of simulations to calculate null disparity-through-time.
-#' @param model A evolutionary model for the simulations (see \code{geiger::sim.char} - default is \code{"BM"}).
-#' @param alternative The H1 alternative (for calculating the p-value). Can be \code{"two-sided"} (default), \code{"greater"} or \code{"lesser"}; see details.
-#' @param scale.time Optional, whether to scale the time (between 0 and 1; \code{TRUE}, default) or not (\code{FALSE}).
-#' @param ... Any other arguments to be passed to \code{geiger::dtt}.
-#' 
-#' @details
-#' See \code{geiger::dtt} for details.
-#' 
-#' @examples
-#'
-#' ## Loading morphological data and a tree
-#' data(BeckLee_mat50)
-#' data(BeckLee_tree)
-#' 
-#' ## The average squared pairwise distance metric (used in geiger::dtt)
-#' average.sq <- function(X) mean(pairwise.dist(X)^2)
-#' 
-#' ## Calculate the disparity of the dataset using dtt.dispRity
-#' dispRity_dtt <- dtt.dispRity(data = BeckLee_mat50, metric = average.sq,
-#'                              tree = BeckLee_tree, nsim = 20)
-#' 
-#' ## Plotting the results
-#' plot(dispRity_dtt)
-#' 
-#' 
-#' @seealso
-#' \code{\link{test.dispRity}}, \code{\link{custom.subsets}}, \code{\link{chrono.subsets}}, \code{\link{plot.dispRity}}.
-#' 
-#' @author Thomas Guillerme
-# @export
-
-# source("sanitizing.R")
-# source("dispRity_fun.R")
-# source("dtt.dispRity_fun.R")
+        
+
+#         ## Disparity is calculated for all nodes as the disparity of the descendants
+#         ## But also as the average disparity from all descendants if there is more data on other branches
+#         ## e.g. for the plot(tree), node 17 is average of node 18 and 23, 18 is the average of 19 and 23, etc.
+#         ## This is some kind of acctran chrono.subsets with no node values? - probably not
+# library(geiger); library(dispRity)
 # geiger_data <- get(data(geospiza))
 # data = geiger_data$dat
 # average.sq <- function(X) mean(pairwise.dist(X)^2)
 # metric = average.sq
 # tree = geiger_data$phy
-# tree$root.time <- 5
-# nsim = 100
-# model = "BM"
-# #relative.disp = TRUE
-# alternative = "two-sided"
-# dots <- list()
+# cleaned_data <- clean.data(data, tree)
+# data <- cleaned_data$data
+# tree <- cleaned_data$tree
+# tree$node.label <- paste0(1:Nnode(tree)+Ntip(tree))
 
-# set.seed(1)
-# dispRity_dtt <- dtt.dispRity(data = geiger_data$dat, metric = average.sq, tree = geiger_data$phy, nsim = 100, alternative = "two-sided", ...)
-# dispRity_dtt$p_value
 
-# average.sq <- function(X) mean(pairwise.dist(X)^2)
-# metric = average.sq
-# tree <- rcoal(5)
-# data <- space.maker(5, 4, rnorm)
-# rownames(data) <- tree$tip.label
-
-
-# Modified version of the geiger::dtt function (https://github.com/mwpennell/geiger-v2/blob/master/R/disparity.R)
-dtt.dispRity.tmp <- function(data, metric, tree, nsim = 0, model = "BM", alternative = "two-sided", scale.time = TRUE, ...) {
-
-
-    match_call <- match.call()
-    dots <- list(...)
-
-    ## SANITIZING
-    ## data
-    data_class <- check.class(data, c("dispRity", "matrix"))
-    if(data_class == "dispRity") {
-        data <- data$matrix
-    } else {
-        data <- list(data)
-    }
-
-    ## metric
-    metrics_list <- get.dispRity.metric.handle(metric, match_call)
-
-    ## Tree
-    check.class(tree, "phylo")
-    #TODO: in the future allow for a list of trees
-
-    ## Match the tree to the data
-    cleaned_data <- clean.data(data, tree)
-    data <- cleaned_data$data
-    tree <- cleaned_data$tree
-
-    if(!is.na(cleaned_data$dropped_tips)) {
-        warning("The following tip(s) was not present in the data: ", paste(cleaned_data$dropped_tips, collapse = ", "), ".")
-    }
-    if(!is.na(cleaned_data$dropped_rows)) {
-        warning("The following element(s) was not present in the tree: ", paste(cleaned_data$dropped_rows, collapse = ", "), ".")
-    }
-
-    ## alternative
-    check.method(alternative, c("two-sided", "greater", "lesser"), msg = "alternative")
-    ## Set up the p-value function
-    get.p.value <- switch(alternative,
-        "two-sided" = {
-            function(sim_disparity, obs_disparity, replicates) {
-                ## Centring the randoms and observed
-                center_random <- abs(sim_disparity - mean(sim_disparity, na.rm = TRUE))
-                center_observed <- abs(obs_disparity - mean(sim_disparity, na.rm = TRUE))
-                ## Getting the p
-                return((sum(center_random >= center_observed))/(length(sim_disparity)))
-            }
-        },
-        greater = {
-           function(sim_disparity, obs_disparity, replicates) {
-                # Getting the p
-                return((sum(sim_disparity >= obs_disparity))/(length(sim_disparity)))
-            }
-        },
-        lesser = {
-            function(sim_disparity, obs_disparity, replicates) {
-                # Getting the p
-                return((sum(sim_disparity <= obs_disparity))/(length(sim_disparity)))
-            }
-        }
-    )
-
-    ## nsim
-    check.class(nsim, c("numeric", "integer"))
-    if(nsim > 0) {
-        ## Check the model
-        model_default <- c("BM", "OU")
-    }
-
-
-
-
-
-
-    ## New version
-
-    # Step 1 get clade groups through time
-    clade_groups <- custom.subsets(data, group = tree)
-        # If simulation step
-        ## map the trait on the tree n_times
-
-    get.nodes.disparity <- function(data, tree, metric) {
+
+#         ## Basic disparity
+#         disparity_per_node <- unlist(get.disparity(dispRity(custom.subsets(data, group = tree), metric = average.sq)))
+
+#         plot(tree, show.tip.label = FALSE)
+#         nodelabels()
+#         axisPhylo()
+#         # dev.new()
+
+
+#         # plot(NULL, xlim = c(0.6,0), ylim = c(0,1))
+
+#         node_ages <- tree.age(tree)[-c(1:Ntip(tree)), 1]
+#         names(node_ages) <- tree.age(tree)[-c(1:Ntip(tree)), 2]
+#         # node_ages <- (1-rev(nodes_ages))
+
+#         abline(v = max(node_ages) - node_ages, lwd = 0.5, col = "grey")
+#         ordered_nodes <- sort(node_ages, decreasing = TRUE)
+#         plot_nodes <- (max(node_ages) - node_ages)[names(ordered_nodes)]
+#         axis(2)
+
+#         ## Basic sub-clade diversity
+#         disparity <- unlist(get.disparity(dispRity(custom.subsets(data, tree), metric = metric)))
+#         relative_disparity <- disparity/max(disparity)
+#         relative_disparity <- relative_disparity * 12
+#         points(y = relative_disparity[names(ordered_nodes)], x = plot_nodes, col = "blue")
+
+
+
+#     phy <- tree
+#     phy$node.label <- NULL
+#     td <- list("phy" = phy, "data" = data) #TG: data is already cleaned prior to the function
+#     phy2 <- td$phy
+#     ## Converting to old phylo format
+#     phy <- new2old.phylo(td$phy)
+
+#     ## Getting the node depth from ape if no attributes
+#     if(is.null(phy2$root.time)){
+#         node.depth <- branching.times(phy2)
+#     } else {
+#         node.depth <- tree.age(phy2)$ages[-c(1:Ntip(phy2))]
+#         names(node.depth) <- 1:length(node.depth) + Ntip(phy2)
+#     }
+
+#     stem.depth <- numeric()
+#     stem.depth[1] <- node.depth[1]
+#     for(i in 2:phy2$Nnode) {
+#         anc <- which(as.numeric(phy$edge[,2]) == -i)
+#         stem.depth[i] <- node.depth[names(node.depth) == phy2$edge[anc,1]]
+#     }
+
+#     ## Lineages through time
+#     ltt <- sort(node.depth, decreasing = TRUE)
+
+#     ## Scaling by lineage through time
+#     # node.depth <- node.depth/max(ltt)
+#     # stem.depth <- stem.depth/max(ltt)
+#     # ltt <- ltt/max(ltt)
+
+#     result <- numeric()
+
+#     ## By matrix
+#     if(length(dim(td$data)) == 2) {
+
+#         ## Calculate disparity per clade
+#         disparity <- as.vector(summary(dispRity(custom.subsets(td$data, phy2), metric = metric), digits = 10)$obs)
+#         names(disparity) <- 1:length(node.depth) + Ntip(phy2)
         
-        ## Basic disparity
-        disparity_per_node <- unlist(get.disparity(dispRity(custom.subsets(data, group = tree), metric = average.sq)))
+#         ## Disparity at the root
+#         result[1] <- disparity[1]
+
+#         ## Disparity for the other nodes (average per time slice)
+#         for(i in 2:length(ltt)) {
+#             x <- disparity[stem.depth >= ltt[i-1] & node.depth < ltt[i-1]]
+#             ## Checks: which stem.depth greater= than previous lineage through time
+#             ##         and   node.depth smaller  than previous lineage through time
+
+#             #result[i] <- ifelse(length(x) == 0, 0, mean(x))                             #TODO: check 0 for non-fossil tree
+#             if(length(x) == 0) {
+#                 result[i] <- 0
+#             } else {
+#                 result[i] <- mean(x)
+#             }
+#         }
+
+#         result[length(ltt)+1] <- 0                                                      #TODO: check 0 for non-fossil tree
+#         names(result) <- names(disparity)
+
+#         # ## Scaling the results
+#         # if(result[1] > 0){                                                              #TODO: check 0 for non-fossil tree
+#         #     result <- result/result[1]
+#         # }
+
+#     } 
+
+#     result_plot <- result/max(disparity)*12
+
+
+#     lines(y = result_plot[-length(result)][names(ordered_nodes)], x = plot_nodes, col = "orange")
+
+
+#     mean(c(disparity["14"])) == result["14"]
+#     mean(c(disparity["15"])) == result["15"]
+#     mean(c(disparity["16"])) == result["16"]
+#     mean(c(disparity["17"])) == result["17"]
+#     mean(c(disparity["18"], disparity["23"])) == result["18"]
+#     mean(c(disparity["19"], disparity["23"])) == result["19"]
+#     mean(c(disparity["20"], disparity["23"])) == result["20"]
+#     mean(c(disparity["21"], disparity["23"])) == result["21"]
+#     mean(c(disparity["22"], disparity["23"])) == result["22"]
+#     mean(c(disparity["22"], disparity["24"])) == result["23"]
+#     mean(c(disparity["24"])) == result["24"]
+#     mean(c(disparity["25"])) == result["25"]
+
+
+#     ## TODO: do an option like: @param use.average logical, whether to use the average clade disparity (\code{TRUE}; as used in \code{geiger::dtt}) or just the clade disparity (\code{FALSE}). See details.
+
+#     #@details When using \code{use.average = FALSE}, disparity is calculated as the disparity of each clade (i.e. the disparity for each node and their unique descendants). When using \code{use.average = TRUE}, disparity is calculated 
+
+
+
+
+
+#         ## Combining the tree and the data
+#         phy <- tree
+
+#         ## Get the node depth
+#         node.depth <- tree.age(phy2)$ages[-c(1:Ntip(phy2))]
+#         names(node.depth) <- 1:length(node.depth) + Ntip(phy2)
+    
+    
+#         stem.depth <- numeric()
+#         for(one_node in 1:Nnode(phy)) {
+#             ## Get the ancestor
+#             anc <- which(phy$edge[,2] == one_node+Ntip(phy)+1)
+#             stem.depth[one_node] <- node.depth[names(node.depth) == phy2$edge[anc,1]]
+#         }
+
+#         ## Lineages through time
+#         ltt <- sort(node.depth, decreasing = TRUE)
+
+#         ## Scaling by lineage through time
+#         node.depth <- node.depth/max(ltt)
+#         stem.depth <- stem.depth/max(ltt)
+#         ltt <- ltt/max(ltt)
+
+#         result <- numeric()
+
+#         ## Calculate disparity per clade
+
+
+
 
-        ## Disparity is calculated for all nodes as the disparity of the descendants
-        ## But also as the average disparity from all descendants if there is more data on other branches
-        ## e.g. for the plot(tree), node 17 is average of node 18 and 23, 18 is the average of 19 and 23, etc.
-        ## This is some kind of acctran chrono.subsets with no node values? - probably not
 
+#         ## Disparity at the root
+#         result[1] <- disparity[1]
+#         backup_check <- list()
+#         backup_check[[1]] <- disparity[1]
 
+#         ## Disparity for the other nodes (average per time slice)
+#         for(i in 2:length(ltt)) {
+#             x <- disparity[stem.depth >= ltt[i-1] & node.depth < ltt[i-1]]
+#             #result[i] <- ifelse(length(x) == 0, 0, mean(x))                             #TODO: check 0 for non-fossil tree
+#             backup_check[[i]] <- x
+#             if(length(x) == 0) {
+#                 result[i] <- 0
+#             } else {
+#                 result[i] <- mean(x)
+#             }
+#         }
 
 
-        if(!node.clade.disparity) {
-            ## Basic disparity
-            return(disparity_per_node)
-        } else {
-            ## Detect each node that has a pair of nodes as descendants
-            which(tree$edge[, ]
 
-            tree$edge[which(!(tree$edge[,2] %in% 1:Ntip(tree))), ]
 
-            ## Get the disparity for all nodes
-            disparity_per_node 
-        }
-    }
+# #        14        15        16        17        18        19        20        21 
+# # 0.7198663 0.6992722 0.4539000 0.4677565 0.5668206 0.5991425 0.7182857 0.9987066 
+# #        22        23        24        25 
+# # 0.4368961 0.1490411 0.1522438 0.2928359 
+
+# # [[1]]
+# #        14 
+# # 0.7198663 
+
+# # [[2]]
+# # named numeric(0)
+
+# # [[3]]
+# # named numeric(0)
+
+# # [[4]]
+# # named numeric(0)
+
+# # [[5]]
+# #        15 
+# # 0.6992722 
+
+# # [[6]]
+# #        17 
+# # 0.4677565 
+
+# # [[7]]
+# #        18        23 
+# # 0.5668206 0.1490411 # this is node 18 + 23
+
+# # [[8]]
+# #        18        23 
+# # 0.5668206 0.1490411 
+
+# # [[9]]
+# #        20        23 
+# # 0.7182857 0.1490411 
+
+# # [[10]]
+# #        21        23 
+# # 0.9987066 0.1490411 
+
+# # [[11]]
+# #        22        23 
+# # 0.4368961 0.1490411 
+
+# # [[12]]
+# #        22        24 
+# # 0.4368961 0.1522438 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#         if(!node.clade.disparity) {
+#             ## Basic disparity
+#             return(disparity_per_node)
+#         } else {
+#             ## Detect each node that has a pair of nodes as descendants
+#             which(tree$edge[, ]
+
+#             tree$edge[which(!(tree$edge[,2] %in% 1:Ntip(tree))), ]
+
+#             ## Get the disparity for all nodes
+#             disparity_per_node 
+#         }
+#     }
+
+
+# geiger_data <- get(data(geospiza))
+# data = geiger_data$dat
+# average.sq <- function(X) mean(pairwise.dist(X)^2)
+# metric = average.sq
+# tree = geiger_data$phy
+# tree$root.time <- 5
 
 
-    # Step 2 measure disparity for that group through time
-    clade_disparity <- dispRity(clade_groups, metric = average.sq)
-    disp_clade <- unlist(get.disparity(clade_disparity))
-        # If simulation step
-        ## measure disparity on the simulations
+#     # Step 2 measure disparity for that group through time
+#     clade_disparity <- dispRity(clade_groups, metric = average.sq)
+#     disp_clade <- unlist(get.disparity(clade_disparity))
+#         # If simulation step
+#         ## measure disparity on the simulations
 
-# To investigate patterns of disparity through time, we moved up the phylogeny from the root. At each divergence event (i.e., each node), we calculated the mean relative disparity for that point in time as the average of the relative disparities of all subclades whose ancestral lineages were present at that time. Values near 0 imply that subclades contain relatively little of the variation present within the taxon as a whole and that, consequently, most variation is partitioned as among-subclade differences; conversely, values near 1 imply that subclades contain a substantial proportion of the total variation and thus are likely to overlap extensively, indicating that subclades have independently evolved to occupy similar regions of morphological space. In this way, we plotted an average disparity through time analogous to the plots of lineage accumulation through time
+# # To investigate patterns of disparity through time, we moved up the phylogeny from the root. At each divergence event (i.e., each node), we calculated the mean relative disparity for that point in time as the average of the relative disparities of all subclades whose ancestral lineages were present at that time. Values near 0 imply that subclades contain relatively little of the variation present within the taxon as a whole and that, consequently, most variation is partitioned as among-subclade differences; conversely, values near 1 imply that subclades contain a substantial proportion of the total variation and thus are likely to overlap extensively, indicating that subclades have independently evolved to occupy similar regions of morphological space. In this way, we plotted an average disparity through time analogous to the plots of lineage accumulation through time
 
-# To calculate how much mean disparity differed from that expected under a null hypothesis of character evolution by unconstrained Brownian motion, we conducted 1000 simulations of morphological diversification on each taxon's phylogeny (25). [(We also conducted simulations with a speciational model of character evolution; results were qualitatively unchanged (fig. S5).] Variances of each morphological axis among species in simulations were set equal to the actual variances of the principal components for each taxon. The morphological disparity index (MDI), the overall difference in relative disparity of a clade compared with that expected under the null hypothesis, was calculated as the area contained between the line connecting observed relative disparity points versus the line connecting median relative disparity points of the simulations; areas in which observed values were above expected were given positive values, whereas those below expected were given negative values. Results from this analysis were corroborated with the δ statistic, which tests for acceleration or deceleration of morphological change through time against a null model of constant rates of character evolution (26).
+# # To calculate how much mean disparity differed from that expected under a null hypothesis of character evolution by unconstrained Brownian motion, we conducted 1000 simulations of morphological diversification on each taxon's phylogeny (25). [(We also conducted simulations with a speciational model of character evolution; results were qualitatively unchanged (fig. S5).] Variances of each morphological axis among species in simulations were set equal to the actual variances of the principal components for each taxon. The morphological disparity index (MDI), the overall difference in relative disparity of a clade compared with that expected under the null hypothesis, was calculated as the area contained between the line connecting observed relative disparity points versus the line connecting median relative disparity points of the simulations; areas in which observed values were above expected were given positive values, whereas those below expected were given negative values. Results from this analysis were corroborated with the δ statistic, which tests for acceleration or deceleration of morphological change through time against a null model of constant rates of character evolution (26).
 
-    # If simulation step
-        ## measure the alternative
+#     # If simulation step
+#         ## measure the alternative
 
 
 
 
-    ## Create the subsets for each clade and calculate disparity
-clade_disparity <- dispRity(clade_groups, metric = average.sq)
-disp_clade <- unlist(get.disparity(clade_disparity))
-disp_clade <- disp_clade/max(disp_clade) 
-lines(x = dtt_disparity$times[-c(1)], y = disp_clade, col = "blue", lwd = 2)
+#     ## Create the subsets for each clade and calculate disparity
+# clade_disparity <- dispRity(clade_groups, metric = average.sq)
+# disp_clade <- unlist(get.disparity(clade_disparity))
+# disp_clade <- disp_clade/max(disp_clade) 
+# lines(x = dtt_disparity$times[-c(1)], y = disp_clade, col = "blue", lwd = 2)
 
 
 
@@ -200,69 +316,69 @@ lines(x = dtt_disparity$times[-c(1)], y = disp_clade, col = "blue", lwd = 2)
 
 
 
-    ## mdi.range
-    if(is.null(dots$mdi.range)) {
-        dots$mdi.range <- c(0,1)
-    }
+#     ## mdi.range
+#     if(is.null(dots$mdi.range)) {
+#         dots$mdi.range <- c(0,1)
+#     }
 
-    ## Get the scaled disparity through time
-    disparity_through_time1 <- geiger.dtt.dispRity(tree, data, metric)
-    disparity_through_time2 <- c(unname(unlist(get.disparity(scale.dispRity(dispRity(custom.subsets(data, tree), metric = metric))))), 0)
+#     ## Get the scaled disparity through time
+#     disparity_through_time1 <- geiger.dtt.dispRity(tree, data, metric)
+#     disparity_through_time2 <- c(unname(unlist(get.disparity(scale.dispRity(dispRity(custom.subsets(data, tree), metric = metric))))), 0)
 
 
-    ## Get the lineages through time
-    lineage_through_time <- sort(branching.times(tree), decreasing = TRUE)
-    if(scale.time) {
-        lineage_through_time <- c(0, (max(lineage_through_time)-lineage_through_time)/max(lineage_through_time))
-    }
+#     ## Get the lineages through time
+#     lineage_through_time <- sort(branching.times(tree), decreasing = TRUE)
+#     if(scale.time) {
+#         lineage_through_time <- c(0, (max(lineage_through_time)-lineage_through_time)/max(lineage_through_time))
+#     }
 
-    ## Simulating the null disparity through time
-    if(is.numeric(nsim) && nsim > 0){
+#     ## Simulating the null disparity through time
+#     if(is.numeric(nsim) && nsim > 0){
 
-        ## Calculating the rate matrix
-        rate_matrix <- geiger.ratematrix(tree, data)
+#         ## Calculating the rate matrix
+#         rate_matrix <- geiger.ratematrix(tree, data)
 
-        ## Simulate the data
-        #TODO: Eventually replace this with dads!
-        simulated_data <- geiger.sim.char(tree, rate_matrix, nsim, model = model)
+#         ## Simulate the data
+#         #TODO: Eventually replace this with dads!
+#         simulated_data <- geiger.sim.char(tree, rate_matrix, nsim, model = model)
 
-        disparity_through_time_sim <- geiger.dtt.dispRity(tree, simulated_data, metric)
+#         disparity_through_time_sim <- geiger.dtt.dispRity(tree, simulated_data, metric)
 
-        # ## Convert into a list
-        # simulated_data <- lapply(seq(dim(simulated_data)[3]), function(x) simulated_data[ , , x])
+#         # ## Convert into a list
+#         # simulated_data <- lapply(seq(dim(simulated_data)[3]), function(x) simulated_data[ , , x])
 
-        # ## Calculating the disparity
-        # disparity_through_time_sim <- lapply(simulated_data, function(simulated_data, tree, metric) .dtt.dispRity(tree, simulated_data, metric), tree, metric)
-        # disparity_through_time_sim <- matrix(unlist(disparity_through_time_sim), ncol = nsim, byrow = FALSE)
+#         # ## Calculating the disparity
+#         # disparity_through_time_sim <- lapply(simulated_data, function(simulated_data, tree, metric) .dtt.dispRity(tree, simulated_data, metric), tree, metric)
+#         # disparity_through_time_sim <- matrix(unlist(disparity_through_time_sim), ncol = nsim, byrow = FALSE)
 
-        colnames(disparity_through_time_sim) <- NULL
+#         colnames(disparity_through_time_sim) <- NULL
 
-        ## MDI
-        MDI <- unname(geiger.area.between.curves(lineage_through_time, apply(disparity_through_time_sim, 1, median, na.rm = TRUE), disparity_through_time, sort(dots$mdi.range)))
+#         ## MDI
+#         MDI <- unname(geiger.area.between.curves(lineage_through_time, apply(disparity_through_time_sim, 1, median, na.rm = TRUE), disparity_through_time, sort(dots$mdi.range)))
 
-        ## Get the simulated MDIs
-        sim_MDI <- apply(disparity_through_time_sim, 2, function(X) geiger.area.between.curves(x = lineage_through_time, f1 = X, f2 = disparity_through_time))
+#         ## Get the simulated MDIs
+#         sim_MDI <- apply(disparity_through_time_sim, 2, function(X) geiger.area.between.curves(x = lineage_through_time, f1 = X, f2 = disparity_through_time))
 
-        ## Calculate the p_value
-        p_value <- get.p.value(sim_MDI, MDI)
+#         ## Calculate the p_value
+#         p_value <- get.p.value(sim_MDI, MDI)
 
-        ## Sort the output
-        output <- list(dtt = disparity_through_time, times = lineage_through_time, sim = disparity_through_time_sim, MDI = MDI, sim_MDI = sim_MDI, p_value = p_value, call = match_call)
+#         ## Sort the output
+#         output <- list(dtt = disparity_through_time, times = lineage_through_time, sim = disparity_through_time_sim, MDI = MDI, sim_MDI = sim_MDI, p_value = p_value, call = match_call)
 
-        ## Add the model (if it was used as default)
-        if(is.null(output$call$model)) {
-            output$call$model <- model
-        }
-        if(is.null(output$call$alternative)) {
-            output$call$alternative <- alternative
-        }
+#         ## Add the model (if it was used as default)
+#         if(is.null(output$call$model)) {
+#             output$call$model <- model
+#         }
+#         if(is.null(output$call$alternative)) {
+#             output$call$alternative <- alternative
+#         }
 
-    } else {
+#     } else {
 
-        output <- list(dtt = disparity_through_time, times = lineage_through_time)
-    }
+#         output <- list(dtt = disparity_through_time, times = lineage_through_time)
+#     }
 
-    class(output) <- c("dispRity", "dtt")
+#     class(output) <- c("dispRity", "dtt")
 
-    return(output)
-}
+#     return(output)
+# }
diff --git a/R/dtt.dispRity_fun.R b/R/dtt.dispRity_fun.R
index 163e2d38..710d250c 100755
--- a/R/dtt.dispRity_fun.R
+++ b/R/dtt.dispRity_fun.R
@@ -1,6 +1,6 @@
 
 ## Modified .dtt function from https://github.com/mwpennell/geiger-v2/blob/master/R/disparity.R
-geiger.dtt.dispRity <- function(phy, data, metric, relative){
+geiger.dtt.dispRity <- function(phy, data, metric){
 
     ## Combining the tree and the data
     phy$node.label <- NULL