User Tools

Site Tools


en:div-ind_rscript

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

en:div-ind_rscript [2017/10/11 20:36] (current)
Line 1: Line 1:
 +<code rsplus>
 +# 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 ()
 +
 +</​code>​
 +
 +<code rsplus>
 +# 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 ()
 +</​code>​
en/div-ind_rscript.txt · Last modified: 2017/10/11 20:36 (external edit)