Difference between revisions of "Diversity index"

####NOTE! q.wiki MUST BE DEFINED OR ALPHA DIVERSITIES WILL BE CALCULATED WRONG

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

qD <- function(pi, q = q.wiki){ # q.wiki is the q given by the user.
	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]))
			}
		}

	q[q == 999] <- "∞"
	names(out) <- q
	return(out)
}

###### qDa calculates alpha diversity using gamma diversities of each transect.

qDa <- function(divj, wj = 1, q = q.wiki){
	q. <- rep(q, times = length(wj))
	wj <- rep(wj, each = length(q))

	minqDj <- tapply(divj, q., min) # Should this be the inverse?
	out <- ifelse(q. == 1, wj * log(divj), wj * divj^(1-q.))
	out <- tapply(out, q., sum)
	out <- ifelse(q == 1, exp(out), out^(1 / (1 - q)))
	out[q == 999] <- minqDj[q == 999]
	q[q == 999] <- "∞"
	names(out) <- q
	return(out)
}

####### 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 = q.wiki){

	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, qD)

	qDj <- NULL

	for(i in 1:nrow(wj)){
		qDj <- c(qDj, out[[i]])
	}

	gamma <- qD(pi, q)
	alpha <- qDa(qDj, wj$wj, q)
	dim(qDj) <- c(length(q), length(qDj)/length(q))
	q[q == 999] <- "∞"
	dimnames(qDj) <- list(q, levels(as.factor(pij$Transect)))

	return(list(gamma = gamma, alpha = alpha, tr.gammas = t(qDj)))
}

########## diversity.table makes a data table out of the results.

diversity.table <- function(amount = rep(1,length(species)), species = 1:length(amount), transect = 1, q = q.wiki){

	temp <- diversity(amount, species, transect, q)
	gamma <- temp[["gamma"]]
	alpha <- temp[["alpha"]]
	beta. <- gamma / alpha
	beta.A <- gamma - alpha
	beta.W <- gamma / alpha - 1
	beta.P <- 1 - alpha / gamma
	qt <- gsub("999", "∞", q)

	outlabel <- c(
		paste("Gamma diversity with q=", qt, sep=""),
		paste("Alpha diversity with q=", qt, sep=""),
		paste("Beta diversity (strict) with q=", qt, sep=""),
		paste("Beta diversity (absolute turnover) with q=", qt, sep=""),
		paste("Beta diversity (Whittaker's turnover) with q=", qt, sep=""),
		paste("Beta diversity (proportional turnover) with q=", qt, sep=""))

	outsymbol <- c(
		paste(qt, "D", sep=""),
		paste(qt, "Da", sep=""),
		paste(qt, "Beta.", sep=""),
		paste(qt, "Beta.A", sep=""),
		paste(qt, "Beta.W", sep=""),
		paste(qt, "Beta.P", 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(gamma, alpha, beta., beta.A, beta.W, beta.P, 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")

library(xtable)
q.wiki = 1
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 = q.wiki) else out <- diversity(amount = data, q = q.wiki)
print(xtable(out), type = 'html')

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

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

diversity(amount, species, transect, q.wiki)
print(xtable(diversity(amount, species, transect, q.wiki)[[3]]), type = 'html')
print(xtable(diversity.table(amount, species, transect, q.wiki)), type = 'html')