Doing a community ecology analysis using `npe_community_matrix` and {vegan}

Setup

Load the tidyverse library along with neonPlantEcology. There is a pre-installed dataset which contains the neonUtilities-generated list object for the Jornada Experimental Range (“JORN”) and the Santa Rita Experimental Range (“SRER”) from domain 14. This can be downloaded directly by uncommenting the following code, or just call data("D14").

# D14 <- npe_site_ids(domain = "D14") |>
#   npe_download_plant_div(sites = .)

data("D14")

Wrangle the data

Use npe_community_matrix to get the data into the needed format. By default, it aggregates annually to the (20m x 20m) plot scale. Use npe_cm_metadata to get all of the information from the rownames into a more interpretable format. Plot centroids, along with additional plot-level metadata, can be loaded with data("plot_centroids").

comm <- npe_community_matrix(D14)
#> Loading required namespace: lubridate

comm[1:4, 1:4]
#>                       2PLANT      ALLIO   AMAC       AMFI
#> JORN_001_plot_2017 1.2500000 0.00000000 0.0000 0.00000000
#> JORN_001_plot_2021 0.5833333 0.00000000 0.0000 2.00000000
#> JORN_001_plot_2022 0.0000000 0.08333333 0.0000 0.08333333
#> JORN_001_plot_2016 0.0000000 0.00000000 2.3125 0.00000000


data("plot_centroids")
plot_centroids <- sf::st_set_geometry(plot_centroids, NULL)
metadata <- npe_cm_metadata(comm) |>
  dplyr::left_join(plot_centroids)
#> Joining with `by = join_by(plotID)`
#> Warning in dplyr::left_join(npe_cm_metadata(comm), plot_centroids): Detected an unexpected many-to-many relationship between `x` and `y`.
#> ℹ Row 1 of `x` matches multiple rows in `y`.
#> ℹ Row 1001 of `y` matches multiple rows in `x`.
#> ℹ If a many-to-many relationship is expected, set `relationship =
#>   "many-to-many"` to silence this warning.

`vegan` analysis

first we’ll do a species accumulation curve for each site.


# checking the metadata file to see where JORN stops and SRER begins
which(metadata$site == "JORN") |> range()
#> [1]   1 354

which(metadata$site == "SRER") |> range()
#> [1] 355 729


# performing and plotting the species accumulation curve analysis
sp_jorn <- vegan::specaccum(comm[1:152,])
sp_srer <- vegan::specaccum(comm[152:347,])
plot(sp_srer, col = "blue")
plot(sp_jorn, col = "gold", add=T)

NMDS

The following code reproduces figure 3 in Mahood et al 2024 (reference ).

nmds <- metaMDS(comm,trace = F)

nmds_sites <- nmds$points |>
  as_tibble(rownames = "rowname") |>
  left_join(metadata)
#> Joining with `by = join_by(rowname)`


ggplot(nmds_sites, aes(x=MDS1, y=MDS2, color = eventID, size = elevation, shape = site)) +
  geom_point() +
  theme_classic() +
  scale_color_brewer(palette = "Set1") +
  theme(panel.background = element_rect(fill=NA, color= "black"))

getting summary info by site using `npe_summary`

Species richness by biogeographical origin

di <- npe_summary(D14, scale = "site", timescale = "all")

dj<- di |>
  dplyr::select(site, eventID, starts_with("nspp")) |>
  tidyr::pivot_longer(cols = starts_with("nspp")) |>
  dplyr::filter(!name %in% c("nspp_notexotic", "nspp_total")) |>
  dplyr::transmute(origin = str_remove_all(name, "nspp_"),
            nspp = value,
            site=site)

p1<-ggplot(dj, aes(x=nspp, y=site, fill = origin)) +
  geom_bar(stat = "identity", color = "black", lwd = .2) +
  xlab("Species Richness") +
  ylab("Site")

Relative cover by biogeographical origin

 dk<- di |>
  dplyr::select(site, eventID, starts_with("rel_cover")) |>
  pivot_longer(cols = starts_with("rel_cover")) |>
  filter(!name %in% c("rel_cover_notexotic", "rel_cover_total")) |>
  transmute(origin = str_remove_all(name, "rel_cover_"),
            relative_cover = value,
            site=site)

p2<- ggplot(dk, aes(x=relative_cover, y=site, fill = origin)) +
  geom_bar(stat = "identity", color = "black", lwd = .2) +
  xlab("Relative Cover")

Alpha biogeographical diversity

dl <- di|>
  dplyr::select(site, eventID, starts_with("shannon")) |>
  pivot_longer(cols = starts_with("shannon")) |>
  filter(!name %in% c("shannon_notexotic", "shannon_total", "shannon_family")) |>
  transmute(origin = str_remove_all(name, "shannon_"),
            shannon = value,
            site=site)

p3<- ggplot(dl, aes(x = shannon, y = site, fill = origin)) +
  geom_bar(stat = 'identity', color = "black", lwd = .2) +
  xlab("Shannon-Weiner Diveristy")

Make a multipanel figure

ggpubr::ggarrange(p1, 
                  p2 + theme(axis.title.y = element_blank(), 
                                 axis.text.y = element_blank(), 
                                 axis.ticks.y = element_blank()), 
                  p3 + theme(axis.title.y = element_blank(), 
                                 axis.text.y = element_blank(), 
                                 axis.ticks.y = element_blank()), 
                  nrow=1, ncol=3, 
                  common.legend = TRUE, 
                  widths = c(1.35,1,1))

Family cover through time using `npe_longform`

data("D14")
lf <- npe_longform(D14, scale = "plot", timescale = "annual")

lf_f <- lf |>
  mutate(family = ifelse(!family %in% c("Poaceae", "Fabaceae"), "Other", family)) |>
  group_by(site, plotID, eventID, family) |>
  summarise(cover = sum(cover, na.rm=T)) |>
  ungroup()
#> `summarise()` has grouped output by 'site', 'plotID', 'eventID'. You can
#> override using the `.groups` argument.

  
ggplot(lf_f, aes(x=eventID, y=cover, fill = family)) +
  geom_boxplot(position="dodge") +
  facet_wrap(~site, scales = "free_x") +
  theme_bw() +
  scale_fill_brewer(palette = "Set1", name = "Family") +
  ylab("Cover") +
  xlab("eventID (Annual Timestep)")

Using neonPlantEcology

Doing a community ecology analysis using npe_community_matrix and {vegan}