#' @title dispRity object summary
#'
#' @description Creates a summary of a \code{dispRity} object.
#'
#' @param object A \code{dispRity} object.
#' @param ... Additional arguments to be passed to \code{\link{summary}} or \code{cent.tend}.
#' @param quantiles The quantiles to display (default is \code{quantiles = c(50, 95)}; is ignored if the \code{dispRity} object is not bootstrapped).
#' @param cent.tend A function for summarising the bootstrapped disparity values (default is \code{\link[stats]{median}}).
#' @param recall \code{logical} value specifying whether to recall the \code{dispRity} parameters input (default = \code{FALSE}).
#' @param digits Optional, a value for digits the values in the output table (default = 2).
# ' @param results Optional, in the case of summarising a \code{\link{sequential.test}} which results to display (default = "coefficients")
#'
#' @return
#' A \code{data.frame} with:
#' \item{subsets}{the subset names.}
#' \item{n}{the maximum number of elements in each subset (see details).}
#' \item{observed}{the observed disparity or the the observed central tendency (<cent_tend>) of disparity (\code{obs.<cent_tend>}).}
#' \item{bootstraps...}{if \code{data} is bootstrapped, the bootstrapped disparity's central tendency (\code{bs.<cent_tend>}) and the quantiles of the bootstrapped disparities (or, if \code{data} is not bootstrapped but disparity is calculated as a distribution - see \code{\link[dispRity]{dispRity}}) - the quantiles of the observed disparity are displayed).}
#'
#' @details
#' If the \code{dispRity} object to summarise comes from a \code{\link[dispRity]{chrono.subsets}} using a \code{"multiPhylo"} object, the displayed number of observations (\code{n}) corresponds to the maximum number of observation at the specific time slice (some slices through some trees might have less observations).
#'
#' @examples
#' ## Load the disparity data based on Beck & Lee 2014
#' data(disparity)
#'
#' ## Summarising the results
#' summary(disparity) # default
#' ## Using different options
#' summary(disparity, quantiles = 75, cent.tend = mean, digits = 8,
#' recall = TRUE)
#'
#' @seealso \code{\link{dispRity}}, \code{\link{plot.dispRity}}.
#'
#' @author Thomas Guillerme
## DEBUG
# source("sanitizing.R")
# source("make.metric.R")
# source("summary.dispRity_fun.R")
# data(BeckLee_mat50)
# groups <- as.data.frame(matrix(data = c(rep(1, 12), rep(2, 13), rep(3, 12), rep(4, 13)), dimnames = list(rownames(BeckLee_mat50))), ncol = 1)
# customised_subsets <- custom.subsets(BeckLee_mat50, groups)
# bootstrapped_data <- boot.matrix(customised_subsets, bootstraps = 3, rarefaction = TRUE)
# subsets <- extract.dispRity(sum_of_variances, observed = FALSE, keep.structure = TRUE, concatenate = TRUE)
# data <- sequential.test(subsets, family = gaussian, correction = "hommel")
# data <- dispRity(bootstrapped_data, metric = variances)
# quantiles <- c(50, 95)
# cent.tend <- median
# recall <- FALSE
# match_call <- list() ; match_call$cent.tend <- "median"
summary.dispRity <- function(object, ..., quantiles = c(50, 95), cent.tend = median, recall = FALSE, digits){#, results = "coefficients") {
## Renaming object
data <- object
rm(object)
#----------------------
# SANITIZING
#----------------------
#Get call
match_call <- match.call()
# warning("DEBUG summary.dispRity") ; return(match_call)
#cent.tend
#Must be a function
check.class(cent.tend, "function")
#The function must work
if(make.metric(cent.tend, silent = TRUE)$type != "level1") {
stop.call(match_call$cent.tend, " cannot be used for measuring the central tendency.")
}
## Update match_call if argument is empty
if(is.null(match_call$cent.tend)) match_call$cent.tend <- "median"
#recall
check.class(recall, "logical")
#digits
if(missing(digits)) {
#Set to default (see below)
digits <- "default"
} else {
check.class(digits, c("numeric", "integer"))
}
#DATA
#must be class dispRity
check.class(data, "dispRity")
#Check if it is a bootstrapped dispRity object
if(is.null(data$disparity) && is.na(class(data)[2])) {
stop.call("", "Disparity has not been calculated yet.\nUse the dispRity() function to do so.")
}
#Check quantiles
check.class(quantiles, "numeric", " must be any value between 1 and 100.")
if(any(quantiles < 1) | any(quantiles > 100)) {
stop.call("", "quantiles(s) must be any value between 1 and 100.")
}
#----------------------
# SPECIAL SHORTCUTS (dual classes)
#----------------------
if(length(class(data)) > 1) {
## Model test summary
if(is(data, "model.test")) {
## Extracting the AICs and the log likelihoods
base_results <- cbind(data$aic.models, "log.lik" = sapply(data$full.details, function(x) x$value))
## Extracting the additional parameters
parameters <- lapply(data$full.details, function(x) x$par)
# MP: allow summaries to work on a single model
if(length(parameters) == 1) {
param.tmp <- c(parameters[[1]])
#names(param.tmp) <- rownames(parameters[[1]])
parameters <- list(param.tmp)
}
base_results <- cbind(base_results, "param" = unlist(lapply(parameters, length)))
## Get the full list of parameters
names_list <- lapply(parameters, names)
full_param <- unique(unlist(names_list))
output_table <- cbind(base_results, do.call(rbind, lapply(parameters, match.parameters, full_param)))
## Rounding
summary_results <- digits.fun(output_table, digits, model.test = TRUE)
return(summary_results)
}
## Model sim summary
if(is(data, "model.sim")) {
# if(recall){
# print.dispRity(data)
# }
## Extract the central tendencies
simulation_data_matrix <- sapply(data$simulation.data$sim, function(x) x$central_tendency)
## Get the quantiles
simulation_results <- apply(simulation_data_matrix, 1, get.summary, cent.tend = cent.tend, quantiles = quantiles, ...)
simulation_results <- cbind(do.call(rbind, lapply(simulation_results, function(X) rbind(X$cent_tend[[1]]))),
do.call(rbind, lapply(simulation_results, function(X) rbind(X$quantiles))))
colnames(simulation_results)[1] <- as.character(match_call$cent.tend)
## Output table
output_table <- cbind("subsets" = rev(data$simulation.data$fix$subsets),
"n" = data$simulation.data$fix$sample_size,
"var" = unname(data$simulation.data$fix$variance),
simulation_results)
rownames(output_table) <- seq(1:nrow(output_table))
return(output_table)
}
## test.metric summary
if(is(data, "test.metric")) {
## Summary table
get.cent.tends <- function(results) {
combined <- aggregate(disparity~reduction, data = results, FUN = function(x, cent.tend) cent.tend(x, na.rm = TRUE), cent.tend, na.action = na.pass)
return(combined[order(as.numeric(combined$reduction)),2])
}
all_results <- round(do.call(rbind, lapply(data$results, get.cent.tends)), digits = ifelse(digits == "default", 2, digits))
colnames(all_results) <- paste0(unique(data$results[[1]][,"reduction"]), "%")
## Adding the model results
if(!is.null(data$models)) {
get.model.summary <- function(model) {
## Try to get the model parameters
try_slope <- try.get.from.model(model, "Estimate")
try_p_val <- try.get.from.model(model, "Pr\\(")
try_r_squ <- try.get.from.model(model, "r.squared")
if(!is.null(try_r_squ) && any(names(try_r_squ) == "adj.r.squared")) {
try_r_squ <- try_r_squ$adj.r.squared
is_adjusted <- TRUE
} else {
is_adjusted <- FALSE
}
## Make a results vector
results <- c(
"slope" = ifelse(any(!is.null(try_slope)) && any(!is.na(try_slope)), try_slope["reduction"], NA),
"p_val" = ifelse(any(!is.null(try_p_val)) && any(!is.na(try_p_val)), try_p_val["reduction"], NA),
"r_squ" = ifelse(any(!is.null(try_r_squ)) && any(!is.na(try_r_squ)), try_r_squ, NA))
names(results) <- c("slope", "p_value", ifelse(is_adjusted, "R^2(adj)", "R^2"))
return(results)
}
## Get the model summaries
model_summaries <- do.call(rbind, lapply(data$models, get.model.summary))
## Combine both
all_results <- cbind(all_results, model_summaries)
}
return(all_results)
}
## randtest summary
if(is(data, "randtest")) {
table <- make.randtest.table(data)
return(round(table, digits = ifelse(digits== "default", 2, digits)))
}
## axes summary
if(is(data, "axes")) {
if(is(digits, "character") && digits == "default") {
digits <- 3
}
## Get the column names
if(is.null(data$call$colnames)) {
col_names <- paste0("dim.", 1:length(data$scaled.var[[1]]))
} else {
col_names <- data$call$colnames
}
## Get each variance component
scaled_var <- do.call(rbind, data$scaled.var)
colnames(scaled_var) <- paste0(col_names, ".var")
cum_var <- do.call(rbind, data$cumsum.var)
colnames(cum_var) <- paste0(col_names, ".sum")
## Get the combined variance components
combined_var <- cbind(scaled_var, cum_var)
## Combine them by order
base_col <- c(0, ncol(scaled_var))
col_order <- numeric()
for(i in 1:ncol(scaled_var)) {
col_order <- c(col_order, base_col+i)
}
return(round(combined_var[, col_order, drop = FALSE], digits = digits))
}
## projection summary
if(is(data, "projection")) {
## Summarise each object
return(lapply(data, summary.dispRity, quantiles = quantiles, cent.tend = cent.tend, recall = FALSE))
}
if(is(data, "pgls.dispRity")) {
## SUmmarise pgls.dispRity output into phylolm
return(convert.to.summary.phylolm(data))
}
## No dual class summary available
stop.call("", paste0("No specific summary for combined class \"dispRity\" and \"", class(data)[2], "\"."))
}
#----------------------
# TRANSFORMING THE DATA INTO A TABLE
#----------------------
## Check if disparity is a value or a distribution
is_distribution <- ifelse(length(data$disparity[[1]]$elements) != 1, TRUE, FALSE)
## Check the bootstraps
bootstrapped <- !is.null(data$call$bootstrap) && !(data$call$bootstrap[[2]] == "covar")
## Get the elements per subsets
elements <- lapply(data$subsets, lapply.get.elements, bootstrapped)
nulls <- unlist(lapply(elements, is.null))
if(any(nulls)) {
for(null_elem in which(nulls)) {
elements[[null_elem]] <- nrow(data$subsets[[null_elem]]$elements)
}
}
## Get the names of the subsets
names <- names(data$subsets)
if(is.null(names)) {
names <- seq(1:length(data$subsets))
}
## Get the disparity values
disparity_values <- lapply(data$disparity, lapply.observed)
names(disparity_values) <- NULL
## Change the names and subsets if between groups
if(data$call$disparity$metrics$between.groups) {
names(elements) <- gsub(":", ";", names(elements))
## Groups names
names <- names(data$disparity)
## Remove warning (desired: matching rows with different numbers of elements!)
options(warn = -1)
## Get the table of elements for each pair of groups
elements_matrix <-lapply(as.list(names), function(name, elements)
do.call(cbind, lapply(as.list(strsplit(name, split = ":")[[1]]), function(name, elements) return(unname(elements[[match(name, names(elements))]])), elements)),
elements = elements)
options(warn = 0)
## Initialise the results
summary_results <- data.frame(row.names = NULL,
"subsets" = rep(names, unlist(lapply(elements_matrix, nrow))),
"n_1" = do.call(rbind,elements_matrix)[,1],
"n_2" = do.call(rbind,elements_matrix)[,2])
} else {
## Initialise the results
summary_results <- data.frame(row.names = NULL, "subsets" = rep(names, unlist(lapply(elements, length))), "n" = unlist(elements))
}
## Set the observed colume
if(data$call$disparity$metrics$between.groups) {
obs_col <- 4
elements <- lapply(elements_matrix, function(x) seq(1:nrow(x)))
} else {
obs_col <- 3
}
## Add the observed values
if(is_distribution) {
summary_results <- cbind(summary_results, as.vector(unlist(mapply(mapply.observed, lapply(disparity_values, cent.tend, ...), elements))), row.names = NULL)
names(summary_results)[obs_col] <- paste("obs", as.expression(match_call$cent.tend), sep = ".")
} else {
summary_results <- cbind(summary_results, as.vector(unlist(mapply(mapply.observed, disparity_values, elements))), row.names = NULL)
names(summary_results)[obs_col] <- "obs"
}
if(bootstrapped) {
## Calculate the central tendencies and the quantiles
summary_results <- cbind(summary_results, matrix(unlist(lapply(data$disparity, lapply.summary, cent.tend, quantiles, ...)), byrow = TRUE, ncol = (1+length(quantiles)*2)))
## Adding the labels
names(summary_results)[(obs_col+1):length(summary_results)] <- c(paste("bs", as.expression(match_call$cent.tend), sep = "."), names(quantile(rnorm(5), probs = CI.converter(quantiles))))
} else {
if(is_distribution) {
## Calculate the quantiles
summary_results <- cbind(summary_results, matrix(unlist(lapply(data$disparity, lapply, get.summary, quantiles = quantiles)), byrow = TRUE, ncol = (length(quantiles)*2)))
## Adding the labels
names(summary_results)[(obs_col+1):length(summary_results)] <- c(names(quantile(rnorm(5), probs = CI.converter(quantiles))))
}
}
## Round the results (number of decimals = maximum number of digits in the output)
summary_results <- digits.fun(summary_results, digits)
## If any elements is equal to one, check if not NA
if(data$call$disparity$metrics$between.groups) {
col_n <- c(2,3)
} else {
col_n <- 2
}
if(any(summary_results[, col_n] == 1)) {
## Select the values to check
to_check <- which(summary_results[, col_n] == 1)
## Get their replacement values
check.elements.NA <- function(row, summary_results, data) {
## Check if the subset contains NA elements
ifelse(is.na(data$subsets[[as.character(summary_results[row, 1])]]$elements[1,1]), 0, 1)
}
replace_vals <- sapply(to_check, check.elements.NA, summary_results, data)
## Replace them in the results
summary_results[to_check, col_n] <- replace_vals
}
#----------------------
# OUTPUT
#----------------------
if(recall) print.dispRity(data)
rm(data)
return(summary_results)
#Summary sequential.test shortcut
# if(length(class(data)) == 2) {
# if(class(data)[[1]] == "dispRity" && class(data)[[2]] == "seq.test") {
# #Results sanitizing
# check.class(results, "character")
# #At least one must be "coefficients"
# if(is.na(match("coefficients", results))) {
# stop("At least one of the returned results must be 'coefficients'.")
# }
# #Results must be at least coefficients
# if(is.na(match("coefficients", results))) {
# results <- c(results, "coefficients")
# }
# #Creating the table results
# results_out <- summary.seq.test(data, quantiles, cent.tend, recall, digits, results, match_call)
# #Checking if distribution
# is.distribution <- ifelse(length(data$models[[1]]) == 1, FALSE, TRUE)
# #Rounding the results
# if(is.distribution == FALSE) {
# results_out <- lapply(results_out, digits.fun, digits, seq.test = TRUE)
# } else {
# results_out <- lapply(results_out, lapply, digits.fun, digits, seq.test = TRUE)
# }
# return(results_out)
# }
# if(class(data)[1] == "dispRity" & class(data)[2] == "randtest") {
# #No summary
# return(data)
# }
# }
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.