Difference between revisions of "Diversity index"

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)))
}

true2 <- function(divj, wj){
q <- rep(c(0, 1, 2, 999), each = length(wj))
wj <- rep(wj, times = 4)
#print(divj)
#print(wj)
#print(q)
     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)

#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
qDj <- c(qDj$Zero, qDj$One, qDj$Two, qDj$Infi) 
#qDj
#wj
qDa <- true2(qDj, wj$wj)

out <- data.frame(Name = c("Alpha diversity with q=0", "Alpha diversity with q=1", "Alpha diversity with q=2", "Alpha diversity with q=∞", "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("0Da", "1Da", "2Da", "∞Da", "0D", "1D", "2D", "∞D", "S", "H' or log(1D)", "λ or 1/(2D)", "1/λ or 2D", 
     "1-λ or 1-1/(2D)", "1/(∞D)"),
     Value = c(qDa, 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')

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')