R/hill_func.R
In daijiang/hillR: Diversity Through Hill Numbers
Defines functions
hill_func
Documented in
hill_func
#' Functional diversity through Hill Numbers
#'
#' Calculate functional diversity for each site (alpha diversity).
#'
#' @inheritParams hill_taxa
#' @param traits A data frame of species functional traits data. Species as rows, traits as columns.
#' It can include both continuous and categorical data. It will be transformed into a distance
#' matrix using `FD::gowdis(traits)`. If all traits are numeric, then it will use Euclidean distance.
#' @param traits_as_is if \code{FALSE} (default) traits data frame will be transformed into a distance
#' matrix. Otherwise, will use as is (i.e. traits is a symmetric distance matrix).
#' @param check_data whether to check data first? Default is \code{TRUE}.
#' @param div_by_sp as FD calculated in this way will be highly correlated with taxonomic diversity,
#' one potential simple way to correct this is to divide the results by the number of species.
#' However, a more common way to deal with correlations is to use null models and calculate standardized effect sizes.
#' Therefore, I set the default to be \code{FALSE}.
#' @param ord ord in \code{FD::gowdis}.
#' @param fdis whether to calculate FDis, default is \code{TRUE}
#' @param stand_dij whether to standardize distance matrix to have max value of 1? Default is \code{FALSE}.
#' @export
#' @references Chao, Anne, Chun-Huo Chiu, and Lou Jost. Unifying Species Diversity, Phylogenetic Diversity, Functional Diversity, and Related Similarity and Differentiation Measures Through Hill Numbers. Annual Review of Ecology, Evolution, and Systematics 45, no. 1 (2014): 297–324. <doi:10.1146/annurev-ecolsys-120213-091540>.
#'
#' Chiu, Chun-Huo, and Anne Chao. Distance-Based Functional Diversity Measures and Their Decomposition: A Framework Based on Hill Numbers. PLoS ONE 9, no. 7 (July 7, 2014): e100014. <doi:10.1371/journal.pone.0100014>.
#' @return A matrix, with these information for each site: Q (Rao's Q); D_q (functional hill number,
#' the effective number of equally abundant and functionally equally distinct species);
#' MD_q (mean functional diversity per species, the effective sum of pairwise distances between
#' a fixed species and all other species); FD_q (total functional diversity, the effective total functional
#' distance between species of the assemblage). See Chiu and Chao 2014 page 4 for more information.
#'
#' @examples
#' dummy = FD::dummy
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 0)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 1)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 0.9999)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 2)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 3)
#'
hill_func <- function(comm, traits, traits_as_is = FALSE, q = 0, base = exp(1), check_data = TRUE,
div_by_sp = FALSE, ord = c("podani", "metric"), fdis = TRUE, stand_dij = FALSE) {
if (check_data) {
if (any(comm < 0))
stop("Negative value in comm data")
if (any(colSums(comm) == 0))
warning("Some species in comm data were not observed in any site,\n delete them...")
if (any(rowSums(comm) == 0))
warning("Some sites in comm data do not have any species,\n delete them...")
comm <- comm[rowSums(comm) != 0, colSums(comm) != 0, drop = FALSE]
if(traits_as_is){
if(is.null(attributes(traits)$Labels)) stop("\n Traits distance matrix has no labels\n")
} else {
if (is.null(rownames(traits))) stop("\n Traits have no row names\n")
}
if (is.null(colnames(comm))) {
stop("\n Comm data have no col names\n")
}
}
if(traits_as_is){
trait_sp = attributes(traits)$Labels
} else {
trait_sp = rownames(traits)
}
if (any(!colnames(comm) %in% trait_sp)) {
warning("\n There are species from community data that are not on traits matrix\n
Delete these species from comm data...\n")
if(any(!trait_sp %in% colnames(comm))){
warning("\n There are species in the traits data not in the communit data\n
Delete these species from trait data...\n")
trait_sp = trait_sp[trait_sp %in% colnames(comm)]
}
comm <- comm[, trait_sp]
}
# all(rownames(traits) == names(comm))
if (traits_as_is) {
# traits is already a distance matrix
dij <- traits
if(!inherits(dij, "dist")) stop("`traits` is not a distance object yet `trait_as_is` is TRUE\n")
} else {
# traits is not a distance matrix
traits <- traits[trait_sp, , drop = FALSE]
if (ncol(traits) == 1) {
# only 1 trait
if (any(is.na(traits))) {
warning("Warning: Species with missing trait values have been excluded.",
"\n")
traits <- na.omit(traits)
comm <- comm[, colnames(comm) %in% rownames(traits)]
}
if (is.numeric(traits[, 1])) {
# 1 numeric trait
dij <- dist(traits)
}
if (is.factor(traits[, 1]) | is.character(traits[, 1])) {
# 1 categorical trait
if (is.ordered(traits[, 1])) {
traits2 <- data.frame(rank(traits[, 1]))
rownames(traits2) <- rownames(traits)
names(traits2) <- names(traits)
dij <- dist(traits2)
} else {
traits[, 1] <- as.factor(traits[, 1])
x.f <- as.factor(traits[, 1])
x.dummy <- diag(nlevels(x.f))[x.f, ]
x.dummy.df <- data.frame(x.dummy, row.names = rownames(traits))
dij <- ade4::dist.binary(x.dummy.df, method = 2)
}
}
} else {
# more than 1 trait:
for (i in 1:ncol(traits)) {
if (is.factor(traits[, i]) & nlevels(traits[, i]) == 2) {
traits[, i] <- as.numeric(traits[, i]) - 1 # so to be 0, 1
}
}
if (all(sapply(traits, is.numeric)) & all(!is.na(traits))) {
dij <- dist(scale(traits, center = TRUE, scale = TRUE))
} else {
ord <- match.arg(ord)
dij <- FD::gowdis(x = traits, asym.bin = NULL, ord = ord)
}
# dij = gowdis(x=traits, ...)
}
# if (!is.euclid(dij)) { if (corr == 'lingoes') { dij2 <- lingoes(dij)
# warning('Species x species distance matrix was not Euclidean. Lingoes correction was
# applied.','\n') } if (corr == 'cailliez') { dij2 <- cailliez(dij) warning('Species
# x species distance matrix was not Euclidean. Cailliez correction was
# applied.','\n') } if (corr == 'sqrt') { dij2 <- sqrt(dij) # check if sqrt
# correction actually worked if(!is.euclid(dij2) ) stop('Species x species distance
# matrix was still is not Euclidean after 'sqrt' correction. Use another correction
# method.','\n') if (is.euclid(dij2) ) warning('Species x species distance matrix was
# not Euclidean. 'sqrt' correction was applied.','\n') } if (corr == 'none') { dij2
# <- quasieuclid(dij) warning('Species x species distance was not Euclidean, but no
# correction was applied. Only the PCoA axes with positive eigenvalues were
# kept.','\n') } dij = dij2 }
}
if (fdis) {
# calculate fdis
FDis <- FD::fdisp(d = dij, a = as.matrix(comm))$FDis
}
comm <- as.matrix(comm)
N <- nrow(comm)
S <- ncol(comm)
SR <- rowSums(comm > 0) # species richness of each site
if(q == 0) comm[comm > 0] = 1 # so that quadratic entropy Q does not use pi pj
comm <- sweep(comm, 1, rowSums(comm, na.rm = TRUE), "/") # relative abun
dij <- as.matrix(dij)
if(all(is.na(dij[upper.tri(dij)])))
stop("All pairwise distance is NA, do all species have the same trait values?")
if (stand_dij)
dij <- dij/max(dij)
# inter = comm %*% dij # \sum_i,j_S(p_i * dij) Q = rowSums(sweep(comm,1,inter,'*',
# check.margin = F))/2 # \sum_j_S\sum_i,j_S(p_i * dij)
Q <- vector("numeric", length = N)
names(Q) <- dimnames(comm)[[1]]
for (k in 1:N) {
Q[k] <- (comm[k, ]) %*% dij %*% matrix(comm[k, ], ncol = 1) # /2
# Q[k] = sum(dij * outer(comm[k,], comm[k,], '*'))/2
}
## D_q
FD_q <- MD_q <- D_q <- vector("numeric", length = N)
names(D_q) <- dimnames(comm)[[1]]
names(MD_q) <- dimnames(comm)[[1]]
names(FD_q) <- dimnames(comm)[[1]]
if (q == 0) {
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
D_q[k] <- sum(dis2/Q[k])^0.5
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
} else {
if (q == 1) {
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
D_q[k] <- exp(-0.5 * sum(dis2/Q[k] * outer(df2, df2, FUN = "*") * log(outer(df2,
df2, FUN = "*"), base)))
# exp(-0.5 * (dij/Q) * pi*pj * log(pi*pj)
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
} else {
# q != 0 or 1
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
din <- sum(dis2/Q[k] * (outer(df2, df2, FUN = "*")^q))
if (din == 0) {
D_q[k] <- 0
} else {
D_q[k] <- din^(1/(2 * (1 - q)))
}
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
}
}
if (fdis) {
if (div_by_sp == TRUE) {
return(rbind(Q, FDis, D_q/SR, MD_q/SR, FD_q/choose(SR, 2)))
} else {
return(rbind(Q, FDis, D_q, MD_q, FD_q))
}
} else {
if (div_by_sp == TRUE) {
return(rbind(Q, D_q/SR, MD_q/SR, FD_q/choose(SR, 2)))
} else {
return(rbind(Q, D_q, MD_q, FD_q))
}
}
}
daijiang/hillR documentation built on April 28, 2024, 1:45 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions hill_func
Documented in hill_func
#' Functional diversity through Hill Numbers
#'
#' Calculate functional diversity for each site (alpha diversity).
#'
#' @inheritParams hill_taxa
#' @param traits A data frame of species functional traits data. Species as rows, traits as columns.
#' It can include both continuous and categorical data. It will be transformed into a distance
#' matrix using `FD::gowdis(traits)`. If all traits are numeric, then it will use Euclidean distance.
#' @param traits_as_is if \code{FALSE} (default) traits data frame will be transformed into a distance
#' matrix. Otherwise, will use as is (i.e. traits is a symmetric distance matrix).
#' @param check_data whether to check data first? Default is \code{TRUE}.
#' @param div_by_sp as FD calculated in this way will be highly correlated with taxonomic diversity,
#' one potential simple way to correct this is to divide the results by the number of species.
#' However, a more common way to deal with correlations is to use null models and calculate standardized effect sizes.
#' Therefore, I set the default to be \code{FALSE}.
#' @param ord ord in \code{FD::gowdis}.
#' @param fdis whether to calculate FDis, default is \code{TRUE}
#' @param stand_dij whether to standardize distance matrix to have max value of 1? Default is \code{FALSE}.
#' @export
#' @references Chao, Anne, Chun-Huo Chiu, and Lou Jost. Unifying Species Diversity, Phylogenetic Diversity, Functional Diversity, and Related Similarity and Differentiation Measures Through Hill Numbers. Annual Review of Ecology, Evolution, and Systematics 45, no. 1 (2014): 297–324. <doi:10.1146/annurev-ecolsys-120213-091540>.
#'
#' Chiu, Chun-Huo, and Anne Chao. Distance-Based Functional Diversity Measures and Their Decomposition: A Framework Based on Hill Numbers. PLoS ONE 9, no. 7 (July 7, 2014): e100014. <doi:10.1371/journal.pone.0100014>.
#' @return A matrix, with these information for each site: Q (Rao's Q); D_q (functional hill number,
#' the effective number of equally abundant and functionally equally distinct species);
#' MD_q (mean functional diversity per species, the effective sum of pairwise distances between
#' a fixed species and all other species); FD_q (total functional diversity, the effective total functional
#' distance between species of the assemblage). See Chiu and Chao 2014 page 4 for more information.
#'
#' @examples
#' dummy = FD::dummy
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 0)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 1)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 0.9999)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 2)
#' hill_func(comm = dummy$abun, traits = dummy$trait, q = 3)
#'
hill_func <- function(comm, traits, traits_as_is = FALSE, q = 0, base = exp(1), check_data = TRUE,
div_by_sp = FALSE, ord = c("podani", "metric"), fdis = TRUE, stand_dij = FALSE) {
if (check_data) {
if (any(comm < 0))
stop("Negative value in comm data")
if (any(colSums(comm) == 0))
warning("Some species in comm data were not observed in any site,\n delete them...")
if (any(rowSums(comm) == 0))
warning("Some sites in comm data do not have any species,\n delete them...")
comm <- comm[rowSums(comm) != 0, colSums(comm) != 0, drop = FALSE]
if(traits_as_is){
if(is.null(attributes(traits)$Labels)) stop("\n Traits distance matrix has no labels\n")
} else {
if (is.null(rownames(traits))) stop("\n Traits have no row names\n")
}
if (is.null(colnames(comm))) {
stop("\n Comm data have no col names\n")
}
}
if(traits_as_is){
trait_sp = attributes(traits)$Labels
} else {
trait_sp = rownames(traits)
}
if (any(!colnames(comm) %in% trait_sp)) {
warning("\n There are species from community data that are not on traits matrix\n
Delete these species from comm data...\n")
if(any(!trait_sp %in% colnames(comm))){
warning("\n There are species in the traits data not in the communit data\n
Delete these species from trait data...\n")
trait_sp = trait_sp[trait_sp %in% colnames(comm)]
}
comm <- comm[, trait_sp]
}
# all(rownames(traits) == names(comm))
if (traits_as_is) {
# traits is already a distance matrix
dij <- traits
if(!inherits(dij, "dist")) stop("`traits` is not a distance object yet `trait_as_is` is TRUE\n")
} else {
# traits is not a distance matrix
traits <- traits[trait_sp, , drop = FALSE]
if (ncol(traits) == 1) {
# only 1 trait
if (any(is.na(traits))) {
warning("Warning: Species with missing trait values have been excluded.",
"\n")
traits <- na.omit(traits)
comm <- comm[, colnames(comm) %in% rownames(traits)]
}
if (is.numeric(traits[, 1])) {
# 1 numeric trait
dij <- dist(traits)
}
if (is.factor(traits[, 1]) | is.character(traits[, 1])) {
# 1 categorical trait
if (is.ordered(traits[, 1])) {
traits2 <- data.frame(rank(traits[, 1]))
rownames(traits2) <- rownames(traits)
names(traits2) <- names(traits)
dij <- dist(traits2)
} else {
traits[, 1] <- as.factor(traits[, 1])
x.f <- as.factor(traits[, 1])
x.dummy <- diag(nlevels(x.f))[x.f, ]
x.dummy.df <- data.frame(x.dummy, row.names = rownames(traits))
dij <- ade4::dist.binary(x.dummy.df, method = 2)
}
}
} else {
# more than 1 trait:
for (i in 1:ncol(traits)) {
if (is.factor(traits[, i]) & nlevels(traits[, i]) == 2) {
traits[, i] <- as.numeric(traits[, i]) - 1 # so to be 0, 1
}
}
if (all(sapply(traits, is.numeric)) & all(!is.na(traits))) {
dij <- dist(scale(traits, center = TRUE, scale = TRUE))
} else {
ord <- match.arg(ord)
dij <- FD::gowdis(x = traits, asym.bin = NULL, ord = ord)
}
# dij = gowdis(x=traits, ...)
}
# if (!is.euclid(dij)) { if (corr == 'lingoes') { dij2 <- lingoes(dij)
# warning('Species x species distance matrix was not Euclidean. Lingoes correction was
# applied.','\n') } if (corr == 'cailliez') { dij2 <- cailliez(dij) warning('Species
# x species distance matrix was not Euclidean. Cailliez correction was
# applied.','\n') } if (corr == 'sqrt') { dij2 <- sqrt(dij) # check if sqrt
# correction actually worked if(!is.euclid(dij2) ) stop('Species x species distance
# matrix was still is not Euclidean after 'sqrt' correction. Use another correction
# method.','\n') if (is.euclid(dij2) ) warning('Species x species distance matrix was
# not Euclidean. 'sqrt' correction was applied.','\n') } if (corr == 'none') { dij2
# <- quasieuclid(dij) warning('Species x species distance was not Euclidean, but no
# correction was applied. Only the PCoA axes with positive eigenvalues were
# kept.','\n') } dij = dij2 }
}
if (fdis) {
# calculate fdis
FDis <- FD::fdisp(d = dij, a = as.matrix(comm))$FDis
}
comm <- as.matrix(comm)
N <- nrow(comm)
S <- ncol(comm)
SR <- rowSums(comm > 0) # species richness of each site
if(q == 0) comm[comm > 0] = 1 # so that quadratic entropy Q does not use pi pj
comm <- sweep(comm, 1, rowSums(comm, na.rm = TRUE), "/") # relative abun
dij <- as.matrix(dij)
if(all(is.na(dij[upper.tri(dij)])))
stop("All pairwise distance is NA, do all species have the same trait values?")
if (stand_dij)
dij <- dij/max(dij)
# inter = comm %*% dij # \sum_i,j_S(p_i * dij) Q = rowSums(sweep(comm,1,inter,'*',
# check.margin = F))/2 # \sum_j_S\sum_i,j_S(p_i * dij)
Q <- vector("numeric", length = N)
names(Q) <- dimnames(comm)[[1]]
for (k in 1:N) {
Q[k] <- (comm[k, ]) %*% dij %*% matrix(comm[k, ], ncol = 1) # /2
# Q[k] = sum(dij * outer(comm[k,], comm[k,], '*'))/2
}
## D_q
FD_q <- MD_q <- D_q <- vector("numeric", length = N)
names(D_q) <- dimnames(comm)[[1]]
names(MD_q) <- dimnames(comm)[[1]]
names(FD_q) <- dimnames(comm)[[1]]
if (q == 0) {
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
D_q[k] <- sum(dis2/Q[k])^0.5
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
} else {
if (q == 1) {
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
D_q[k] <- exp(-0.5 * sum(dis2/Q[k] * outer(df2, df2, FUN = "*") * log(outer(df2,
df2, FUN = "*"), base)))
# exp(-0.5 * (dij/Q) * pi*pj * log(pi*pj)
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
} else {
# q != 0 or 1
for (k in 1:N) {
df2 <- comm[k, ][comm[k, ] > 0]
dis2 <- dij[names(df2), names(df2)]
if (Q[k] == 0) {
D_q[k] <- 0
} else {
din <- sum(dis2/Q[k] * (outer(df2, df2, FUN = "*")^q))
if (din == 0) {
D_q[k] <- 0
} else {
D_q[k] <- din^(1/(2 * (1 - q)))
}
}
MD_q[k] <- D_q[k] * Q[k]
FD_q[k] <- (D_q[k])^2 * Q[k]
}
}
}
if (fdis) {
if (div_by_sp == TRUE) {
return(rbind(Q, FDis, D_q/SR, MD_q/SR, FD_q/choose(SR, 2)))
} else {
return(rbind(Q, FDis, D_q, MD_q, FD_q))
}
} else {
if (div_by_sp == TRUE) {
return(rbind(Q, D_q/SR, MD_q/SR, FD_q/choose(SR, 2)))
} else {
return(rbind(Q, D_q, MD_q, FD_q))
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.