Diversity index: Difference between revisions

Revision as of 12:26, 8 January 2012

[show] This page is a knowledge crystal of subtype method. The page identifier is Op_en5455
Moderator:Jouni (see all)
This page is a stub. You may improve it into a full page.
Upload data {{#opasnet_base_link:Op_en5455}}

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

library(OpasnetBaseUtils)
library(xtable)
data <- opasnet.csv("5/54/Tigre-data.csv", wiki = "opasnet_en", sep = ",")
data <- data[data$Species != 1 & data$Individuals != 0, ]
head(data)
species <- data$Species
amount <- data$Individuals
transect <- data$Transect

true <- function(pi, q){
pi <- pi/sum(pi)
qD <- rep(NA, length(q))
for(i in 1:length(q)){
if(q[i] == 1)
     {qD[i] <- exp(-sum(pi * log(pi)))} else
     if( q[i] == 999){
           qD[i] <- 1/max(pi)} else
     {qD[i] <- (sum(pi^q[i]))^(1 / (1 - q[i]))}}
return(qD)
}

truelist <- function(pi){
return(true(pi, c(0, 1, 2, 999)))
}

truelist(c(0.2, 0.2, 0.1, 0.4, 0.1))

true2 <- function(divj, wj){
q <- c(rep(0, length(wj)), rep(1, length(wj)), rep(2, length(wj)), rep(999, length(wj)))
wj <- rep(wj, times = 4)
     qDa <- ifelse(q == 1, exp(sum(wj * log(divj))), (sum(wj * divj^(1-q)))^1/(1 - q))
minqDj <- tapply(qDa, q, min)
qDa <- tapply(qDa, q, sum)
qDa <- qDa^(1-c(0, 1, 2, 3))
qDa[4] <- minqDj[4]
return(qDa)
}

# diversity <- function(amount = rep(1,length(species)), species = 1:length(amount), transect = 1){

     pij <- as.data.frame(as.table(tapply(amount, data.frame(Transect = transect, Species = species), sum)))
     pij[is.na(pij)] <- 0
head(pij)
colnames(pij)[3] <- "pij"

     pij$pij <- pij$pij/sum(pij$pij)
pij <- pij[pij$pij != 0, ]

head(pij)
wj <- as.data.frame(as.table(tapply(pij$pij, pij$Transect, sum)))
head(wj)
colnames(wj) <- c("Transect", "wj")
out <- merge(pij, wj)
sum(out$pij)
sum(tapply(out$wj, out$Transect, max))
head(out)
out <- tapply(out$pij, out$Transect, truelist)

# truelist WORKS NOW!

head(out)
qDj <- data.frame(Zero = 1:nrow(wj), One = NA, Two = NA, Infi = NA)
for(i in 1:nrow(wj)){
qDj[i,] <- out[[i]][1:4]
}
qDj <- c(qDj$Zero, qDj$One, qDj$Two, qDj$Infi) 

true2(qDj, wj$wj)

out <- data.frame(Name = c("True diversity with q=0", "True diversity with q=1", 
     "True diversity with q=2", "True diversity with q=∞", "Richness", "Shannon index", 
     "Simpson index", "Inverse Simpson index", "Gini-Simpson index", "Berger-Parker index"), 
     Symbol = c("0D", "1D", "2D", "∞D", "S", "H' or log(1D)", "λ or 1/(2D)", "1/λ or 2D", 
     "1-λ or 1-1/(2D)", "1/(∞D)"),
     Value = c(true(pi,0), true(pi,1), true(pi,2), 1/max(pi), length(pi), log(true(pi,1)), 
     1/true(pi,2), true(pi,2), 1-1/true(pi,2), max(pi)))
#return(out)
#}

print(xtable(out), type = 'html')

Example 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
if(individual==TRUE) out <- diversity(species = data) else out <- diversity(amount = data)
print(xtable(out), type = 'html')

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.

References

Related files