Difference between revisions of "Diversity index"

####### qD calculates true diversity for abundance pi with exponent q.

qD <- function(pi, q){
print(pi)
print(q)
print(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)
}
library(OpasnetBaseUtils)
data <- c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20)
diversity(amount=data, species = 1:20, transect = c(rep(1, 10), rep(2, 10)), q = c(0, 1, 2))

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

library(OpasnetBaseUtils)
library(xtable)
data <- opasnet.csv(data, wiki = "opasnet_en", sep = ",")
data <- data[data$Species != 1, ]

species <- data$Species
amount <- data$Individuals
transect <- data$Transect

out <- diversity(amount, species, transect, q)
print(xtable(out), type = 'html')