Cleaning species names before you count

One species written six ways becomes six species in your richness count. Clean names with trimws, tolower and a lookup table, and see why fuzzy matching is a trap.

length(unique(species)) is the first line of almost every diversity analysis, and it trusts the spelling completely. Write Festuca rubra one day and festuca rubra the next, add a stray space, let someone abbreviate it to F. rubra, and R counts four species where the field held one. Nobody decided this; it is the sum of small keystrokes across a season, and it inflates richness in exactly the plots where you worked hardest.

This post takes one messy notebook, cleans it in the order that actually works, and stops at the point where string tidying can do no more. The last step is a lookup table, and the tempting shortcut past it, fuzzy matching, is where the real damage happens.

One notebook, many spellings

The notebook below is synthetic: twelve plots, six real species, and a realistic rate of transcription slips. The generator is in the source; what matters is the column it produced, which is what a shared field sheet looks like after a season of hands.

head(nb[, c("plot", "recorded")], 8)

  plot             recorded
1  P01   Trifolium pratense
2  P01        eFstuca rubra
3  P01  Plantago lanceolata
4  P01          Carex flava
5  P01  Plantago lanceolata
6  P02  Plantago lanceolata
7  P02 Achillea millefolium
8  P02  Plantago lanceolata

nrow(nb)

[1] 298

Now the line that starts every analysis, run on the column as typed:

length(unique(nb$recorded))

[1] 49

49 species, from a field that held six. Here they are, and every one is a variant of the same six binomials:

sort(unique(nb$recorded))

 [1] "A. millefolium"          "Achillea  millefolium"  
 [3] "achillea millefolium"    "Achillea millefolium"   
 [5] "Achillea millefolium "   "Achillea millefolium L."
 [7] "C. flacca"               "C. flava"               
 [9] "carex flacca"            "Carex flacca"           
[11] "Carex flacca "           "Carex flacca Schreb."   
[13] "carex flava"             "Carex flava"            
[15] "Carex flava "            "Carex flava L."         
[17] "Carxe flacca"            "Carxe flava"            
[19] "eFstuca rubra"           "F. rubra"               
[21] "Festcua rubra"           "Festuca  rubra"         
[23] "Festuca rubar"           "festuca rubra"          
[25] "Festuca rubra"           "Festuca rubra "         
[27] "Festuca rubra L."        "Fetsuca rubra"          
[29] "Fsetuca rubra"           "lPantago lanceolata"    
[31] "P. lanceolata"           "Plantago  lanceolata"   
[33] "Plantago alnceolata"     "Plantago lanceloata"    
[35] "Plantago lanceolaat"     "plantago lanceolata"    
[37] "Plantago lanceolata"     "Plantago lanceolata "   
[39] "Plantago lanceolata L."  "T. pratense"            
[41] "Trifoilum pratense"      "Trifolium  pratense"    
[43] "Trifolium partense"      "Trifolium praetnse"     
[45] "trifolium pratense"      "Trifolium pratense"     
[47] "Trifolium pratense "     "Trifolium pratense L."  
[49] "Trifoluim pratense"     

Trailing spaces, double spaces, lower case, transposed letters, abbreviations, appended authorities. Each is a different string to R, and unique() has no way to know they are the same plant.

Clean in the order that works

Three string operations remove most of it, and the order matters because each assumes the last has run.

x <- nb$recorded
step <- c(raw = length(unique(x)))

x <- trimws(x)
step["trim whitespace"] <- length(unique(x))

x <- gsub("\\s+", " ", x)
step["squish doubles"] <- length(unique(x))

x <- tolower(x)
step["lower case"] <- length(unique(x))

step

            raw trim whitespace  squish doubles      lower case 
             49              43              39              33 

trimws strips the leading and trailing spaces, gsub("\\s+", " ", x) collapses any run of internal whitespace to one space, and tolower removes the capitalisation differences. Do tolower first and the trailing spaces survive; collapse before you trim and the edge spaces are still there. This order, outside in then case last, is the one that leaves nothing for the next step to trip on.

library(ggplot2)

steps <- data.frame(
  stage = factor(c(names(step), "lookup table"), levels = c(names(step), "lookup table")),
  n = c(step, 6),
  kind = c("string", "string", "string", "string", "decision")
)

ggplot(steps, aes(stage, n, group = 1)) +
  geom_hline(yintercept = 6, linetype = "dashed", colour = "#b5534e") +
  geom_step(colour = "#93a87f", linewidth = 1) +
  geom_point(aes(colour = kind), size = 3.6) +
  geom_text(aes(label = n), vjust = -1, size = 3.4, colour = "#16241d") +
  annotate("text", x = 1.4, y = 6, label = "true richness", vjust = 1.6,
           colour = "#b5534e", size = 3.2, hjust = 0) +
  scale_colour_manual(values = c(string = "#275139", decision = "#b5534e")) +
  scale_y_log10(limits = c(5, 60)) +
  labs(x = NULL, y = "distinct names (log scale)", colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top", panel.grid.minor = element_blank(),
        axis.text.x = element_text(angle = 20, hjust = 1))

Figure 1: The distinct name count after each cleaning operation. Whitespace and case tidying remove a lot, but the count settles well above the six species that were actually present. The last drop, from the plateau down to six, is not a string operation at all; it is a decision about which spellings mean which species.

The chain gets from 49 to 33. That is most of the way, and it is where a lot of analyses stop. It is also still 27 species too many, and no amount of further string tidying closes the gap, because the remaining variants are not typos of layout. They are abbreviations and misspellings, and f. rubra is not mechanically closer to festuca rubra than it is to anything else.

The lookup table is the honest tool

The remaining 33 strings have to be mapped to accepted names by a decision, one row per variant, written down where it can be checked.

clean <- tolower(gsub("\\s+", " ", trimws(nb$recorded)))

lookup <- data.frame(recorded = sort(unique(clean)))
head(lookup, 6)

                 recorded
1          a. millefolium
2    achillea millefolium
3 achillea millefolium l.
4               c. flacca
5                c. flava
6            carex flacca

For a real dataset you write the accepted name beside each variant by hand, or against a taxonomic authority. Since this notebook is synthetic, we already know the intended species for every record, so we can build the map from that and see the result:

key <- tapply(tolower(nb$true_name), clean, function(v) v[1])
accepted <- unname(key[clean])
length(unique(accepted))

[1] 6

Six. The lookup is the step that reached the truth, and the reason it works where the string chain stalled is that it encodes knowledge the strings do not contain: that carxe flacca and c. flacca and carex flacca are the same sedge. A join does the same thing at scale; the key here is the small version of it.

Keep the lookup table as a file in the project. It is the record of every naming decision made, it is auditable, and next season it starts already populated.

Why this is worse where you worked harder

The inflation is not spread evenly. A plot with five records has little room to misspell; a plot with sixty has sixty chances.

per <- do.call(rbind, lapply(plots, function(p) {
  s <- nb[nb$plot == p, ]
  data.frame(plot = p, records = nrow(s),
             apparent = length(unique(s$recorded)),
             true = length(unique(s$true_name)))
}))
per$inflation <- per$apparent - per$true
per

   plot records apparent true inflation
1   P01       5        4    4         0
2   P02       7        5    4         1
3   P03       9        7    5         2
4   P04      12       10    6         4
5   P05      15        8    6         2
6   P06      18        8    4         4
7   P07      22        6    5         1
8   P08      27       16    6        10
9   P09      33       14    6         8
10  P10      40       16    6        10
11  P11      50       20    6        14
12  P12      60       18    6        12

fit <- lm(inflation ~ records, data = per)
c(slope = unname(round(coef(fit)[2], 3)),
  p_value = signif(summary(fit)$coefficients[2, 4], 2),
  cor_apparent = round(cor(per$records, per$apparent), 2),
  cor_true = round(cor(per$records, per$true), 2))

       slope      p_value cor_apparent     cor_true 
     2.5e-01      4.1e-05      9.0e-01      6.2e-01 

Each extra record adds about 0.25 of a spurious species, and the effect is strong: the p-value in the output above is well under 0.001 on twelve plots. The damage is in the last two numbers. Apparent richness correlates with sampling effort at 0.9; true richness correlates at 0.62. The typos have manufactured an effort gradient in richness that is not in the ecology, and effort is rarely random across plots. If the large plots are the productive ones, the misspellings hand you a productivity to diversity relationship for free, and it is entirely an artefact of the keyboard.

long <- rbind(
  data.frame(records = per$records, richness = per$apparent, count = "as typed"),
  data.frame(records = per$records, richness = per$true, count = "cleaned")
)

ggplot(long, aes(records, richness, colour = count)) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.8) +
  geom_point(size = 3) +
  scale_colour_manual(values = c("as typed" = "#b5534e", "cleaned" = "#275139")) +
  labs(x = "records in plot (sampling effort)", y = "species richness", colour = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "top", panel.grid.minor = element_blank())

`geom_smooth()` using formula = 'y ~ x'

Figure 2: Apparent richness, from the names as typed, against sampling effort, with true richness for comparison. The apparent count climbs with effort while the true count is nearly flat. A richness gradient that is really a typo gradient will survive into any model that does not clean the names first.

The trap: fuzzy matching

Writing a lookup table by hand feels like work, and the obvious escape is to let the computer decide which strings are close enough. R measures string distance with adist, the number of single character edits between two strings.

adist("festuca rubra", "fetsuca rubra")

     [,1]
[1,]    2

Two edits for a transposed pair of letters. So match anything within a distance of two to the nearest accepted name and the typos collapse on their own. The problem arrives when two of the accepted names are themselves close:

adist("carex flacca", "carex flava")

     [,1]
[1,]    2

Also two. Carex flacca and Carex flava are two different sedges, and they sit at the same edit distance as a typo and its correct spelling. Any threshold loose enough to fix fetsuca is loose enough to merge them, and there is no threshold that separates the real pair from the mistakes.

Here is what that costs when it happens quietly. Suppose last year’s reference list is used to snap this year’s names, and last year Carex flava was never recorded, so it is not on the list.

reference <- setdiff(tolower(pool), "carex flava")

snap <- function(s, ref, max_edit = 2) {
  d <- adist(s, ref)
  i <- which.min(d)
  if (d[i] <= max_edit) ref[i] else NA_character_
}
matched <- vapply(clean, snap, character(1), ref = reference)

flava_records <- sum(nb$true_name == "Carex flava")
flava_to_flacca <- sum(nb$true_name == "Carex flava" & matched %in% "carex flacca")
c(flava_records = flava_records, silently_relabelled = flava_to_flacca,
  richness_after = length(unique(na.omit(matched))))

      flava_records silently_relabelled      richness_after 
                 24                  18                   5 

Most Carex flava records fall within two edits of carex flacca, so 18 of the 24 get relabelled to it; the rest fail to match anything, and either way the species vanishes from the dataset. Richness drops to 5. Worse, the abundance of the wrong sedge is now inflated:

before <- sum(nb$true_name == "Carex flacca")
after <- sum(matched == "carex flacca", na.rm = TRUE)
c(flacca_before = before, flacca_after = after,
  per_cent_change = round(100 * (after / before - 1)))

  flacca_before    flacca_after per_cent_change 
             22              38              73 

The count for Carex flacca rises 73 per cent, built entirely from another species’ records, and nothing in the output flags it. A hand lookup would have hit the missing name and stopped: there is no row for carex flava, so it demands a decision instead of guessing. Fuzzy matching guesses, and the guess is silent.

Fuzzy matching has a place as a way to draft the lookup table, listing near matches for a human to confirm. It has no place applying itself unseen, because the one case it cannot handle, two real names a few letters apart, is common in exactly the genera where careful identification matters most.

What to carry to the next dataset

Richness and every diversity index built on it inherit the spelling of the names, so the names come first, before unique, before diversity, before any count. The workable order is trimws, then collapse internal whitespace, then tolower, and then a lookup table for everything that is left, kept as a file. Read the plateau on the funnel figure as the honest signal: when string cleaning stops removing names well above your known species total, the rest is identity, not layout, and it is a decision to be written down rather than a function to be called.

Where to go next

The cleaned names feed straight into Diversity indices in R and Rarefaction and accumulation curves, which is the honest way to compare richness across the uneven effort this post exposed. The lookup step is a join, covered in Reading field data into R.

Related tutorials

• Reading field data into R
• Wide and long: reshaping species data
• Diversity indices in R
• Rarefaction and accumulation curves

References

• Wickham H 2014 Journal of Statistical Software 59(10):1-23 (10.18637/jss.v059.i10)
• Boakes EH, Gliozzo G, Seymour V, Harvey M, Smith C, Roy DB, Haklay M 2016 Scientific Reports 6:33051 (10.1038/srep33051)
• Chao A, Gotelli NJ, Hsieh TC, Sander EL, Ma KH, Colwell RK, Ellison AM 2014 Ecological Monographs 84(1):45-67 (10.1890/13-0133.1)