Question
How to calculate diversity indices?
Answer
Upload your data to Opasnet Base. Use the function diversity to calculate the most common indices.
Actual function diversity
+ Show code- Hide code
####### qD calculates true diversity for abundance pi with exponent q.
qD <- function(pi, q = q.div){ # q.div is the q used by the diversity function.
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.div){
q. <- rep(q, each = length(wj))
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.div: exponent for the diversity calculation
diversity <- function(amount = rep(1,length(species)), species = 1:length(amount), transect = 1, q.div = 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, qD)
qDj <- rep(NA, nrow(wj)*length(q.div))
for(i in 1:nrow(wj)){
for(j in 1:length(q.div)){
qDj[i + (j - 1) * nrow(wj)] <- out[[i]][j]
}
}
out <- qDa(qDj, wj$wj, q.div)
out2 <- qD(pi, q.div)
qt <- gsub("999", "∞", q.div)
outlabel <- c(paste("Alpha diversity with q=", qt, sep=""), paste("Gamma diversity with q=", qt, sep=""))
outsymbol <- c(paste(qt, "Da", sep=""), paste(qt, "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)
}
cat("Functions loaded.\n")
| |
Examples
Example 1 to use function
+ Show code- Hide code
library(xtable)
if(is.null(data)) {data <- c(1,4,6,9,3,4,5,6,5,4,3,3,5,5,7,5,5,4,3,4,5,6,8,9,5,4,5,4,3,3,4,9,6,6,4,5,3,2,1,1,2,3,4,3,2,3,4,5,6,7,7)}
data
diversity(data)
if(individual==TRUE) out <- diversity(species = data, q = 1) else out <- diversity(amount = data, q = 1)
print(xtable(out), type = 'html')
| |
Example 2
TODO: {{#todo:Mikä ihme tässä koodissa on vikana? q.diviä ei muka löydy mutten ymmärrä miksei. Ongelma on diversity-funktiossa, joka on määritelty sivun ensimmäisessä rcodessa. --
Jouni 00:36, 9 January 2012 (EET)|Teemu Rintala|}}
The data should be given in R format as a list of values in parenthesis, beginning with c:
c(3,5,3,5,2,1,3,3,4,2) or equivalently c(0.1,0.2,0.4,0.1,0.2)
where the values are either
- identifiers of the species 1,2,3... in which the individuals belong (one entry per individual), or
- abundancies of species, i.e. proportions of individuals belonging to each species among the whole population (one entry per species).
Rationale
Diversity indices are thoroughly described in Wikipedia.
See also
References
Related files
<mfanonymousfilelist></mfanonymousfilelist>
