+ Show code- Hide code
####### qD calculates true diversity for abundance pi with exponent q.
qD <- function(pi, q){
pi <- pi/sum(pi)
out <- rep(NA, length(q))
for(i in 1:length(q)){
if(q[i] == 1)
{out[i] <- exp(-sum(pi * log(pi)))} else
if( q[i] == 999){
out[i] <- 1/max(pi)} else
{out[i] <- (sum(pi^q[i]))^(1 / (1 - q[i]))
}
}
return(out)
}
##### qD4tapply is the same as the qD but using only a single parameter, because that is what is required for tapply.
qD4tapply <- function(pi){
return(qD(pi, q))
}
###### qDa calculates alpha diversity using gamma diversities of each transect.
qDa <- function(divj, wj = 1, q){
q. <- rep(q, each = length(wj))
# wj <- rep(wj, times = length(q))
out <- ifelse(q. == 1, exp(sum(wj * log(divj))), (sum(wj * divj^(1-q.)))^(1/(1 - q.)))
minqDj <- tapply(out, q., min)
out <- tapply(out, q., sum)^(1-q)
out[q == 999] <- minqDj[q == 999]
return(out)
}
## There seems to be an error in qDa. What is it?
####### diversity is the function for the user. It calculates several diversity
####### indices. Parameters:
## amount: the number of individuals of a species in a transect.
## species: an identifier for a species
## transect: an identifier for a transect
## q: exponent for the diversity calculation
diversity <- function(amount = rep(1,length(species)), species = 1:length(amount), transect = 1, q = 0){
pij <- as.data.frame(as.table(tapply(amount, data.frame(Transect = transect, Species = species), sum)))
colnames(pij)[3] <- "pij"
pij <- dropall(pij[!is.na(pij$pij) & pij$pij != 0, ])
pij$pij <- pij$pij/sum(pij$pij)
wj <- as.data.frame(as.table(tapply(pij$pij, pij$Transect, sum)))
pi <- as.data.frame(as.table(tapply(pij$pij, pij$Species, sum)))[, 2]
colnames(wj) <- c("Transect", "wj")
out <- merge(pij, wj)
out <- tapply(out$pij, out$Transect, qD4tapply)
qDj <- rep(NA, nrow(wj)*length(q))
for(i in 1:nrow(wj)){
for(j in 1:length(q)){
qDj[i + (j - 1) * nrow(wj)] <- out[[i]][j]
}
}
out <- qDa(qDj, wj$wj, q)
out2 <- qD(pi, q)
q <- gsub("999", "∞", q)
outlabel <- c(paste("Alpha diversity with q=", q, sep=""), paste("Gamma diversity with q=", q, sep=""))
outsymbol <- c(paste(q, "Da", sep=""), paste(q, "D", sep=""))
out <- data.frame(Name = c(outlabel, "Richness", "Shannon index",
"Simpson index", "Inverse Simpson index", "Gini-Simpson index", "Berger-Parker index"),
Symbol = c(outsymbol, "S", "H' or log(1D)", "λ or 1/(2D)", "1/λ or 2D",
"1-λ or 1-1/(2D)", "1/(∞D)"),
Value = c(out, out2, length(pi), log(qD(pi,1)),
1/qD(pi,2), qD(pi,2), 1-1/qD(pi,2), 1/qD(pi,999)))
return(out)
}
| |