Three functions are included: random.drift generates the time sequence of changes in a community caused by random drift; plot.rd and barplot.rd visualise these changes. Simplified version (with fewer parameters to modify) of this function is also published as a Shiny application.

To define the functions directly from this website, use the following script:

source ('http://www.davidzeleny.net/anadat-r/doku.php/en:customized_functions:random.drift?do=export_code&codeblock=1')

• no.spe Number of species. Integer.
• no.ind Number of individuals. Integer.
• no.gen Number of generations. Integer.
• ratio Ratio between species abundances (could be relative or absolute numbers). Default is that all species have the same abundance.
• set.seed Should we fix the current permutation? Logical value.
• seed Seed for the generator of pseudorandom values.
• replicates Number of simulations to plot in one plot as panels.
• draw.plot Should the function random.drift directly plot the result? Logical value.
• x Object returned by random.drift function.
• to How many generations should be ploted? Maximum is no.gen.
• col Color palette of the species. Defaults is rainbow (no.spe). Vector of colors of the length equal to no.spe.
• sep Should be the species trends visually separated by a line? Logical value, default TRUE.
• col.sep Color of the line separating the species. Vector of the length equal to one. Default white.

random.drift.r

random.drift <- function (no.spe = 10, no.ind = 100, no.gen = 1000, ratio = rep (1/no.spe, no.spe), set.seed = FALSE, seed = 1234, replicates = 1, draw.plot = TRUE, ...)
{
  replicates <- as.numeric (replicates)
  ratio <- ratio/sum (ratio)*no.ind
 
  com.res.out <- list ()
  for (i in seq (1, replicates))
  {
    if (set.seed) set.seed (seed+i-1)
    com0 <- ordered (unlist (lapply (1:no.spe, FUN = function (n) rep (n, ceiling (ratio[n])))), levels = 1:no.spe)
    com <- sample (com0, no.ind)
    com.res <- matrix (NA, nrow = no.spe, ncol = no.gen)
    com.res[,1] <- as.vector (table (com))
    for (n in seq (2, no.gen))
    {
      die.reproduce <- sample (1:no.ind, 2)
      com [die.reproduce[1]] <- com[die.reproduce[2]]
      com.res[,n] <- as.vector (table (com))
    }
    com.res.out[[i]] <- com.res
  }
  com.res.out$pars <- list (no.spe = no.spe, no.ind = no.ind, no.gen = no.gen, ratio = ratio, replicates = replicates)
  if (draw.plot) plot.rd (com.res.out, ...)
  return (com.res.out)
}
 
plot.rd <- function (x, to = x$pars$no.gen, col = rainbow (x$pars$no.spe), sep = TRUE, col.sep = "white", max.no.bars = 100, panel.letter.add = FALSE, title.add = TRUE)
{
  if (to > x$pars$no.gen) stop ("Argument 'to' cannot be larger than the number of generations in 'x'.")
  pars <- x$pars
  if (pars$replicates > 1) par (mfrow = c(ceiling (sqrt (pars$replicates)), ceiling (sqrt (pars$replicates))))
  for (i in seq (1, pars$replicates))
  {
    com.res <- x[[i]][, 1:to]
    no.gen.temp <- ncol (com.res)
    if (no.gen.temp > max.no.bars) com.res.draw <- com.res[, ceiling (seq (1, no.gen.temp, len = max.no.bars))] else com.res.draw <- com.res
    barplot (com.res.draw, space = 0, border = NA, col = col, xlab = 'Number of generations', ylab = 'Number of individuals')
    if (sep) matplot (x = seq (0, ncol (com.res.draw)), rbind (cumsum (com.res.draw[,1]), t(apply (com.res.draw, 2, cumsum))), type = 'S', col = col.sep, lty = 'solid', add = T)
    if (title.add) title (main = list (paste (pars$no.ind, 'individuals,', no.gen.temp, 'generations\nStart: ', sum (com.res[,1]>0), 'species, end: ', sum (com.res[,no.gen.temp]>0), 'species'), font = 1))
    if (no.gen.temp > max.no.bars) {coef <- no.gen.temp/max.no.bars; pretty.tick <- pretty (1:no.gen.temp); axis (1, at = pretty.tick/coef, labels = pretty.tick )} else axis (1)
    if (panel.letter.add) legend ('topright', LETTERS[i], bty = 'n', adj = 1, text.font = 2)
  }
 
}
 
barplot.rd <- function (x, to = x$pars$no.gen, col1 = "grey", col2 = rainbow (x$pars$no.spe), first = 1, last = to, sort.by = 'none', panel.letter.add = FALSE, title.add = TRUE, alpha = 0.5)
{
  if (to > x$pars$no.gen) stop ("Argument 'to' cannot be larger than the number of generations in 'x'.")
  pars <- x$pars
  if (pars$replicates > 1) par (mfrow = c(ceiling (sqrt (pars$replicates)), ceiling (sqrt (pars$replicates))))
  for (i in seq (1, pars$replicates))
  {
    com.res <- x[[i]][, 1:to]
    no.gen.temp <- ncol (com.res)
    if (sort.by == 'first') sorting <- order (com.res[,first], decreasing = T)
    if (sort.by == 'last') sorting <- order (com.res[,last], decreasing = T)
    if (sort.by == 'none') sorting <- 1:nrow (com.res)
    barplot (com.res[sorting,first], col = col1, ylim = c(0, max (com.res[, c(first, last)])), xlab = 'Species', ylab = 'No individuals')
    barplot (com.res[sorting,last], add = T, col = scales::alpha (col2[sorting], alpha = alpha))
    if (title.add) title (main = list (paste (pars$no.ind, 'individuals,', no.gen.temp, 'generations\nStart: ', sum (com.res[,1]>0), 'species, end: ', sum (com.res[,no.gen.temp]>0), 'species'), font = 1))
    if (panel.letter.add) legend ('topright', LETTERS[i], bty = 'n', adj = 1, text.font = 2)
  }
}