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