|
|
| Line 18: |
Line 18: |
| | | | |
| | qD <- function(pi, q){ | | qD <- function(pi, q){ |
| − | print(pi)
| |
| − | print(q)
| |
| − | print(q)
| |
| | pi <- pi/sum(pi) | | pi <- pi/sum(pi) |
| | out <- rep(NA, length(q)) | | out <- rep(NA, length(q)) |
| Line 99: |
Line 96: |
| | return(out) | | return(out) |
| | } | | } |
| − | library(OpasnetBaseUtils)
| + | |
| − | data <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
| + | cat("Functions loaded.\n") |
| − | diversity(amount=data, species = 1:20, transect = c(rep(1, 10), rep(2, 10)), q = c(0, 1, 2))
| |
| | | | |
| | </rcode> | | </rcode> |
| Line 121: |
Line 117: |
| | '''Example 2 | | '''Example 2 |
| | | | |
| − | <rcode include="page:Diversity_index|name:answer|page:OpasnetBaseUtils|name:generic" variables=" | + | <rcode include="page:OpasnetBaseUtils|name:generic|page:Diversity_index|name:answer" variables=" |
| − | name:data|description:Select your data|type:selection|options:'5/54/Tigre-data.csv';Tigre|default:'5/54/Tigre-data.csv'| | + | name:wikidata|description:Select your data|type:selection|options:'5/54/Tigre-data.csv';Tigre|default:'5/54/Tigre-data.csv'| |
| − | name:q|description:Which q values you want to calculate.|type:checkbox|options:0;0;1;1;2;2;3;3;999;∞|default:0;1;2;999 | + | name:wikiq|description:Which q values you want to calculate.|type:checkbox|options:0;0;1;1;2;2;3;3;999;∞|default:0;1;2;999 |
| | "> | | "> |
| | library(OpasnetBaseUtils) | | library(OpasnetBaseUtils) |
| | library(xtable) | | library(xtable) |
| − | data <- opasnet.csv(data, wiki = "opasnet_en", sep = ",") | + | data <- opasnet.csv(wikidata, wiki = "opasnet_en", sep = ",") |
| | data <- data[data$Species != 1, ] | | data <- data[data$Species != 1, ] |
| | | | |
| Line 134: |
Line 130: |
| | transect <- data$Transect | | transect <- data$Transect |
| | | | |
| − | out <- diversity(amount, species, transect, q) | + | out <- diversity(amount, species, transect, wikiq) |
| | print(xtable(out), type = 'html') | | print(xtable(out), type = 'html') |
| | | | |
| | + | diversity(amount=1:20, species = 1:20, transect = c(rep(1, 10), rep(2, 10)), q = c(0, 1, 2)) |
| | | | |
| | </rcode> | | </rcode> |
Revision as of 19:29, 8 January 2012
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){
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)
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")
| |
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
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>
