Working with Neotoma data in R

R workshop — neotoma2 · Bonn, 26.–27. 8. 2026 · COST PalaeOpen — PalaeoNoma 2026

Abstract

A complete walkthrough of the workshop scripts. You will install and load the toolchain, search the Neotoma Paleoecology Database for testate amoebae records across Europe, download a multi-proxy record from Prášilské jezero, convert counts to proportions, draw stratigraphic diagrams with cluster zonation, and relate a charcoal trend to diatom compositional change through ordination. Every code block from the scripts is reproduced verbatim and then explained.

1 How to use this handout

The workshop consists of five scripts that are meant to be run in order. Later scripts depend on objects created by earlier ones, so skipping one leaves you with object not found errors.

#	Script	What it does	Run it
0	___Init_project___.R	Installs every package, including two from GitHub	Once, before the workshop
1	00_setup.R	Attaches the packages	At the start of every session
2	01_example1.R	Multi-site search, filtering, mapping	Independent
3	01_example2.R	Single site, proportions, stratigraphic diagrams	Independent
4	01_example3.R	Smoothing, PCA, envfit, DCA	Requires objects from example 2

Example 3 continues example 2

01_example3.R uses site_charcoal_wide and site_diat_wide, which are built in 01_example2.R. Run example 2 first in the same session, or save its results with saveRDS() and load them back.

Rendering this document without running the code

Every code chunk here is live, and the Neotoma queries take time. To produce a code-only version of the handout, add this to the YAML header:

execute:
  eval: false

2 Part 0 — ___Init_project___.R: installing the toolchain

2.1 Listing the packages

package_list <- c(
  "dplyr",
  "ggplot2",
  "ggrepel",
  "ggvegan",
  "neotoma2",
  "remotes",
  "rioja",
  "vegan"
)

Keeping the names in a vector rather than writing eight install.packages() calls means the list can be reused, printed, and checked. What each one is for:

Package	Role in the workshop
dplyr	All data manipulation: filter(), group_by(), summarise(), mutate()
ggplot2	The ordination biplot in example 3
ggrepel	Non-overlapping species labels on that biplot
ggvegan	fortify() methods that turn vegan objects into data frames
neotoma2	The database interface itself
remotes	Installs the two packages that are not on CRAN
rioja	Stratigraphically constrained cluster analysis (chclust, bstick)
vegan	rda, decorana, envfit, vegdist, decostand

2.2 The installer function

install_packages <-
  function(pkgs_list) {
    # install all packages in the list from CRAN
    sapply(pkgs_list, utils::install.packages, character.only = TRUE)

    # install Github packages
    remotes::install_github("nsj3/riojaPlot")
    remotes::install_github("petrkunes/PolEco")
  }

install_packages(package_list)
rm(install_packages, package_list)

Two packages are not available from CRAN and must come from GitHub:

• riojaPlot (nsj3/riojaPlot) — Steve Juggins’ modern replacement for rioja::strat.plot(), used for both stratigraphic diagrams in example 2.
• PolEco (petrkunes/PolEco) — the workshop’s own helper package.

The final rm() removes the function and the vector from the global environment, keeping it clean. This is a good habit in teaching scripts: what remains in Environment is only what students actually need.

Two details worth knowing

sapply() installs one package at a time. A single vectorised call is faster and lets R resolve the whole dependency tree at once:

install.packages(package_list)     # one call, one dependency resolution

Only install what is missing if you want the script to be re-runnable without re-downloading everything:

missing <- setdiff(package_list, rownames(installed.packages()))
if (length(missing)) install.packages(missing)

GitHub installation needs build tools

remotes::install_github() compiles from source. Windows users need Rtools and macOS users need the Xcode command line tools (xcode-select --install). Without them the two GitHub installs fail while the CRAN ones succeed — and the failure appears only later, as there is no package called 'riojaPlot' in example 2.

If a GitHub install refuses to work, nsj3/riojaPlot also publishes binaries on R-universe, which need no compiler:

install.packages("riojaPlot",
                 repos = c("https://nsj3.r-universe.dev",
                           "https://cloud.r-project.org"))

3 Part 1 — 00_setup.R: loading the packages

package_list <- c(
  "dplyr",
  "ggplot2",
  "ggrepel",
  "ggvegan",
  "leaflet",
  "neotoma2",
  "rioja",
  "riojaPlot",
  "tidyr",
  "vegan"
)

invisible(sapply(package_list, library, character.only = TRUE))

rm(package_list)

library() normally takes an unquoted name (library(dplyr)). Inside sapply() the names arrive as character strings, which is exactly what character.only = TRUE permits — the argument that was decorative in ___Init_project___.R is essential here. sapply(package_list, library, ...) prints a long, noisy list of every attached package for each element. Wrap it with invisible() to stay quiet.

Note the differences from the installation list: remotes is gone (needed only once, at install time) and riojaPlot has appeared (installed from GitHub, so it could not be in the CRAN vector).

Load order decides which filter() you get

dplyr is loaded before neotoma2, and both export filter(). The package loaded later wins, so after this script filter() means neotoma2::filter().

That is why the example scripts write dplyr::filter() explicitly whenever they operate on a data frame. Get into the habit of qualifying both:

• dplyr::filter() — for data frames, including anything from samples()
• neotoma2::filter() — for sites, datasets and download objects

The symptom of getting it wrong:

Error in UseMethod("filter") :
  no applicable method for 'filter' applied to an object of class "sites"

vegan also masks a few names, and select() is contested between dplyr and several other packages. When in doubt, run conflicts(detail = TRUE) or simply qualify the call.

4 Part 2 — 01_example1.R: getting multiple datasets

4.1 Defining the search area

europe <- list(geoJSON = '{"type": "Polygon",
        "coordinates": [[
            [-30, 30],
            [70, 30],
            [70, 90],
            [-30, 90],
            [-30, 30]]]}')

Neotoma accepts a spatial constraint as GeoJSON. This is a rectangle covering Europe: longitude −30° to 70°, latitude 30° to 90°.

• Coordinates are [longitude, latitude] — longitude first, the opposite of the order most people say out loud.
• The ring must be closed: the last pair repeats the first. Omit it and the API rejects the query.
• list(geoJSON = ...) is a convention from the Neotoma tutorials, where a matching sf object is stored alongside as europe$sf. A bare string works just as well.

For anything more complicated than a box, use a real spatial object — sf objects are accepted directly:

eu_sf <- geojsonsf::geojson_sf(europe$geoJSON)   # string -> sf
my_area <- sf::st_read("data/study_area.shp")    # your own polygon
get_sites(loc = my_area, all_data = TRUE)

4.2 Finding sites

sites_eu <- get_sites(datasettype = "testate amoebae",
                      loc = europe$geoJSON,
                      all_data = TRUE)

This queries site-level metadata only — names, coordinates, elevation. No fossil data is transferred, which is why it is fast.

Argument	Meaning
datasettype	One kind of proxy. Also valid: "pollen", "diatom", "charcoal", "geochronologic", "vertebrate fauna", "plant macrofossil"
loc	Spatial filter, as GeoJSON text or an sf object
all_data	TRUE returns every match; the default returns only the first 25

Warning

all_data = TRUE is easy to forget and the default limit of 25 truncates your study region silently. Nothing in the output warns you that records are missing.

4.3 Inspecting what came back

neotoma2::summary(head(sites_eu, 20))

A sites object is a list of site objects with custom methods for printing, plotting and exporting. summary() gives one row per collection unit, showing which datasets each site holds and the age range covered. head(sites_eu, 20) keeps the table readable.

length(sites_eu)          # how many sites were found
as.data.frame(sites_eu)   # flat metadata table
sites_eu[[1]]             # the first site object
datasets(sites_eu)        # dataset metadata, once retrieved

4.4 From sites to datasets

datasets_eu <- sites_eu %>% get_datasets(all_data = TRUE)

The second call enriches the object with dataset-level metadata: dataset IDs, types, age ranges, analysts. Sites → datasets → downloads is the idiomatic order, because each step narrows the query before the expensive one.

You can already filter on metadata, before downloading anything:

datasets_eu %>% neotoma2::filter(!is.na(age_range_young))
datasets_eu %>% neotoma2::filter(lat > 50 & lat < 60)

Supported fields are siteid, sitename, lat, long, elev, altitude and datasettype.

4.5 A first map

neotoma2::plotLeaflet(datasets_eu)

plotLeaflet() returns an ordinary leaflet object, so anything the leaflet package can do can be piped onto it. Clicking a marker shows the site name, ID and a link to Neotoma Explorer. plot(datasets_eu) also works offline, but without any geographic context.

4.6 Downloading the data

sites_data    <- sites_eu %>% get_downloads()
sites_samples <- sites_data %>% samples()

get_downloads() is the expensive call — it retrieves every sample, count and chronology for the whole set of sites. Run it once, never inside a loop.

samples() flattens the nested result into a long-format data frame: one row per taxon per sample, with all context repeated on every row. That is why it is roughly 37 columns wide.

dim(sites_samples)
#> [1] 5024   39
names(sites_samples)
#>  [1] "age"             "agetype"         "ageolder"        "ageyounger"     
#>  [5] "chronologyid"    "chronologyname"  "units"           "value"          
#>  [9] "context"         "element"         "taxonid"         "symmetry"       
#> [13] "taxongroup"      "elementtype"     "variablename"    "ecologicalgroup"
#> [17] "analysisunitid"  "sampleanalyst"   "sampleid"        "depth"          
#> [21] "thickness"       "samplename"      "datasetid"       "database"       
#> [25] "datasettype"     "age_range_old"   "age_range_young" "age_units"      
#> [29] "recdatecreated"  "datasetnotes"    "siteid"          "sitename"       
#> [33] "lat"             "long"            "area"            "sitenotes"      
#> [37] "description"     "elev"            "collunitid"

The columns you will actually use:

Column	Contents
siteid, sitename	Site identity
datasetid	One dataset = one proxy from one collection unit
sampleid, depth, age	Which level in the core
variablename	Taxon name
value	Count or measurement
units	NISP for counts, but also number, ml, % for laboratory variables
ecologicalgroup	Functional group code — TRSH, UPHE, TEAM, ALGA, …
elementtype	What was counted — pollen, spore, test, stomate
taxongroup	Broad taxonomic group — Diatoms, Vascular plants, …
datasettype	Which proxy the row belongs to

Long format is not a species matrix

Absences are not zeros; they are missing rows. Any analysis needing a site-by-taxon matrix must reshape with an explicit fill — see Section 5.8.

4.7 Keeping the richer datasets

sites_samples <- sites_samples %>%
  group_by(datasetid) %>%
  filter(n_distinct(variablename) > 5) %>%
  ungroup()

Datasets with a handful of taxa are usually surface samples, pilot counts or training-set entries, and they distort comparisons. This block removes them.

1. group_by(datasetid) splits the table into one group per dataset.
2. filter() on a grouped data frame evaluates its condition per group and keeps or discards whole groups.
3. n_distinct(variablename) counts unique taxa within the current group.
4. ungroup() releases the grouping. Forgetting this is one of the most common sources of baffling dplyr results later in a script.

Look at the distribution before choosing a threshold

richness <- sites_samples %>%
  distinct(datasetid, variablename) %>%
  count(datasetid, name = "n_taxa")

summary(richness$n_taxa)
hist(richness$n_taxa, breaks = 30,
     main = "Taxa per dataset", xlab = "number of taxa")

Consider dplyr::filter(value > 0) first, so taxa recorded as present-but-zero do not count towards the total.

Overwriting the object

The result is assigned back to sites_samples, replacing the original. That saves memory but means the discarded records are gone — recovering them requires downloading again. In teaching code prefer a new name, e.g. samples_rich.

4.8 Subsetting the sites object

sites_filter <- sites_samples %>%
  distinct(siteid) %>%
  pull(siteid)

sites_filter <- neotoma2::filter(sites_data, siteid %in% sites_filter)

The filtering above happened on a data frame, so sites_data still holds every site. These lines carry the decision back to the sites object, which is what plotLeaflet() needs.

• distinct(siteid) reduces the long table to unique IDs.
• pull(siteid) converts that one-column data frame into a plain vector — without it, %in% would be comparing against a data frame.
• neotoma2::filter() applies the vector to the nested object; dplyr::filter() cannot handle a sites object.

The same name used for two different things

sites_filter is first a numeric vector and then a sites object. It runs, because R evaluates the right-hand side before assigning, but the block cannot be re-executed line by line and it is hard to read. Two names are clearer:

keep_ids     <- sites_samples %>% distinct(siteid) %>% pull(siteid)
sites_filter <- neotoma2::filter(sites_data, siteid %in% keep_ids)

neotoma2::plotLeaflet(sites_filter)

How many sites did the richness criterion remove?

length(sites_data)
#> [1] 38
length(sites_filter)
#> [1] 6

4.9 Counting taxa and scaling the markers

sites_rich <- sites_samples %>%
  group_by(siteid, sitename, lat, long) %>%
  summarise(n_taxa = n_distinct(variablename), .groups = "drop")

head(sites_rich)

Grouping by all four columns keeps sitename, lat and long in the output even though the aggregation is really over siteid — they are constant within a site, so this is safe and avoids a later join. .groups = "drop" returns an ungrouped table and suppresses dplyr’s message.

neotoma2::plotLeaflet(sites_filter) %>%
  leaflet::addCircleMarkers(
    data = sites_rich, lng = ~long, lat = ~lat,
    radius = ~sqrt(n_taxa) * 2, color = "firebrick",
    stroke = FALSE, fillOpacity = 0.7,
    label = ~paste0(sitename, ": ", n_taxa, " taxa")
  )

Reading the arguments:

• lng / lat — coordinates, taken from data = sites_rich.
• ~ — leaflet’s formula notation: “look this name up in data”.
• radius = ~sqrt(n_taxa) * 2 — the square root, so that circle area rather than radius is proportional to richness. Scaling radius linearly exaggerates large values badly.
• stroke = FALSE, fillOpacity = 0.7 — no outline, semi-transparent, so overlapping sites stay readable.
• label — text on hover. Use popup for click-activated text; popup accepts HTML, e.g. popup = ~paste0("<b>", sitename, "</b>").

Adding the search polygon

eu_sf <- geojsonsf::geojson_sf(europe$geoJSON)

neotoma2::plotLeaflet(sites_filter) %>%
  leaflet::addPolygons(data = eu_sf, color = "green", weight = 1, fill = FALSE)

Do not write addPolygons(map = ., data = ...) inside a pipe. Supplying the map both positionally and by name gives Error in dispatch(map, method, ...) : Invalid map parameter.

4.10 Exercises

1. Restrict the polygon to Scandinavia (longitude 4–32°, latitude 55–72°). How many testate amoebae sites remain?
2. Re-run with datasettype = "pollen" for one small country. What happens to the running time of get_downloads(), and why?
3. Raise the threshold from 5 to 15 taxa. Plot the distribution of taxa per dataset to justify a defensible cut-off.
4. Add elevation to the hover label: "Site name: 12 taxa, 340 m a.s.l.".
5. Colour the markers by richness with leaflet::colorNumeric() and add addLegend().
6. Which single taxon occurs at the most sites? Hint: count() on distinct(siteid, variablename).

5 Part 3 — 01_example2.R: stratigraphic plotting

5.1 Downloading one site

site_ds <- get_datasets(get_sites(sitename = "Prášilské jezero"), all_data = TRUE)

plotLeaflet(site_ds)

Searching by name instead of by area. sitename matching is case-insensitive and accepts % as a wildcard, so get_sites(sitename = "Prášil%") would also work — useful when you are unsure about diacritics.

The two calls are nested rather than piped; get_sites(...) %>% get_datasets(all_data = TRUE) is equivalent and reads left to right.

datasets(site_ds)
#>  datasetid                            database       datasettype age_range_old
#>      47517            European Pollen Database    geochronologic            NA
#>      47518 European Plant Macrofossil Database plant macrofossil       11985.0
#>      47519            European Pollen Database          charcoal       11954.0
#>      47520            European Pollen Database            pollen       11954.0
#>      55586            European Pollen Database        chironomid       11451.0
#>      55585            European Pollen Database            diatom       11426.0
#>      47521            European Pollen Database    geochronologic            NA
#>      47522 European Plant Macrofossil Database plant macrofossil       10641.1
#>      47523            European Pollen Database    geochronologic            NA
#>      47524 European Plant Macrofossil Database plant macrofossil        9679.1
#>  age_range_young age_units recdatecreated notes
#>               NA      <NA>     2020-03-25  <NA>
#>            -58.0      <NA>     2020-03-25  <NA>
#>            -58.0      <NA>     2020-03-25  <NA>
#>            -58.0      <NA>     2020-03-25  <NA>
#>            -60.0      <NA>     2022-11-10  <NA>
#>            -60.0      <NA>     2022-11-09  <NA>
#>               NA      <NA>     2020-03-28  <NA>
#>            -60.5      <NA>     2020-03-28  <NA>
#>               NA      <NA>     2020-03-30  <NA>
#>            -58.4      <NA>     2020-03-30  <NA>

datasets() lists what this site actually offers. Prášilské jezero is a multi-proxy record, which is why it was chosen: pollen, diatoms and charcoal from the same core, plus geochronological data.

5.2 Filtering to the proxies we want

site_ds <- site_ds %>% neotoma2::filter(datasettype %in% c("diatom", "charcoal", "pollen"))

neotoma2::filter() on the nested object, with the namespace written out because 00_setup.R left dplyr’s filter() masked. Filtering before downloading matters: it avoids transferring geochronological and other datasets you will not use.

5.3 Downloading and flattening

site_rec  <- site_ds %>% get_downloads()
site_samp <- site_rec %>% samples()

site_samp now holds all three proxies stacked in one long table. Everything that follows is about pulling the proxies apart again and putting each on its own appropriate scale.

site_samp %>%
  count(datasettype, taxongroup, units) %>%
  head(15)

5.4 Selecting charcoal

site_charcoal <- samples(site_rec) %>% dplyr::filter(datasettype == "charcoal") %>%
  dplyr::select(depth, age, variablename, elementtype, units, value) %>%
  dplyr::mutate(variable = paste(variablename, elementtype, sep = " "), .keep = "unused") %>%
  dplyr::group_by(depth, age) %>%
  dplyr::arrange(depth)

Step by step:

• dplyr::filter(datasettype == "charcoal") — keep charcoal rows only.
• dplyr::select(...) — drop the ~30 columns that are constant for a single site, leaving only what is needed.
• paste(variablename, elementtype, sep = " ") — charcoal is reported as several distinct measurements (area, count, concentration) that share a variablename; pasting elementtype on makes each one uniquely identifiable. This is what produces the column name `Charcoal area` used in example 3.
• .keep = "unused" — a mutate() argument that removes the input columns once they have been consumed, here dropping variablename and elementtype. Handy, but easy to miss when reading the code.
• dplyr::arrange(depth) — order from the top of the core downwards, which stratigraphic plotting functions expect.

Note

group_by(depth, age) here has no aggregation after it, so it changes nothing except leaving the result grouped. It is harmless, but ungroup() before passing the table on would be tidier.

5.5 Selecting diatoms and pollen

site_samp_diat <- samples(site_rec) %>%
  dplyr::filter(taxongroup == "Diatoms")

site_samp_pollen <- samples(site_rec) %>%
  dplyr::filter(datasettype == "pollen")

Note the two different criteria. Diatoms are selected by taxongroup and pollen by datasettype — both work, but they are not equivalent:

• datasettype == "pollen" returns everything in the pollen dataset, which includes spores, algae, charcoal counted on the pollen slide and laboratory variables such as Lycopodium spike volumes.
• taxongroup == "Diatoms" returns only diatom taxa, excluding any laboratory variables in the diatom dataset.

The extra rows in the pollen selection are exactly why the pollen sum needs careful definition below.

5.6 Diatom proportions

site_diat_perc <- site_samp_diat %>%
  dplyr::group_by(depth, age) %>%
  dplyr::mutate(diatcount = sum(value, na.rm = TRUE)) %>%
  dplyr::group_by(variablename) %>%
  dplyr::mutate(prop = value / diatcount) %>%
  dplyr::arrange(desc(age))

For diatoms the denominator is simple: the total number of valves counted in that sample. mutate() after group_by(depth, age) adds the group total to every row of the group without collapsing the table — that is the difference between mutate() and summarise().

The second group_by(variablename) does not affect the arithmetic, because value / diatcount is an element-wise division that needs no grouping at all. It does leave the result grouped by taxon, which can surprise you later; a plain ungroup() would be equivalent and clearer.

arrange(desc(age)) sorts oldest-first.

Tip

prop is a proportion (0–1), not a percentage. Both diagrams below multiply by 100 at plotting time. Pick one convention and stick to it; mixing them is a classic source of axes labelled 0–1 that should read 0–100.

5.7 Pollen proportions

pollen_sum <- site_samp_pollen %>%
  filter(ecologicalgroup %in% c("TRSH", "UPHE"),
                    elementtype == "pollen",
                    units == "NISP") %>%
  group_by(sampleid) %>%
  summarise(pollen_sum = sum(value), .groups = "drop")

summary(pollen_sum$pollen_sum)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   382.0   508.0   534.5   551.5   582.2   775.0

The pollen sum is the palynological convention: percentages are expressed relative to a defined terrestrial sum, not to everything on the slide. Here that sum is trees and shrubs (TRSH) plus upland herbs (UPHE).

Each of the three conditions removes a specific hazard:

Condition	Why it is needed
ecologicalgroup %in% c("TRSH", "UPHE")	Excludes aquatics, spores, algae and fungal remains, whose abundance reflects local conditions rather than regional vegetation
elementtype == "pollen"	Excludes stomata. Conifer stomata are recorded under the same variablename as the pollen (Pinus, Picea), so without this filter those taxa are counted twice
units == "NISP"	Excludes laboratory variables (CHAR, LABO) measured in number or ml, which are not counts of palynomorphs

site_pollen_perc <- site_samp_pollen %>%
  filter(units == "NISP") %>%                 # drops CHAR and LABO
  group_by(sampleid, depth, age, variablename,
           ecologicalgroup, elementtype) %>%
  summarise(count = sum(value), .groups = "drop") %>%
  left_join(pollen_sum, by = "sampleid") %>%
  mutate(prop = count / pollen_sum)

The percentages are then calculated for all taxa — including spores, algae and aquatics — against that terrestrial sum. That is deliberate and standard: non-pollen palynomorphs are reported relative to the pollen sum but are not part of it, so their percentages can exceed 100% in total.

The group_by() / summarise(count = sum(value)) step collapses any repeated rows for the same taxon–element combination within a sample, which can occur after taxonomic harmonisation. The left_join() then attaches one sum per sample.

# terrestrial pollen must total exactly 100% in every sample
site_pollen_perc %>%
  filter(ecologicalgroup %in% c("TRSH", "UPHE"), elementtype == "pollen") %>%
  group_by(sampleid) %>%
  summarise(total = sum(prop) * 100, .groups = "drop") %>%
  summarise(min = min(total), max = max(total))

Always run this check

If the minimum and maximum are not both 100, your sum definition and your percentage calculation disagree. Do this before plotting, not after.

A judgement call

Whether indeterminate pollen (UNID) belongs in the sum is a matter of convention. Most European Holocene studies exclude it from the denominator but still report it as a percentage, which is what this code does.

5.8 Reshaping to wide format

site_diat_perc <- site_diat_perc %>%
  dplyr::select(depth, age, variablename, prop) %>%
  dplyr::mutate(prop = as.numeric(prop))

site_diat_wide <- tidyr::pivot_wider(site_diat_perc,
                                     id_cols = c(depth, age),
                                     names_from = variablename,
                                     values_from = prop,
                                     values_fill = 0)

site_pollen_wide <- site_pollen_perc %>%
  filter(elementtype %in% c("pollen", "spore")) %>%
  tidyr::pivot_wider(id_cols = c(depth, age), names_from = variablename,
                                       values_from = prop,
                                       values_fill = 0)

site_charcoal_wide <- tidyr::pivot_wider(site_charcoal,
                                         id_cols = c(depth, age),
                                         names_from = variable,
                                         values_from = value,
                                         values_fill = 0)

Every downstream tool — riojaPlot, vegdist, rda, decorana — needs a samples-by-taxa matrix, so all three proxies are pivoted.

Argument	Role
id_cols = c(depth, age)	The columns that identify a row; everything else is dropped
names_from	Which column supplies the new column names
values_from	Which column fills the cells
values_fill = 0	Essential. Absences are missing rows in long format; without this they become NA and every distance calculation fails

The pollen table additionally keeps only pollen and spore elements, so stomata are excluded from the diagram while remaining available in site_pollen_perc if you want them.

values_fill = 0 versus a real NA

Filling with zero is right for counts of taxa that were looked for and not found. It would be wrong for a genuinely unmeasured variable. Since Neotoma stores only observed occurrences, zero is correct here.

5.9 Selecting taxa to plot

site_plot_taxa <- taxa(site_rec) %>% dplyr::filter(taxongroup == "Diatoms") %>%
  dplyr::arrange(desc(samples)) %>%
  head(n = 15)

site_plot_taxa

taxa() is a convenience function on a download object: it returns the unique taxa plus two extra columns, sites and samples, counting how many sites and how many samples each taxon appears in. Sorting on samples and taking the top 15 gives the most consistently present diatoms — a reasonable, reproducible way to keep a diagram readable.

Most frequent is not most abundant

samples counts occurrences, not counts. A taxon present at 1% in every sample ranks above one that dominates three samples and is absent elsewhere. If you want abundance instead:

site_diat_wide[, -1:-2] %>%
  colSums() %>%
  sort(decreasing = TRUE) %>%
  head(15) %>%
  names()

5.10 Constrained cluster analysis

site_clust <- chclust(vegdist(site_diat_wide[,-1:-2], "chord"), method = "coniss")
plot(site_clust, hang = -1)

bstick(site_clust)

Three things happen here.

site_diat_wide[,-1:-2] drops the first two columns — depth and age — so that only species data enters the distance calculation. The idiom is compact but depends entirely on column order; if you ever reorder or add an identifier column, it silently includes the wrong data. A safer equivalent:

site_diat_spec <- site_diat_wide %>% dplyr::select(-depth, -age)

vegdist(..., "chord") computes chord distance, appropriate for compositional data because it downweights the influence of the most abundant taxa and is not sensitive to differences in total count.

chclust(..., method = "coniss") is CONISS — stratigraphically constrained incremental sum of squares. Unlike ordinary clustering it may only merge adjacent samples, which is what makes the result interpretable as biostratigraphic zones. plot(..., hang = -1) aligns the leaves at the bottom.

bstick() compares the variance explained by each successive split with a broken-stick model. Where the observed curve drops below the broken-stick curve, additional zones are no longer distinguishable from random — that crossing point is how you choose the number of zones rather than eyeballing the dendrogram.

5.11 The diatom diagram

riojaPlot(site_diat_wide[,-1:-2]*100, site_diat_wide[,1:2],
          selVars = site_plot_taxa$variablename,
          scale.percent = TRUE,
          sec.yvar.name="age",
          plot.sec.axis = TRUE,
          srt.xlabel = 60,
          xRight = 0.85) |>
  addRPClust(site_clust) |>
  addRPClustZone(site_clust, col = "red")

riojaPlot() takes two tables: the species data first, then the depth and age scale second. The *100 converts proportions to percentages at the last possible moment.

Argument	Effect
selVars	Which taxa to draw, and in what order — the top 15 from taxa()
scale.percent = TRUE	Gives every curve the same percentage scale, so bar widths are comparable between taxa
sec.yvar.name = "age"	The second column of the scale table to use as a secondary axis
plot.sec.axis = TRUE	Actually draws that age axis on the right
srt.xlabel = 60	Rotates taxon names by 60°, which is how you fit long names in
xRight = 0.85	Reserves the right 15% of the figure — here, for the dendrogram

addRPClust() appends the CONISS dendrogram as an extra panel, and addRPClustZone() draws horizontal zone boundaries across all curves.

Two pipes in one expression

This block uses R’s native pipe |> while the rest of the script uses magrittr’s %>%. They are interchangeable for simple chaining. The difference that matters: |> has no . placeholder, so x |> f(data = .) is an error where x %>% f(data = .) works. riojaPlot’s helper functions are designed to take the plot object as their first argument, so |> is fine here.

5.12 The pollen diagram

site_ecol_groups <- site_pollen_perc %>%
  select(variablename, ecologicalgroup) %>%
  rename(Group = ecologicalgroup, Name = variablename) %>% mutate(Group = factor(Group, levels = c("TRSH", "UPHE"))) %>%
  distinct()

head(site_ecol_groups)

riojaPlot can group taxa into blocks — trees and shrubs, then herbs — but it expects a lookup table with columns named exactly Name and Group, hence the rename().

factor(Group, levels = c("TRSH", "UPHE")) does two jobs: it fixes the display order, so trees come before herbs as in a conventional diagram, and it converts every other group to NA, effectively excluding aquatics and spores from the grouping.

site_plot_taxa <- taxa(site_rec) %>% dplyr::filter(taxongroup == "Vascular plants") %>%
  dplyr::arrange(desc(samples)) %>%
  head(n = 15)

site_plot_taxa is reused

This overwrites the diatom selection made earlier. If you go back and re-run the diatom diagram now, it will try to plot vascular plant names against diatom data and fail. Use distinct names — taxa_diat and taxa_pollen — when teaching.

riojaPlot(site_pollen_wide[,-1:-2]*100, site_pollen_wide[,1:2],
          selVars = site_plot_taxa$variablename,
          scale.percent = TRUE,
          sec.yvar.name="age",
          plot.sec.axis = TRUE,
          plot.groups = TRUE,
          groups = site_ecol_groups,
          plot.cumul = TRUE,
          srt.xlabel = 60
          )

Two new arguments relative to the diatom diagram: plot.groups = TRUE with groups = site_ecol_groups turns on the grouped layout, and plot.cumul = TRUE adds a cumulative summary curve showing the total percentage of each group — the familiar filled silhouette on the left of a pollen diagram.

5.13 Exercises

1. Add the charcoal record as a further panel. site_charcoal_wide is already in wide format — what has to match for riojaPlot to accept it alongside the pollen?
2. Recalculate the pollen percentages without the elementtype == "pollen" filter in pollen_sum. How much do Pinus and Picea percentages change, and why?
3. Run bstick(site_clust) and decide how many zones are statistically supported. Redraw the diagram with only that many zone lines.
4. Replace the “most frequent” taxon selection with “most abundant” using the colSums() approach. Which taxa enter and leave the diagram?
5. Repeat the whole workflow for a different Neotoma site of your choice.

6 Part 4 — 01_example3.R: numerical analysis

Important

This script continues directly from example 2 and uses site_charcoal_wide and site_diat_wide. Run example 2 first.

6.1 Comparing smoothers on the charcoal record

plot(site_charcoal_wide$age, log(site_charcoal_wide$`Charcoal area`+1), ylim=c(0, 6), xlab = "Age cal BP", ylab = "CHAR")
lines(site_charcoal_wide$age, lowess(site_charcoal_wide$`Charcoal area`, f=0.1)$y, col="red")
lines(site_charcoal_wide$age, loess(site_charcoal_wide$`Charcoal area`~site_charcoal_wide$age, span = 0.2)$fitted, col="blue")
lines(site_charcoal_wide$age, smooth.spline(site_charcoal_wide$`Charcoal area`, df = 25)$y, col="green")
legend(1,6, legend = c("LOWESS", "LOESS", "smooth spline"), col = c("red", "blue", "green"), lty=1, cex=0.8)

The purpose is didactic: three ways of extracting a trend from a noisy accumulation record, drawn on top of each other for comparison.

• LOWESS (lowess) — locally weighted scatterplot smoothing, robust to outliers. f = 0.1 means each local fit uses 10% of the data.
• LOESS (loess) — the modern formula interface to the same idea, with span playing the role of f. It accepts a model formula, so covariates and weights are possible.
• Smoothing spline (smooth.spline) — a penalised spline whose flexibility is set by df; 25 degrees of freedom is fairly wiggly.

Backticks are needed around `Charcoal area` because the column name contains a space — a consequence of the paste() step in example 2.

Three inconsistencies in this block, worth discussing with students

The points and the lines are on different scales. The points are plotted as log(area + 1), but all three smoothers are fitted to and drawn on the raw values. With ylim = c(0, 6) the lines will mostly sit outside the visible range or hug the axis. Fix by smoothing the transformed variable:

char_age <- site_charcoal_wide$age
char_log <- log(site_charcoal_wide$`Charcoal area` + 1)

plot(char_age, char_log, ylim = c(0, 6),
     xlab = "Age cal BP", ylab = "log CHAR")
lines(char_age, lowess(char_age, char_log, f = 0.1)$y, col = "red")
lines(char_age, loess(char_log ~ char_age, span = 0.2)$fitted, col = "blue")
lines(char_age, smooth.spline(char_age, char_log, df = 25)$y, col = "green")
legend("topright", legend = c("LOWESS", "LOESS", "smoothing spline"),
       col = c("red", "blue", "green"), lty = 1, cex = 0.8, bty = "n")

lowess() was given no x variable. lowess(y, f = 0.1) smooths y against its own index, i.e. against sample number, whereas loess was given age. With unevenly spaced samples the two curves are therefore not comparable. Pass x explicitly, as above.

legend(1, 6, ...) positions by data coordinates. Since the x axis is age in years, x = 1 is the extreme left edge and the legend may be clipped. Use a keyword — legend("topright", ...) — which always works.

6.2 Interpolating charcoal onto the diatom depths

site_char_loess <- loess(site_charcoal_wide$`Charcoal area`~site_charcoal_wide$depth, span = 0.20)
site_char_pred <- stats::predict(site_char_loess, site_diat_wide$depth)

summary(site_char_pred)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.01370 0.09561 0.13259 0.35027 0.29165 2.04463

Charcoal and diatoms were counted on different levels of the same core, so they cannot be related directly. Fitting a LOESS to charcoal and predicting it at the diatom depths puts both proxies on a common set of samples.

Note that this model uses depth, not age, as predictor — even though the plot above used age. Depth is the right choice: it is measured directly, whereas age comes from a chronological model with its own uncertainty, and both proxies share the same depth scale exactly.

Warning

LOESS does not extrapolate. Diatom depths outside the charcoal depth range come back as NA, which will propagate into envfit(). Check with sum(is.na(site_char_pred)) and decide explicitly whether to drop those samples.

6.3 Hellinger PCA of the diatom data

site_pca <- rda(decostand(site_diat_wide[,-1:-2], method = "hellinger"))
site_pca
#> 
#> Call: rda(X = decostand(site_diat_wide[, -1:-2], method = "hellinger"))
#> 
#>               Inertia Rank
#> Total          0.3296     
#> Unconstrained  0.3296   35
#> 
#> Inertia is variance
#> 
#> Eigenvalues for unconstrained axes:
#>     PC1     PC2     PC3     PC4     PC5     PC6     PC7     PC8 
#> 0.09200 0.03992 0.03158 0.01973 0.01794 0.01523 0.01231 0.00895 
#> (Showing 8 of 35 unconstrained eigenvalues)

rda() with no constraining variables is a PCA. Applying it to Hellinger-transformed data is the standard “transformation-based” approach: the Hellinger transformation (square root of relative abundance) makes Euclidean distance behave like a chord distance, so the ordination respects the compositional nature of the data while retaining PCA’s linear framework.

site_PCA_species <- site_diat_wide[,-1:-2] %>%
  select(order(colSums(.), decreasing = T))

This reorders the columns by total abundance, most abundant first, purely so that the 15 most abundant taxa can be picked off with head(colnames(...), 15) below. The . is magrittr’s placeholder — this line requires %>% and would fail with the native |> pipe.

6.4 Extracting scores

site_scores <- fortify(site_pca)
site_scores_species <- subset(site_scores, score == "species" & label %in% head(colnames(site_PCA_species), 15))

site_PCA_env <- envfit(site_pca, data.frame(CHAR = site_char_pred))
site_PCA_env_scores <- fortify(site_PCA_env)

fortify() from ggvegan converts a vegan object into a long data frame that ggplot2 can consume, with a score column distinguishing species from sites rows and one column per ordination axis.

envfit() fits the interpolated charcoal values onto the ordination as a vector, and reports how much of its variation the ordination axes explain. Print it to get the goodness of fit and a permutation p-value:

site_PCA_env
#> 
#> ***VECTORS
#> 
#>           PC1      PC2     r2 Pr(>r)   
#> CHAR -0.96836 -0.24957 0.3748  0.002 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> Permutation: free
#> Number of permutations: 999

Check the column names once

ggvegan reached CRAN as version 0.2.1 in February 2026 (CRAN), and that release changed fortify() to return tibbles with tidied column names. The script mixes lower case (pc1, pc2 for species) and upper case (PC1, PC2 for the envfit result), which suggests the two methods do not agree. Before running the plot, look:

names(site_scores)
names(site_PCA_env_scores)
head(site_scores)

and adjust the aes() calls to match. A version-proof alternative that avoids ggvegan entirely:

sp <- as.data.frame(vegan::scores(site_pca, display = "species", choices = 1:2))
sp$label <- rownames(sp)

ar <- as.data.frame(vegan::scores(site_PCA_env, display = "vectors"))
ar$label <- rownames(ar)

6.5 The biplot

ggplot() +
  geom_point(data = site_scores_species, aes(x=pc1, y=pc2), color = "#117733", size = 2) +
  geom_text_repel(data = site_scores_species, aes(x = pc1, y = pc2, label = label), color = "#117733", size = 5) +
  geom_segment(data = site_PCA_env_scores, aes(x = 0, y = 0, xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.2, "cm")), color = "#D55E00") +
  geom_text_repel(data = site_PCA_env_scores, aes(x = PC1, y = PC2, label = label), color = "#D55E00", size = 5) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkgray") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkgray") +
  theme_classic() +
  theme(
    panel.border = element_rect(color = "black", fill = NA, linewidth = 1),
    axis.line = element_blank(),
    axis.text = element_text(size = 10)
  ) +
  labs(x = "PC1", y = "PC2")

The plot is built from two data sets rather than one, which is why ggplot() is called empty and each layer supplies its own data.

• geom_point + geom_text_repel — the 15 most abundant diatom taxa. ggrepel pushes labels apart so they stay legible; plain geom_text would overplot badly.
• geom_segment with arrow() — the fitted charcoal vector, drawn from the origin. Its direction indicates the compositional gradient most strongly associated with charcoal; its length reflects the strength of that association.
• geom_vline / geom_hline at zero — the ordination origin, essential for reading a biplot.
• theme_classic() plus an explicit panel.border and axis.line = element_blank() gives a full box rather than two axis lines.

The colours #117733 (green) and #D55E00 (vermillion) are from colourblind-safe palettes — a good habit for published figures.

Two caveats when interpreting this figure

Arrow lengths are not on the species scale. envfit vectors are scaled for display only; comparing arrow length with species distances is meaningless. The usual practice is to multiply the arrow coordinates by a constant chosen to fill the plot, and to say so in the caption.

The samples are missing. Only species and the environmental vector are plotted. Adding the sample scores, coloured by CONISS zone from example 2, would connect this figure back to the stratigraphic diagram — see the exercises.

6.6 Checking gradient length with DCA

site_dca <- decorana(decostand(site_diat_wide[,-1:-2], method = "hellinger"))
site_dca
#> 
#> Call:
#> decorana(veg = decostand(site_diat_wide[, -1:-2], method = "hellinger")) 
#> 
#> Detrended correspondence analysis with 26 segments.
#> Rescaling of axes with 4 iterations.
#> Total inertia (scaled Chi-square): 2.0524 
#> 
#>                        DCA1   DCA2    DCA3    DCA4
#> Eigenvalues          0.2985 0.1344 0.09716 0.07879
#> Additive Eigenvalues 0.2985 0.1343 0.09604 0.07763
#> Decorana values      0.3067 0.1425 0.11246 0.06256
#> Axis lengths         1.6684 1.9264 1.78166 1.31981
plot(site_dca, display = "sites")

Detrended correspondence analysis reports the length of the first axis in standard deviation units of species turnover, which is the conventional diagnostic for choosing between linear and unimodal methods:

Axis 1 length	Interpretation
< 3 SD	Short gradient; linear methods (PCA, RDA) are appropriate
3–4 SD	Ambiguous; either family works
> 4 SD	Long gradient with much turnover; unimodal methods (CA, CCA) preferred

Running this after the PCA is a bit back-to-front pedagogically — the check justifies the choice of method, so in your own work do it first. Note also that decorana is normally applied to untransformed or square-root-transformed abundances; combining it with a Hellinger transformation is unusual and will shorten the reported gradient, so interpret the number with that in mind.

6.7 Exercises

1. Fix the smoother plot so that points and lines are on the same scale, and pass x to lowess(). Do the three smoothers now agree?
2. Vary span in loess() between 0.05 and 0.6. What happens to the trend, and how would you choose a value defensibly?
3. Add the sample scores to the biplot (score == "sites"), coloured by the CONISS zone each sample belongs to.
4. Is the charcoal vector significant? Report the envfit \(r^2\)$r^2$ and p-value, and repeat with permutations = 9999.
5. Run decorana() on untransformed diatom percentages and compare the axis lengths with the Hellinger version.
6. Replace the unconstrained PCA with a constrained rda(species ~ CHAR). How much variance does charcoal explain on its own?

7 Troubleshooting

Symptom	Cause and fix
no applicable method for 'filter' applied to an object of class "sites"	dplyr’s filter() is masking neotoma2’s. Write neotoma2::filter().
there is no package called 'riojaPlot'	The GitHub install in Init_project.R failed, usually for lack of build tools. Install from R-universe instead.
Error in dispatch(map, method, ...) : Invalid map parameter	The first argument was not a leaflet map — an empty sites object, or map = . used together with a pipe. Check length() and class().
Only 25 sites returned	all_data = TRUE is missing.
object 'site_diat_wide' not found	Example 3 was run without example 2 in the same session.
NA/NaN/Inf in foreign function call in vegdist	values_fill = 0 was omitted in pivot_wider(), leaving NAs.
Percentages do not total 100	The pollen sum and the percentage calculation use different filters. Re-run the check in Section 5.7.
object 'pc1' not found in the biplot	fortify() column names differ in your ggvegan version. See Section 6.4.
summarise() prints a regrouping message	Harmless; add .groups = "drop".
Results change after a group_by()	A grouping was never released. Add ungroup().
Column name with a space rejected	Wrap it in backticks: `Charcoal area`.

8 Session information

Include this in any handout or supplement so results can be reproduced against the same package versions.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Tahoe 26.5.2
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Europe/Warsaw
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] vegan_2.7-5      permute_0.9-10   tidyr_1.3.2      riojaPlot_0.1-25
#>  [5] rioja_1.0-7      neotoma2_1.0.12  leaflet_2.2.3    ggvegan_0.2.1   
#>  [9] ggrepel_0.9.8    ggplot2_4.0.3    dplyr_1.2.1     
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6       xfun_0.60          htmlwidgets_1.6.4  insight_1.5.2     
#>  [5] lattice_0.22-9     vctrs_0.7.3        tools_4.6.1        crosstalk_1.2.2   
#>  [9] generics_0.1.4     curl_7.1.0         parallel_4.6.1     tibble_3.3.1      
#> [13] proxy_0.4-29       cluster_2.1.8.2    pkgconfig_2.0.3    Matrix_1.7-5      
#> [17] KernSmooth_2.23-26 RColorBrewer_1.1-3 S7_0.2.2           assertthat_0.2.1  
#> [21] uuid_1.2-2         lifecycle_1.0.5    compiler_4.6.1     farver_2.1.2      
#> [25] stringr_1.6.0      sjmisc_2.8.11      progress_1.2.3     codetools_0.2-20  
#> [29] htmltools_0.5.9    class_7.3-23       yaml_2.3.12        jquerylib_0.1.4   
#> [33] pillar_1.11.1      crayon_1.5.3       MASS_7.3-65        classInt_0.4-11   
#> [37] nlme_3.1-169       sjlabelled_1.2.0   tidyselect_1.2.1   digest_0.6.39     
#> [41] stringi_1.8.9      sf_1.1-2           purrr_1.2.2        labeling_0.4.3    
#> [45] forcats_1.0.1      splines_4.6.1      fastmap_1.2.0      grid_4.6.1        
#> [49] cli_3.6.6          magrittr_2.0.5     e1071_1.7-17       withr_3.0.3       
#> [53] prettyunits_1.2.0  scales_1.4.0       rmarkdown_2.31     httr_1.4.8        
#> [57] hms_1.1.4          evaluate_1.0.5     knitr_1.51         mgcv_1.9-4        
#> [61] rlang_1.3.0        Rcpp_1.1.2         glue_1.8.1         DBI_1.3.0         
#> [65] geojsonsf_2.0.5    rstudioapi_0.19.0  jsonlite_2.0.0     R6_2.6.1          
#> [69] units_1.0-1