User Tools

Site Tools


en:div-ind_rscript
# Three communities - even, moderately uneven and highly uneven
 
image.bubble <- function (mat, col = NULL, cex = 1, jitter.factor = 1, main = 'Community')
{
  if (is.null (col)) col <- c('#ffffff', heat.colors (max (mat)))
  long <- as.data.frame.table (mat)
  opar <- par ('mar')
  par (mar = c(2,2,4,2))
  plot (jitter (as.numeric (long$Var1), factor = jitter.factor) ~ jitter (as.numeric (long$Var2), factor = jitter.factor), col = col[long$Freq], pch = 16, cex = cex, axes = F, asp = 1, main = main, xlab = '', ylab = '')
  par (mar = opar)
}
 
 
prob.even <- rep (1/12, 12)
prob.mod <- dgeom (1:12, 0.35)
prob.mod <- prob.mod/sum (prob.mod)
prob.high <- dgeom (1:12, 0.7)
prob.high <- prob.high/sum (prob.high)
 
mat.even <- replicate (20, sample (1:12, 20, prob = rep (1, 12), replace = T))
mat.mod <- replicate (20, sample (1:12, 20, prob = prob.mod, replace = T))
mat.high <- replicate (20, sample (1:12, 20, prob = prob.high, replace = T))
set.seed (2344)
cols <- sample (RColorBrewer:::brewer.pal(n = 12, name = 'Paired'))
 
 
# Plotting even, moderately uneven and highly uneven distribution
 
png ('bubbles-sad-even.png', width = 8, height = 4, units = 'in', res = 600)
par (mfrow = c(1,2))
mat <- mat.even
prob <- prob.even
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community A\n (perfectly even)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .8)
legend ('topright', bty = 'n', legend = c('Species richness = 12', 'Shannon entropy = 2.48', 'Gini-Simpson index = 0.92'))
dev.off ()
 
png ('bubbles-sad-mod_uneven.png', width = 8, height = 4, units = 'in', res = 600)
par (mfrow = c(1,2))
mat <- mat.mod
prob <- prob.mod
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community B\n (moderately uneven)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .8)
legend ('topright', bty = 'n', legend = c('Species richness = 12.00', 'Shannon entropy = 1.81', 'Gini-Simpson index = 0.79'))
dev.off ()
 
png ('bubbles-sad-high_uneven.png', width = 8, height = 4, units = 'in', res = 600)
par (mfrow = c(1,2))
mat <- mat.high
prob <- prob.high
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community C\n (highly uneven)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .8)
legend ('topright', bty = 'n', legend = c('Species richness = 12.00', 'Shannon entropy = 0.87', 'Gini-Simpson index = 0.46'))
dev.off ()
 
# three at one
png ('bubbles-sad-3in1.png', width = 9, height = 6, units = 'in', res = 600)
par (mfcol = c(2,3))
mat <- mat.even
prob <- prob.even
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community A\n (perfectly even)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .7)
legend ('topright', bty = 'n', legend = c('Species richness = 12', 'Shannon entropy = 2.48', 'Gini-Simpson index = 0.92'))
 
mat <- mat.mod
prob <- prob.mod
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community B\n (moderately uneven)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .7)
legend ('topright', bty = 'n', legend = c('Species richness = 12.00', 'Shannon entropy = 1.81', 'Gini-Simpson index = 0.79'))
 
mat <- mat.high
prob <- prob.high
image.bubble (mat, col = cols, cex = 2, jit = 2, main = 'Community C\n (highly uneven)')
barplot (prob*100, col = cols, ylim = c(0, 70), main = 'Species abundance distribution', las = 1, ylab = 'Number of individuals', xlab = 'Species ID', names.arg = 1:12, cex.names = .7)
legend ('topright', bty = 'n', legend = c('Species richness = 12.00', 'Shannon entropy = 0.87', 'Gini-Simpson index = 0.46'))
dev.off ()
 
 
diverprof.prob <- function (prob, q = seq (0, 5, by = 0.1))
{
  abund <- table (sample (1:length (prob), 1000, replace = T))
  cbind (q, unlist (lapply (seq (0,5, by = 0.1), FUN = function (q) vegetarian:::d(t(as.matrix (abund)), q = q))))
}
 
diverprof <- function (mat, q = seq (0, 5, by = 0.1))
  cbind (q, unlist (lapply (q, FUN = function (q) vegetarian:::d (mat, q = q))))
 
set.seed (123)
mat.even <- rep (100000, 12)
prof.even <- diverprof (mat.even)
 
mat.mod <- as.vector (table (sample (1:12, 1200000, prob = prob.mod, replace = T)))
prof.mod <- diverprof (mat.mod)
 
mat.high <- as.vector (table (sample (1:12, 1200000, prob = prob.high, replace = T)))
prof.high <- diverprof (mat.high)
 
 
png ('diversity-profile-even-mod-high.png', height = 5, width = 8, res = 600, units = 'in')
par (mar = c(5.1, 4.1, 4.1, 20))
plot (prof.even, ylim = c (0,15), type = 'l', las = 1, lwd = 2, ylab = 'Effective number of species', main = 'Diversity profiles for three communities', xlab = expression (italic (q)))
lines (prof.mod, lty = 'dashed', col = 'red', lwd = 2)
lines (prof.high, lty = 'dotted', col = 'blue', lwd = 2)
points (x = c(0), y = c(12), cex = 1.7, pch = (c(22)), bg = 'white')
points (x = c(1,1,1), y = c(12, prof.mod [prof.mod[,'q'] == 1, 2], prof.high [prof.high[,'q'] == 1, 2]), cex = 1.7, pch = c(21), bg = 'white')
points (x = c(2,2,2), y = c(12, prof.mod [prof.mod[,'q'] == 2, 2], prof.high [prof.high[,'q'] == 2, 2]), cex = 1.7, pch = c(24), bg = 'white')
par (xpd = T)
legend (6,15, title = 'Species abundance distribution', legend = c('even', 'moderately uneven', 'highly uneven'), bty = 'n', lwd = 2, lty = c('solid', 'dashed', 'dotted'), col = c('black', 'red', 'blue'), title.adj = 0)
legend (6, 10, title = 'Diversity index', pch = c(22,21,24), pt.cex = 1.7, legend = c('species richness', 'Shannon diversity', 'Simpson diversity'), bty = 'n', title.adj = 0)
dev.off ()
# Relationship of sp.richness, Shannon and Simpson on number of species and (un)evenness
 
library (vegetarian)
library (vegan)
n <- 1:100
p <- seq (0.000001, .7, length = 100)
np <- expand.grid (n = n, p = p)
 
q012 <- apply (np, 1, FUN = function (x) 
  {
  n <- x[1]
  p <- x[2]
  prob <- (dgeom (1:n, p))
  prob <- prob/sum (prob)
  res <- c(q0 = d (prob, q = 0), q1 = d (prob, q = 1), q2 = d (prob, q = 2))
})
 
 
contour (x = p, y = n, matrix (q012[1,], ncol = 100, byrow = T), method = 'edge', labcex = 1)
contour (x = p, y = n, matrix (q012[2,], ncol = 100, byrow = T), add = T, col = 'green', drawlabels = F)
contour (x = p, y = n, matrix (q012[3,], ncol = 100, byrow = T), add = T, col = 'red', method = 'edge')
 
opar <- par ()
png ('rich-shan-simp-effect-n-sp.png', width = 9, height = 3, units = 'in', res = 600)
par (mfrow = c(1,3))
contour (x = p, y = n, matrix (q012[1,], ncol = 100, byrow = T), method = 'edge', labcex = 1, las = 1, levels = seq (0, 100, 1), log = 'y', drawlabels = F, col = 'grey90', main = 'Species richness', xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (q012[1,], ncol = 100, byrow = T), method = 'edge', labcex = 0.7, las = 1, levels = c(1,2,3,4,5,10, 20, 50, 100), log = 'y', add = T)
 
contour (x = p, y = n, matrix (q012[2,], ncol = 100, byrow = T), method = 'edge', labcex = 1, las = 1, levels = seq (0, 100, 1), log = 'y', drawlabels = F, col = 'grey90', main = "Shannon effective number of species (exp H)", xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (q012[2,], ncol = 100, byrow = T), method = 'edge', labcex = 0.7, levels = c(1,2,3,4,5,10, 20, 50, 100), log = 'y', add = T)
 
contour (x = p, y = n, matrix (q012[3,], ncol = 100, byrow = T), method = 'edge', labcex = 1, las = 1, levels = seq (0, 100, 1), log = 'y', drawlabels = F, col = 'grey90', main = "Simpson's effective number of species (1/D)", xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (q012[3,], ncol = 100, byrow = T), method = 'edge', labcex = 0.7, levels = c(1,2,3,4,5,10, 20, 50, 100), log = 'y', add = T)
dev.off ()
 
 
 
png ('rich-shan-simp.png', width = 9, height = 3, units = 'in', res = 600)
par (mfrow = c(1,3))
contour (x = p, y = n, matrix (q012[1,], ncol = 100, byrow = T), method = 'edge', labcex = 1, las = 1, log = 'y', main = 'Species richness', nlevels = 50, drawlabels = F, col = 'grey90', xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (q012[1,], ncol = 100, byrow = T), method = 'edge', labcex = .7, add = T)
 
contour (x = p, y = n, matrix (log (q012[2,]), ncol = 100, byrow = T), drawlabels = F, las = 1, log = 'y', main = 'Shannon entropy', nlevels = 50, col = 'grey90', method = 'edge', labcex = 1, xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (log (q012[2,]), ncol = 100, byrow = T), drawlabels = T, add = T, method = 'edge', labcex = .7)
 
contour (x = p, y = n, matrix (1-1/q012[3,], ncol = 100, byrow = T), method = 'edge', labcex = 1, las = 1, log = 'y', main = 'Gini-Simpson index (1-Simpson)', drawlabels = F, col = 'grey90', nlevels = 50, xlab = list ('Uneveness', cex = 1.5), ylab = list ('Number of species', cex = 1.5))
points (x = c(0, 0.35, 0.7), y = rep (12, 3), cex = 3, pch = 21, bg = 'white', col = 'darkgrey')
text (x = c(0, 0.35, 0.7), y = rep (12, 3), labels = c('A', 'B', 'C'), col = 'darkgrey')
contour (x = p, y = n, matrix (1-1/q012[3,], ncol = 100, byrow = T), method = 'edge', labcex = .7, add = T)
dev.off ()
 
# Species richness, Shannon and Simpson for perfectly even communities
div <- t(sapply (seq (1, 100), FUN = function (sp)
{
  H <- log (sp)
  D <- 1 - 1/sp
  c(sp = sp, H=H, D=D)
}))
 
png ('Shannon-and-Simpson-on-sp-richness.png', height = 6, width = 6, res = 600, units = 'in')
par (mar = c(5,5,4,5))
plot (H ~ sp, div, type = 'l', las = 1, ann = F, axes = F, lwd = 2)
axis (1)
axis (2, las = 1)
box (bty = 'l')
plot.window (xlim = c(1, 100), ylim = c(0, 1))
lines (D ~ sp, div, col = 'red', lwd = 2, lty = 'dashed')
axis (4, las = 1, col = 'red', lwd = 2, lty = 'dashed')
mtext ('Gini-Simpson index (D = 1-(1/S))', side = 4, line = 3)
title (xlab = 'Species richness (S)', ylab = 'Shannon index (H = log S)', main = 'Dependence of Shannon and Gini-Simpson index\n on number of species in perfectly even community')
legend ('bottomright', legend = c('Shannon entropy index', 'Gini-Simpson index'), lty = c('solid', 'dashed'), lwd = 2, col = c('black', 'red'), bty = 'n')
dev.off ()
en/div-ind_rscript.txt · Last modified: 2017/10/11 20:36 (external edit)