Sampling effort/effectiveness rarefaction curve for set number of quadrats at two sites

Question

Goal: produce a plot of rarefaction curves for two sites, showing how species diversity increases with the number of quadrats sampled, using data that is currently in a dataframe.

Problem: A lot of tutorials take data that has already been transformed into a list or matrix, and it is unclear how the authors got there from data which may be collected as a dataframe. Additionally, many rarefaction curves have number of individuals along the x-axis, whereas I want the number of quadrats sampled at the site. It is unclear to me how to go from a dataframe like this, to the format required for rarefaction in most packages, and then how to select the number of sampling units as the x-axis.

Here are some example data:

['''

dat_rare =
structure(list(SITE = c("FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", 
"FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", 
"FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", 
"FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", 
"PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", 
"PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", 
"PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", 
"PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS"), QUADRAT = c("Q1", 
"Q1", "Q1", "Q1", "Q1", "Q10", "Q10", "Q10", "Q10", "Q11", "Q11", 
"Q11", "Q11", "Q12", "Q12", "Q6", "Q6", "Q7", "Q7", "Q7", "Q7", 
"Q8", "Q8", "Q8", "Q8", "Q8", "Q9", "Q9", "Q9", "Q9"), SPECIES = c("PERONS_LIMPET", 
"PURPLE_BARNACLE", "SERPENT_SKIN_CHITON", "STAR_LIMPET", "SURF_BARNACLE", 
"DARWINS_BARNACLE", "PERONS_LIMPET", "PIED_LIMPET", "SURF_BARNACLE", 
"PERONS_LIMPET", "PIED_LIMPET", "SERPENT_SKIN_CHITON", "SURF_BARNACLE", 
"GREEN_CHITON", "PERONS_LIMPET", "PERONS_LIMPET", "TURBAN_SHELL", 
"BLACK_LIMPET", "PERONS_LIMPET", "PIED_LIMPET", "SURF_BARNACLE", 
"BLACK_LIMPET", "PERONS_LIMPET", "PIED_LIMPET", "RIBBED_LIMPET", 
"SURF_BARNACLE", "BLACK_LIMPET", "PERONS_LIMPET", "PIED_LIMPET", 
"SERPENT_SKIN_CHITON"), n = c(836L, 3L, 1L, 108L, 4L, 4L, 60L, 
34L, 24L, 36L, 11L, 2L, 228L, 8L, 20L, 25L, 26L, 20L, 15L, 300L, 
6L, 15L, 60L, 520L, 1L, 3L, 385L, 420L, 210L, 4L)), class = c("grouped_df", 
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -30L), groups = structure(list(
    SITE = c("FOSSIL_CLIFFS", "FOSSIL_CLIFFS", "FOSSIL_CLIFFS", 
    "FOSSIL_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", "PAINTED_CLIFFS", 
    "PAINTED_CLIFFS"), QUADRAT = c("Q1", "Q10", "Q11", "Q12", 
    "Q6", "Q7", "Q8", "Q9"), .rows = structure(list(1:5, 6:9, 
        10:13, 14:15, 16:17, 18:21, 22:26, 27:30), ptype = integer(0), class = c("vctrs_list_of", 
    "vctrs_vctr", "list"))), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -8L), .drop = TRUE))

''']

I have attempted to transform the data and produce rarefaction curves, importantly, taking the two sites and creating individual datasets, then plotting independently. I am not sure if there is a way to combine them and do it all at once. OR whether there is a quicker and neater way to get at the end result.

# Make "wide"
dat_rare.w <-
  dat_rare %>% 
  pivot_wider(names_from = SPECIES, values_from = n, values_fill = 0)
dat_rare.w

dat_rare.w.foss = 
  dat_rare.w %>% 
  filter(SITE=="FOSSIL_CLIFFS") %>% 
  ungroup() %>% 
  select(!SITE)


dat_rare.w.foss = 
  dat_rare.w %>% 
  filter(SITE=="PAINTED_CLIFFS") %>% 
  ungroup() %>% 
  select(!SITE)


check = data.matrix(dat_rare.w.foss)
check

check2 <- check[,-1]
rownames(check2) <- check[,1]

#total number of species at each site (row of data)
S <- specnumber(check2)
# Number of INDIVIDULS per site (?)
raremax <- min(rowSums(check2)) # = 340;


# rarefy, w/ raremax as input (?)
Srare <- rarefy(check2, raremax)




# method = "rarefaction"
plot(specaccum(check2, method = "rarefaction"),main = "method = rarefaction")
legend("bottomleft", legend = "Note x-axis is # of plots")

# method = "rarefaction"
plot(specaccum(check2, 
               method = "rarefaction"), 
     xvar = "individuals",
     main = "rarefaction, xvar = individuals")
legend("bottomleft", legend = "Note x-axis has large range")

Answer 1

Have you tried using the rarecurve function in vegan? Does this do what you were wanting?

library(vegan)
library(data.table)
setDT(dat_rare)

# Aggregate at the site level
dat_aggregated <- dat_rare[, .(n = sum(n)), by = .(SITE, SPECIES)]

# make data wide
dat_wide <- dcast(dat_aggregated, SITE ~ SPECIES, value.var = "n", fill = 0)

# create a matrix
dat_matrix <- as.matrix(dat_wide [,-1, with = FALSE])

# Create rarefaction curves
rarecurve(dat_matrix, xlab = "Number of quadrats", ylab = "Species diversity", col = c("blue", "red"), lty = 1, lwd = 2, label = TRUE)