Plotting Rarefaction Curves in R (UPDATE)

In this video I demonstrate how to do rarefaction using the rarefy () function in the vegan package in R and plot rarefaction curves.
Here's the code I ran in this video, but NOTE that TH-cam does not allow angled brackets in video descriptions, so I could not type the assignment operator ("less than sign" plus "hyphen"):
Example code for generating and plotting a rarefaction curve
John Sakulich
5 August 2024
Read in a data set
These data are from Colorado Breading Bird Survey
for all individuals counted in the BBS in 1978 and 2018
birds -read.csv("Birds!.csv")
Calculate Species richness for each survey year (1978 and 2018)
abundance.2018 - sum(birds$abundance_2018)
abundance.1978 - sum(birds$abundance_1978)
richness.2018 - sum(birds$abundance_2018 != 0)
richness.1978 - sum(birds$abundance_1978 != 0)
Load in the vegan package
library (vegan)
rarefied.2018 - rarefy(birds$abundance_2018, 15905, se = TRUE)
Generate values for 2018 BBS data (221 total species)
rare.2018 - rarefy(birds$abundance_2018,seq(1, 63594, by = 250) , se = TRUE)
trare.2018 - t(rare.2018)
rare.2018.df - as.data.frame(trare.2018)
rare.2018.df$se.high - rare.2018.df$.S + rare.2018.df$.se
rare.2018.df$se.low - rare.2018.df$.S - rare.2018.df$.se
rare.2018.df$N - seq(1, 63594, by = 250)
Generate values for 1978 BBS data (143 total species)
rare.1978 - rarefy(birds$abundance_1978,seq (1, 15905, by = 64) , se = TRUE)
trare.1978 - t(rare.1978)
rare.1978.df - as.data.frame(trare.1978)
rare.1978.df$se.high - rare.1978.df$.S + rare.1978.df$.se
rare.1978.df$se.low - rare.1978.df$.S - rare.1978.df$.se
rare.1978.df$N - seq (1, 15905, by = 64)
Plot the 2018 BBS rarefaction curve with standard errors
plot (rare.2018.df$N, rare.2018.df$.S, type = "l", xlim = c(1, 10000), lwd = 2, col = "firebrick4", cex = 2, xlab = "abundance level", ylab = "species richness", main = "Colorado Breeding Bird Survey 1978 and 2018")
lines (rare.2018.df$N, rare.2018.df$se.high, type = "l", lty = "dashed", col = "firebrick4")
lines (rare.2018.df$N, rare.2018.df$se.low, type = "l", lty = "dashed", col = "firebrick4")
Add the 1978 BBS data to the plot
lines (rare.1978.df$N, rare.1978.df$.S, type = "l", lwd = 2, lty = "solid", cex = 2, col = "dodgerblue")
lines (rare.1978.df$N, rare.1978.df$se.high, type = "l", lty = "dashed", col = "dodgerblue")
lines (rare.1978.df$N, rare.1978.df$se.low, type = "l", lty = "dashed", col = "dodgerblue")
legend (40000, 75, legend = c("BBS 2018", "BBS 1978"), col = c("firebrick4","dodgerblue"), lty=1:1, lwd = 2)