How do I colour lines separately in rarecurve (vegan package)

Question

I'm using rarecurve (vegan) to produce rarefaction curves for nine samples, but I want them to be coloured in groups of three.

The parameters for rarecurve are:

rarecurve(x, step = 1, sample, xlab = "Sample Size", ylab = "Species", label = TRUE, ...)

With the ... passing arguments to 'plot'. However, when I replace the ellipsis with col=c(rep("blue",3), rep("red",3), rep("darkgreen",3)), all lines appear as blue. How do I colour the lines individually?

It takes almost three hours to compute each graph, so trial and error testing is a bit laborious!

Answer 1

## example from ?vegan::rarecurve
library(vegan)
data(BCI)
S <- specnumber(BCI)
(raremax <- min(rowSums(BCI)))
Srare <- rarefy(BCI, raremax)
plot(S, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
abline(0, 1)
rarecurve(BCI, step = 20, sample = raremax, col = "blue", cex = 0.6)

# using new function
plot(S, Srare, xlab = "Observed No. of Species", ylab = "Rarefied No. of Species")
abline(0, 1)
rarec(BCI, step = 20, sample = raremax, cex = 0.6)

The problem is in these lines in vegan::rarecurve

for (ln in seq_len(length(out))) {
  N <- attr(out[[ln]], "Subsample")
  lines(N, out[[ln]], ...)

where each line is made individually by lines which in turn is only taking the first color it sees in the color argument passed by ..., which is blue in your case. After applying a simple hack to this loop:

for (ln in seq_len(length(out))) {
  N <- attr(out[[ln]], "Subsample")
  lines(N, out[[ln]], col = cols[ln], ...)

and specifying a new argument, cols, in the rarecurve function rather than passing col onto plot and lines:

cols = c(rep('red', nrow(x) / 2), rep('blue', nrow(x) / 2))

Here is the new function

rarec <- function (x, step = 1, sample, xlab = "Sample Size", ylab = "Species", 
          label = TRUE, cols = c(rep('red', nrow(x) / 2), rep('blue', nrow(x) / 2)), ...) {
  tot <- rowSums(x)
  S <- specnumber(x)
  nr <- nrow(x)
  out <- lapply(seq_len(nr), function(i) {
    n <- seq(1, tot[i], by = step)
    if (n[length(n)] != tot[i]) 
      n <- c(n, tot[i])
    drop(rarefy(x[i, ], n))
  })
  Nmax <- sapply(out, function(x) max(attr(x, "Subsample")))
  Smax <- sapply(out, max)
  plot(c(1, max(Nmax)), c(1, max(Smax)), xlab = xlab, ylab = ylab, 
       type = "n", ...)
  if (!missing(sample)) {
    abline(v = sample)
    rare <- sapply(out, function(z) approx(x = attr(z, "Subsample"), 
                                           y = z, xout = sample, rule = 1)$y)
    abline(h = rare, lwd = 0.5)
  }
  for (ln in seq_len(length(out))) {
    N <- attr(out[[ln]], "Subsample")
    lines(N, out[[ln]], col = cols[ln], ...)
  }
  if (label) {
    ordilabel(cbind(tot, S), labels = rownames(x), ...)
  }
  invisible(out)
}