+ Show code- Hide code
####NOTE! q.wiki MUST BE DEFINED OR ALPHA DIVERSITIES WILL BE CALCULATED WRONG
####### qD calculates true i.e. gamma diversity for abundance pi with exponent q.
qD <- function(pi, q = q.wiki){ # q.wiki is the q given by the user.
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)
}
###### qDa calculates alpha diversity using gamma diversities of each transect.
qDa <- function(divj, wj = 1, q = q.wiki){
q. <- rep(q, times = length(wj))
wj <- rep(wj, each = length(q))
minqDj <- tapply(divj, q., min) # Should this be the inverse?
out <- ifelse(q. == 1, wj * log(divj), wj * divj^(1-q.))
out <- tapply(out, q., sum)
out <- ifelse(q == 1, exp(out), out^(1 / (1 - q)))
out[q == 999] <- minqDj[q == 999]
return(out)
}
####### 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 = q.wiki){
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, qD)
qDj <- NULL
for(i in 1:nrow(wj)){
qDj <- c(qDj, out[[i]])
}
gamma <- qD(pi, q)
alpha <- qDa(qDj, wj$wj, q)
dim(qDj) <- c(length(q), length(qDj)/length(q))
return(list(gamma = gamma, alpha = alpha, tr.gammas = t(qDj)))
}
########## diversity.table makes a data table out of the results.
diversity.table <- function(amount = rep(1,length(species)), species = 1:length(amount), transect = 1, q = q.wiki){
temp <- diversity(amount, species, transect, q)
gamma <- temp[["gamma"]]
alpha <- temp[["alpha"]]
beta. <- gamma / alpha
beta.A <- gamma - alpha
beta.W <- gamma / alpha - 1
beta.P <- 1 - alpha / gamma
qt <- gsub("999", "∞", q)
outlabel <- c(
paste("Gamma diversity with q=", qt, sep=""),
paste("Alpha diversity with q=", qt, sep=""),
paste("Beta diversity (strict) with q=", qt, sep=""),
paste("Beta diversity (absolute turnover) with q=", qt, sep=""),
paste("Beta diversity (Whittaker's turnover) with q=", qt, sep=""),
paste("Beta diversity (proportional turnover) with q=", qt, sep=""))
outsymbol <- c(
paste(qt, "D", sep=""),
paste(qt, "Da", sep=""),
paste(qt, "Beta.", sep=""),
paste(qt, "Beta.A", sep=""),
paste(qt, "Beta.W", sep=""),
paste(qt, "Beta.P", 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(gamma, alpha, beta., beta.A, beta.W, beta.P, length(pi), log(qD(pi,1)),
1/qD(pi,2), qD(pi,2), 1-1/qD(pi,2), 1/qD(pi,999)))
return(out)
}
cat("Functions loaded.\n")
| |
