diff --git a/DESCRIPTION b/DESCRIPTION index 49167a6..e17dd75 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: divent Type: Package Title: Entropy Partitioning to Measure Diversity -Version: 0.4-99.9013 +Version: 0.4-99.9014 Authors@R: c( person("Eric", "Marcon", email="eric.marcon@agroparistech.fr", role=c("aut", "cre"), comment=c(ORCID="0000-0002-5249-321X")) ) diff --git a/NEWS.md b/NEWS.md index bb2af51..3e55066 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# divent 0.4-99.9013 +# divent 0.4-99.9014 ## Features diff --git a/R/accum_sp.R b/R/accum_sp.R index b9f3521..2ec4f9c 100644 --- a/R/accum_sp.R +++ b/R/accum_sp.R @@ -333,3 +333,84 @@ plot_map <- function( # Return the image for further processing return(invisible(the_image)) } + + +#' Helper to prepare parameters for `accum_sp` plot and autoplot. +#' +#' @param x the `accum_sp` object to plot. +#' @param q the order of diversity. +#' @param main the title of the plot. +#' @param xlab X-axis label. +#' @param ylab Y-axis label. +#' @param ylim Y-axis limits +#' +#' @returns a vector of parameters for the plots +#' @noRd +#' +accum_sp_plot_helper <- function(x, q, main, xlab, ylab, ylim) { + + # Find the row in the accumulation table + q_row <- which(dimnames(x$accumulation)$q == q) + if (length(q_row) != 1) { + stop("The value of q does not correspond to any accumulation curve.") + } + + if (is.null(ylim)) { + # Evaluate ylim if not set by an argument + ymin <- min(x$accumulation[q_row, , ]) + ymax <- max(x$accumulation[q_row, , ]) + } else { + ymin <- ylim[1] + ymax <- ylim[2] + } + + if (main == "Accumulation of ...") { + # Prepare the main title + if (inherits(x, "accum_sp_entropy")) { + main <- paste("Accumulation of Entropy of order", q) + } + if (inherits(x, "accum_sp_diversity")) { + if (q == 0) { + main <- "Species Accumulation Curve" + } else { + main <- paste("Accumulation of Diversity of order", q) + } + } + if (inherits(x, "accum_sp_mixing")) main <- paste("Mixing index of order", q) + } + + if (xlab == "Sample size...") { + if (names(dimnames(x$accumulation)[2]) == "n") { + xlab <- "Number of individuals" + } else { + xlab <- "Distance from the central individual" + } + } + + if (ylab == "Diversity...") { + # Prepare Y-axis + if (inherits(x, "accum_sp_entropy")) { + ylab <- "Diversity" + } + if (inherits(x, "accum_sp_diversity")) { + if (q == 0) { + ylab <- "Richness" + } else { + ylab <- "Diversity" + } + if (inherits(x, "accum_sp_mixing")) { + ylab <- "Mixing index" + } + } + return( + list( + q_row = q_row, + ymin = ymin, + ymax = ymax, + main = main, + xlab = xlab, + ylab = ylab + ) + ) + } +} diff --git a/R/div_hurlbert.R b/R/div_hurlbert.R index e72171e..feadd8b 100644 --- a/R/div_hurlbert.R +++ b/R/div_hurlbert.R @@ -117,3 +117,35 @@ div_hurlbert.species_distribution <- function( return(the_diversity) } + + +#' Compute Hurlbert's diversity from its entropy +#' +#' Find the effective number of species numerically +#' +#' @param hurlbert_entropy The entropy. +#' @param k The order of entropy. +#' +#' @returns Hurlbert's effective number of species. +#' @noRd +#' +hurlbert_ent2div <- function(hurlbert_entropy, k) { + # Relation between diversity and entropy + # (D for diversity, S for entropy, k is the parameter) + f <- function(D, S, k) {D * (1 - (1 - 1 / D)^k) - S} + # Minimize it + return( + vapply( + hurlbert_entropy, + FUN = function(S) { + stats::uniroot( + f = f, + interval = c(1, 1E+7), + S = S, + k = k + )$root + }, + FUN.VALUE = 0 + ) + ) +} diff --git a/R/ent_phylo.R b/R/ent_phylo.R index a339f70..2b9ef42 100644 --- a/R/ent_phylo.R +++ b/R/ent_phylo.R @@ -205,3 +205,109 @@ ent_phylo.species_distribution <- function( } } } + + +#' Phylogenetic entropies +#' +#' Calculate entropies of a list of phylogenetic abundances (obtained by +#' [phylo_abd]). +#' Each item of the list corresponds to a phylogenetic group, i.e. an interval +#' of the tree (where the species do not change). +#' +#' @param phylo_abd A list of matrices of abundance (caution: lines are species, +#' columns are communities). +#' +#' @returns A vector. Each item is the entropy of a community. +#' +#' @noRd +#' +phylo_entropy.phylo_abd <- function( + # Allow lapply along q + q, + phylo_abd, + tree, + normalize, + # Other arguments for ent_tsallis.numeric + estimator, + level, + probability_estimator, + unveiling, + richness_estimator, + jack_alpha, + jack_max, + coverage_estimator, + gamma) { + + if (gamma) { + # Calculate gamma entropy of each group. + # simplify2array() makes a vector with the list of numbers. + phylo_entropies <- simplify2array( + lapply( + # Calculate entropy in each item of the list, i.e. group. + # Obtain a list. + phylo_abd, + FUN = function(group) { + ent_tsallis.species_distribution( + as_abundances(t(group)), + q = q, + estimator = estimator, + level = level, + probability_estimator = probability_estimator, + unveiling = unveiling, + richness_estimator = richness_estimator, + jack_alpha = jack_alpha, + jack_max = jack_max, + coverage_estimator = coverage_estimator, + gamma = TRUE, + # Obtain a vector. + as_numeric = TRUE, + check_arguments = FALSE + ) + } + ) + ) + + } else { + # Calculate entropy of each community in each group. + # simplify2array() makes a matrix with the list of vectors. + phylo_entropies <- simplify2array( + lapply( + # Calculate entropy in each item of the list, i.e. group. + # Obtain a list. + phylo_abd, + FUN = function(group) { + apply( + group, + # Calculate entropy of each column of the matrix, i.e. community. + MARGIN = 2, + FUN = ent_tsallis.numeric, + # Arguments + q = q, + estimator = estimator, + level = level, + probability_estimator = probability_estimator, + unveiling = unveiling, + richness_estimator = richness_estimator, + jack_alpha = jack_alpha, + jack_max = jack_max, + coverage_estimator = coverage_estimator, + # Obtain a vector. + as_numeric = TRUE, + check_arguments = FALSE + ) + } + ) + ) + } + # Should be a matrix, but simplify2array() makes a vector instead of a 1-col + # matrix and gamma entropy is a vector. Force a matrix. + if (is.vector(phylo_entropies)) { + phylo_entropies <- t(phylo_entropies) + } + + # Calculate the weighted mean of entropy and normalize + the_entropy <- as.numeric(tree$intervals %*% t(phylo_entropies)) + if (normalize) the_entropy <- the_entropy / sum(tree$intervals) + + return(the_entropy) +} diff --git a/R/ent_similarity.R b/R/ent_similarity.R index 61b14e5..6444ec6 100644 --- a/R/ent_similarity.R +++ b/R/ent_similarity.R @@ -442,3 +442,81 @@ S_v <- function( v_used <- seq_len(sample_size - abd[species_index]) return(sum(w_v[v_used] * p_V_Ns[v_used, species_index])) } + + +#' Similarity-Based Gamma entropy of a metacommunity +#' +#' Build the metacommunity and check that abundances are integers. +#' +#' See [ent_gamma_tsallis] for details. +#' @param species_distribution An object of class [species_distribution]. +#' +#' @returns A tibble with the estimator used and the estimated entropy. +#' @noRd +#' +ent_gamma_similarity <- function( + species_distribution, + similarities, + q, + estimator, + probability_estimator, + unveiling, + jack_alpha, + jack_max, + coverage_estimator, + as_numeric) { + + # Build the metacommunity + abd <- metacommunity.abundances( + species_distribution, + as_numeric = TRUE, + check_arguments = FALSE + ) + if (is_integer_values(abd)) { + # Sample coverage is useless + sample_coverage <- NULL + } else { + # Non-integer values in the metacommunity. + # Calculate the sample coverage and change the estimator. + sample_coverage <- coverage.numeric( + colSums( + species_distribution[ + , !colnames(species_distribution) %in% non_species_columns + ] + ), + estimator = coverage_estimator, + as_numeric = TRUE, + check_arguments = FALSE + ) + if (!estimator %in% c("Marcon", "ChaoShen")) { + estimator <- "Marcon" + } + } + + # Compute the entropy. + the_entropy <- ent_similarity.numeric( + abd, + similarities = similarities, + q = q, + estimator = estimator, + probability_estimator = probability_estimator, + unveiling = unveiling, + jack_alpha = jack_alpha, + jack_max = jack_max, + sample_coverage = sample_coverage, + as_numeric = as_numeric, + check_arguments = FALSE + ) + + # Return + if (as_numeric) { + return(the_entropy) + } else { + return( + dplyr::bind_cols( + site = "Metacommunity", + the_entropy + ) + ) + } +} diff --git a/R/ent_sp_simpson.R b/R/ent_sp_simpson.R index d80a80a..dc64cf6 100644 --- a/R/ent_sp_simpson.R +++ b/R/ent_sp_simpson.R @@ -189,3 +189,23 @@ ent_sp_simpsonEnvelope <- function( # Return the envelope return(the_envelope) } + + +#' Extract a column from an fv object +#' according to an edge-effect correction +#' +#' @param fv the function value object, see [spatstat.explore::fv.object]. +#' @param correction the edge-effect correction: +#' "isotropic", "translate" or "none" +#' +#' @returns a vector with the function values +#' @noRd +#' +correction_fv <- function(fv, correction) { + switch( + correction, + "isotropic" = fv$iso, + "translate" = fv$trans, + "none" = fv$un + ) +} diff --git a/R/package.R b/R/package.R index 13667e5..7598efe 100644 --- a/R/package.R +++ b/R/package.R @@ -992,128 +992,6 @@ check_divent_args <- function( } -#' Helper to prepare parameters for `accum_sp` plot and autoplot. -#' -#' @param x the `accum_sp` object to plot. -#' @param q the order of diversity. -#' @param main the title of the plot. -#' @param xlab X-axis label. -#' @param ylab Y-axis label. -#' @param ylim Y-axis limits -#' -#' @returns a vector of parameters for the plots -#' @noRd -#' -accum_sp_plot_helper <- function(x, q, main, xlab, ylab, ylim) { - - # Find the row in the accumulation table - q_row <- which(dimnames(x$accumulation)$q == q) - if (length(q_row) != 1) { - stop("The value of q does not correspond to any accumulation curve.") - } - - if (is.null(ylim)) { - # Evaluate ylim if not set by an argument - ymin <- min(x$accumulation[q_row, , ]) - ymax <- max(x$accumulation[q_row, , ]) - } else { - ymin <- ylim[1] - ymax <- ylim[2] - } - - if (main == "Accumulation of ...") { - # Prepare the main title - if (inherits(x, "accum_sp_entropy")) { - main <- paste("Accumulation of Entropy of order", q) - } - if (inherits(x, "accum_sp_diversity")) { - if (q == 0) { - main <- "Species Accumulation Curve" - } else { - main <- paste("Accumulation of Diversity of order", q) - } - } - if (inherits(x, "accum_sp_mixing")) main <- paste("Mixing index of order", q) - } - - if (xlab == "Sample size...") { - if (names(dimnames(x$accumulation)[2]) == "n") { - xlab <- "Number of individuals" - } else { - xlab <- "Distance from the central individual" - } - } - - if (ylab == "Diversity...") { - # Prepare Y-axis - if (inherits(x, "accum_sp_entropy")) { - ylab <- "Diversity" - } - if (inherits(x, "accum_sp_diversity")) { - if (q == 0) { - ylab <- "Richness" - } else { - ylab <- "Diversity" - } - if (inherits(x, "accum_sp_mixing")) { - ylab <- "Mixing index" - } - } - return( - list( - q_row = q_row, - ymin = ymin, - ymax = ymax, - main = main, - xlab = xlab, - ylab = ylab - ) - ) - } -} - - - -#' as_named_vector.character -#' -#' Counts the number of points of a `character` vector and returns a named vector. -#' Names are the items of the character vector. -#' This is equivalent to `as.numeric(table(x))` but `table()` -#' looses the names. -#' -#' @param x a character vector. -#' -#' @returns A named vector with the number of items by name. -#' @noRd -#' -as_named_vector.character <- function(x){ - # Count the number of items. Returns a 1D array, not a vector. - the_array <- tapply(x, x, length) - the_vector <- as.vector(the_array) - # Add the names - names(the_vector) <- names(the_array) - return(the_vector) -} - - -#' as_named_vector.wmppp -#' -#' Counts the number of points of a `wmppp` object and returns a named vector. -#' Names are the point types. -#' This is equivalent to `as.numeric(table(X$marks$PointType))` but `table()` -#' looses the names. -#' -#' @param X a [dbmss::wmppp] object, i.e. a weighted, marked planar point pattern. -#' -#' @returns A named vector with the number of points by type. -#' @noRd -#' -as_named_vector.wmppp <- function(X){ - # Count the number of points by type - return(as_named_vector.character(spatstat.geom::marks(X)$PointType)) -} - - #' Chao's A #' #' Helper for Chao's estimators. @@ -1351,84 +1229,6 @@ ent_gamma_tsallis <- function( } -#' Similarity-Based Gamma entropy of a metacommunity -#' -#' Build the metacommunity and check that abundances are integers. -#' -#' See [ent_gamma_tsallis] for details. -#' @param species_distribution An object of class [species_distribution]. -#' -#' @returns A tibble with the estimator used and the estimated entropy. -#' @noRd -#' -ent_gamma_similarity <- function( - species_distribution, - similarities, - q, - estimator, - probability_estimator, - unveiling, - jack_alpha, - jack_max, - coverage_estimator, - as_numeric) { - - # Build the metacommunity - abd <- metacommunity.abundances( - species_distribution, - as_numeric = TRUE, - check_arguments = FALSE - ) - if (is_integer_values(abd)) { - # Sample coverage is useless - sample_coverage <- NULL - } else { - # Non-integer values in the metacommunity. - # Calculate the sample coverage and change the estimator. - sample_coverage <- coverage.numeric( - colSums( - species_distribution[ - , !colnames(species_distribution) %in% non_species_columns - ] - ), - estimator = coverage_estimator, - as_numeric = TRUE, - check_arguments = FALSE - ) - if (!estimator %in% c("Marcon", "ChaoShen")) { - estimator <- "Marcon" - } - } - - # Compute the entropy. - the_entropy <- ent_similarity.numeric( - abd, - similarities = similarities, - q = q, - estimator = estimator, - probability_estimator = probability_estimator, - unveiling = unveiling, - jack_alpha = jack_alpha, - jack_max = jack_max, - sample_coverage = sample_coverage, - as_numeric = as_numeric, - check_arguments = FALSE - ) - - # Return - if (as_numeric) { - return(the_entropy) - } else { - return( - dplyr::bind_cols( - site = "Metacommunity", - the_entropy - ) - ) - } -} - - #' Error message #' #' Utility for [check_divent_args] @@ -1448,58 +1248,6 @@ error_message <- function(message, argument, parent_function) { } -#' Extract a column from an fv object -#' according to an edge-effect correction -#' -#' @param fv the function value object, see [spatstat.explore::fv.object]. -#' @param correction the edge-effect correction: -#' "isotropic", "translate" or "none" -#' -#' @returns a vector with the function values -#' @noRd -#' -correction_fv <- function(fv, correction) { - switch( - correction, - "isotropic" = fv$iso, - "translate" = fv$trans, - "none" = fv$un - ) -} - - -#' Compute Hurlbert's diversity from its entropy -#' -#' Find the effective number of species numerically -#' -#' @param hurlbert_entropy The entropy. -#' @param k The order of entropy. -#' -#' @returns Hurlbert's effective number of species. -#' @noRd -#' -hurlbert_ent2div <- function(hurlbert_entropy, k) { - # Relation between diversity and entropy - # (D for diversity, S for entropy, k is the parameter) - f <- function(D, S, k) {D * (1 - (1 - 1 / D)^k) - S} - # Minimize it - return( - vapply( - hurlbert_entropy, - FUN = function(S) { - stats::uniroot( - f = f, - interval = c(1, 1E+7), - S = S, - k = k - )$root - }, - FUN.VALUE = 0 - ) - ) -} - - #' Check Integers #' #' Check that the values of a vector are integer, whatever their type. @@ -1552,109 +1300,3 @@ phylo_abd <- function( simplify = FALSE ) } - - -#' Phylogenetic entropies -#' -#' Calculate entropies of a list of phylogenetic abundances (obtained by -#' [phylo_abd]). -#' Each item of the list corresponds to a phylogenetic group, i.e. an interval -#' of the tree (where the species do not change). -#' -#' @param phylo_abd A list of matrices of abundance (caution: lines are species, -#' columns are communities). -#' -#' @returns A vector. Each item is the entropy of a community. -#' -#' @noRd -#' -phylo_entropy.phylo_abd <- function( - # Allow lapply along q - q, - phylo_abd, - tree, - normalize, - # Other arguments for ent_tsallis.numeric - estimator, - level, - probability_estimator, - unveiling, - richness_estimator, - jack_alpha, - jack_max, - coverage_estimator, - gamma) { - - if (gamma) { - # Calculate gamma entropy of each group. - # simplify2array() makes a vector with the list of numbers. - phylo_entropies <- simplify2array( - lapply( - # Calculate entropy in each item of the list, i.e. group. - # Obtain a list. - phylo_abd, - FUN = function(group) { - ent_tsallis.species_distribution( - as_abundances(t(group)), - q = q, - estimator = estimator, - level = level, - probability_estimator = probability_estimator, - unveiling = unveiling, - richness_estimator = richness_estimator, - jack_alpha = jack_alpha, - jack_max = jack_max, - coverage_estimator = coverage_estimator, - gamma = TRUE, - # Obtain a vector. - as_numeric = TRUE, - check_arguments = FALSE - ) - } - ) - ) - - } else { - # Calculate entropy of each community in each group. - # simplify2array() makes a matrix with the list of vectors. - phylo_entropies <- simplify2array( - lapply( - # Calculate entropy in each item of the list, i.e. group. - # Obtain a list. - phylo_abd, - FUN = function(group) { - apply( - group, - # Calculate entropy of each column of the matrix, i.e. community. - MARGIN = 2, - FUN = ent_tsallis.numeric, - # Arguments - q = q, - estimator = estimator, - level = level, - probability_estimator = probability_estimator, - unveiling = unveiling, - richness_estimator = richness_estimator, - jack_alpha = jack_alpha, - jack_max = jack_max, - coverage_estimator = coverage_estimator, - # Obtain a vector. - as_numeric = TRUE, - check_arguments = FALSE - ) - } - ) - ) - } - # Should be a matrix, but simplify2array() makes a vector instead of a 1-col - # matrix and gamma entropy is a vector. Force a matrix. - if (is.vector(phylo_entropies)) { - phylo_entropies <- t(phylo_entropies) - } - - # Calculate the weighted mean of entropy and normalize - the_entropy <- as.numeric(tree$intervals %*% t(phylo_entropies)) - if (normalize) the_entropy <- the_entropy / sum(tree$intervals) - - return(the_entropy) -} diff --git a/R/species_distribution.R b/R/species_distribution.R index 40e5a03..fe9690f 100644 --- a/R/species_distribution.R +++ b/R/species_distribution.R @@ -699,3 +699,45 @@ as.double.species_distribution <- function(x, use.names = TRUE, ...) { as.numeric.species_distribution <- function(x, use.names = TRUE, ...) { return(as.double.species_distribution(x, use.names, ...)) } + + +# Utilities ---- + +#' as_named_vector.character +#' +#' Counts the number of points of a `character` vector and returns a named vector. +#' Names are the items of the character vector. +#' This is equivalent to `as.numeric(table(x))` but `table()` +#' looses the names. +#' +#' @param x a character vector. +#' +#' @returns A named vector with the number of items by name. +#' @noRd +#' +as_named_vector.character <- function(x){ + # Count the number of items. Returns a 1D array, not a vector. + the_array <- tapply(x, x, length) + the_vector <- as.vector(the_array) + # Add the names + names(the_vector) <- names(the_array) + return(the_vector) +} + + +#' as_named_vector.wmppp +#' +#' Counts the number of points of a `wmppp` object and returns a named vector. +#' Names are the point types. +#' This is equivalent to `as.numeric(table(X$marks$PointType))` but `table()` +#' looses the names. +#' +#' @param X a [dbmss::wmppp] object, i.e. a weighted, marked planar point pattern. +#' +#' @returns A named vector with the number of points by type. +#' @noRd +#' +as_named_vector.wmppp <- function(X){ + # Count the number of points by type + return(as_named_vector.character(spatstat.geom::marks(X)$PointType)) +} diff --git a/man/div_similarity.Rd b/man/div_similarity.Rd index 56ae063..6d12bf3 100644 --- a/man/div_similarity.Rd +++ b/man/div_similarity.Rd @@ -37,6 +37,7 @@ div_similarity(x, similarities, q = 1, ...) jack_max = 10, coverage_estimator = c("ZhangHuang", "Chao", "Turing", "Good"), gamma = FALSE, + as_numeric = FALSE, ..., check_arguments = TRUE ) diff --git a/man/ent_rao.Rd b/man/ent_rao.Rd index 96403a4..5a4b8a3 100644 --- a/man/ent_rao.Rd +++ b/man/ent_rao.Rd @@ -26,6 +26,7 @@ ent_rao(x, ...) normalize = TRUE, estimator = c("Lande", "naive"), gamma = FALSE, + as_numeric = FALSE, ..., check_arguments = TRUE ) diff --git a/man/ent_similarity.Rd b/man/ent_similarity.Rd index 8a1fb52..cfb6be6 100644 --- a/man/ent_similarity.Rd +++ b/man/ent_similarity.Rd @@ -37,6 +37,7 @@ ent_similarity(x, similarities, q = 1, ...) jack_max = 10, coverage_estimator = c("ZhangHuang", "Chao", "Turing", "Good"), gamma = FALSE, + as_numeric = FALSE, ..., check_arguments = TRUE )