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
####NOTE! q.wiki MUST BE DEFINED OR ALPHA DIVERSITIES WILL BE CALCULATED WRONG
####### qD calculates true diversity for abundance pi with exponent q.
qD <- function(pi, q = q.wiki){ # q.wiki is the q given by the user.
q. <- rep(q, times = length(pi))
pi <- rep(pi, each = length(q))
pi <- pi/sum(pi)
# out <- ifelse(q. == 1,
# exp(-sum(pi * log(pi))),
# ifelse(q. == 999,
# 1/max(pi),
# tapply(pi^q., q., sum)^(1 / (1 - q))))
#
maxqD <- 1/tapply(pi, q., max)
out <- ifelse(q. == 1, pi * log(pi), pi^q.)
out <- tapply(out, q., sum)
out <- ifelse(q == 1, exp(-out), out^(1 / (1 - q)))
out[q == 999] <- maxqD[q == 999]
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))
print(divj)
print(wj)
print(q.)
print(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)
print(out)
qDj <- NULL
for(i in 1:nrow(wj)){
qDj <- c(qDj, out[[i]])
}
out <- qDa(qDj, wj$wj, q)
out2 <- qD(pi, q)
qt <- gsub("999", "∞", q)
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)
q.wiki = 1
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 = q.wiki) else out <- diversity(amount = data, q = q.wiki)
print(xtable(out), type = 'html')
| |
Example 2
--#: 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) to Teemu Rintala
- --# : Ilmeisesti funktioparametrien defaulteiksi ei saa laittaa ennalta määrittelemättömiä muuttujia. Vaihdoin defaultiksi "q = 0", "q = q.div" sijaan. --Teemu R 16:23, 9 January 2012 (EET)
- ←# : OK. Luulin, että funktion sisältä kutsuttava funktio voi saada toisen funktion sisällä olevia oletusarvoja, muttei näin välttämättä ole. Oletusarvot ovat kuitenkin tärkeitä, koska qD-funktiota käytetään tapplyssa, jossa voi antaa vain yhden parametrin ja loput tulevat oletusarvoista. Toteutin tämän niin, että käytetään q.wikiä, joka yleensä kysytään käyttäjärajapinnassa tai muuten. --Jouni 20:15, 9 January 2012 (EET)
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>