Sea Monsters that Lost their Way

Graphic by Carl Goodwin
library(tidyverse)
library(tidymodels)
library(probably)
library(finetune)
library(textrecipes)
library(stopwords)
library(wesanderson)
library(kableExtra)
library(clock)
library(glue)
library(janitor)
library(vip)
library(ggrepel)
library(tictoc)
library(doParallel)

registerDoParallel(cores = 6)
theme_set(theme_bw())

(cols <- wes_palette(name = "Darjeeling2"))

The Natural History Museum began recording cetacean (whales, dolphins and porpoises) strandings in 1913. I’d like to explore this 1913-1989 dataset.

strandings_df <- read_csv("strandings.csv") %>%
  clean_names() %>% 
  mutate(
    date = date_parse(date, format = "%d/%m/%Y"),
    length = parse_number(length_et),
    species_lumped = fct_lump_n(species, 20),
    across(ends_with("_val"), as.integer)
  )

# glimpse(strandings_df)

Exploratory

Some of the species labels contain a question mark or forward slash. This indicates uncertainty, so it might be fun to see if a machine learning model (multiclass classification) could learn from the known species and suggest an appropriate species where it’s uncertain.

strandings_df2 <- 
  strandings_df %>% 
  mutate(species_uncertainty =
      if_else(str_detect(species, "[?/]"), "Uncertain", "Known"))

strandings_df2 %>% 
  filter(species_uncertainty == "Uncertain") %>% 
  count(species, sort = TRUE, name = "Count") %>% 
  slice_head(n = 6) %>% 
  kbl() %>% 
  kable_material()
species Count
delphis/coeruleoalba 48
phocoena? 42
melaena? 20
delphis? 18
truncatus? 18
acutorostrata? 12

The date variable has many NA’s. Fortunately, the components to construct many of these are in the year_val, month_val and day_val variables. With a little wrangling and imputation, we can coalesce these variables into a new date. This will be useful since plots of sample species by year, month and week of stranding suggest a de-constructed date could be a useful predictor.

strandings_df2 %>% 
  select(date, year_val, month_val, day_val) %>% 
  summary()
##       date               year_val      month_val         day_val     
##  Min.   :1913-02-13   Min.   :   0   Min.   : 0.000   Min.   : 0.00  
##  1st Qu.:1933-09-09   1st Qu.:1933   1st Qu.: 4.000   1st Qu.: 9.00  
##  Median :1954-04-13   Median :1955   Median : 7.000   Median :16.00  
##  Mean   :1955-01-08   Mean   :1954   Mean   : 6.766   Mean   :15.66  
##  3rd Qu.:1979-03-21   3rd Qu.:1979   3rd Qu.:10.000   3rd Qu.:22.00  
##  Max.   :1989-12-25   Max.   :1989   Max.   :12.000   Max.   :31.00  
##  NA's   :121
strandings_df3 <- strandings_df2 %>%
  group_by(species) %>% 
  mutate(
    month_val = if_else(month_val == 0, mean(month_val) %>% 
                          as.integer(), month_val),
    day_val = if_else(day_val == 0, mean(day_val) %>% as.integer(), day_val),
    day_val = if_else(day_val == 0, 1L, day_val),
    date2 = date_build(year_val, month_val, day_val, invalid = "NA"),
  ) %>% 
  ungroup() %>% 
  mutate(date3 = coalesce(date, date2)) %>% 
  arrange(id) %>% 
  mutate(date = if_else(is.na(date), lag(date3), date3)) %>% 
  select(-date2, -date3, -ends_with("_val"))

example_species <-
  c("musculus", "melas", "crassidens", "borealis", "coeruleoalba")

known_species <- strandings_df3 %>% 
  filter(species_uncertainty == "Known")

plot_date_feature <- function(var) {
  known_species %>%
    mutate(
      year = get_year(date),
      month = get_month(date),
      week = as_iso_year_week_day(date) %>% get_week()
    ) %>%
    filter(species %in% example_species) %>%
    count(species, {{ var }}) %>%
    ggplot(aes(species, {{ var }})) +
    geom_violin(
      alpha = 0.7,
      fill = cols[3],
      show.legend = FALSE
    ) +
    coord_flip() +
    labs(
      title = glue("Variation in {str_to_title(as.character(var))}",
                   " of Stranding for Known Species"),
      x = NULL, y = glue("{str_to_title(as.character(var))}")
    )
}

c("year", "month", "week") %>% 
  map(sym) %>% 
  map(plot_date_feature)
## [[1]]

## 
## [[2]]

## 
## [[3]]

Do latitude and longitude carry useful predictive information?

A geospatial visualisation of strandings shows some species do gravitate towards particular stretches of coastline, e.g. “acutus” and “albirostris” to the east, and “coeruleoalba” to the south-west.

Some species may also be more prone to mass stranding, so something that indicates whether a species has such a history (in these data) may be worth including in the mix.

uki <- map_data("world", region = c("uk", "ireland"))

labels <- c("Mass", "Single")

uki %>% 
  ggplot() +
  geom_map(aes(long, lat, map_id = region), map = uki, 
           colour = "black", fill = "grey90", size = 0.1) +
  geom_jitter(aes(longitude, latitude, colour = mass_single, 
                  size = mass_single), 
              alpha = 0.5, data = known_species) +
  facet_wrap(~ species_lumped, nrow = 3) +
  coord_map("mollweide") +
  scale_size_manual(values = c(1, 0.5), labels = labels) +
  scale_colour_manual(values = cols[c(3, 2)], labels = labels) +
  theme_void() +
  theme(legend.position = "top", 
        strip.text = element_text(colour = "grey50")) +
  labs(title = "Strandings by Species", 
       colour = NULL, size = NULL)

# Add history of mass stranding
strandings_df4 <- strandings_df3 %>% 
  group_by(species) %>%
  mutate(mass_possible = min(mass_single, na.rm = TRUE)) %>%
  ungroup()

Some records are missing the length measurement of the mammal. Nonetheless, where present, this is likely to be predictive, and may help, for example, separate species labelled as “delphis/coeruleoalba” where the length is at the extreme ends of the “delphis” range as we see below. And the range of length may differ by mammal sex.

known_species %>%
  mutate(sex = case_when(
    sex == "M" ~ "Male",
    sex == "F" ~ "Female",
    TRUE       ~ "Unknown"
  )) %>% 
  filter(species_lumped != "Other") %>% 
  count(species_lumped, length, sex) %>% 
  mutate(species_lumped = fct_reorder(species_lumped, 
                                      desc(length), min, na.rm = TRUE)) %>% 
  ggplot(aes(species_lumped, length)) + 
  geom_violin(aes(fill = if_else(str_detect(species_lumped, "^de|^co"), 
                                 TRUE, FALSE)), show.legend = FALSE) +
  facet_wrap(~ sex) +
  scale_fill_manual(values = cols[c(1, 5)]) +
  coord_flip() +
  labs(title = "Variation in Species Length by Sex", 
       x = NULL, y = "Length (metres)")

With map coordinates not always available, county could be, with a little string cleaning, a useful additional predictor.

strandings_df4 %>% 
  count(county) %>% 
  filter(str_detect(county, "Shet|Northumberland")) %>% 
  kbl() %>% 
  kable_material()
county n
Fair Isle, Shetland Isles 1
Northumberland 89
Northumberland. 1
Shetland Islands, Scotland 232
Shetland Isles, Scotland 35
Shetland, Scotland 1
Shetlands, Scotland 1
regex_pattern <-
  c("[,/].*$",
    "(?<!Che|Hamp|Lanca|North York)-?shire",
    " Isl.*$",
    " &.*$",
    "[-.]$")

strandings_df5 <- strandings_df4 %>%
  mutate(
    county = str_remove_all(county, str_c(regex_pattern, collapse = "|")),
    county = recode(
      county,
      "Carnarvon" = "Caernarvon",
      "E.Lothian" = "East Lothian",
      "Shetlands" = "Shetland",
      "W.Glamorgan" = "West Glamorgan",
      "S.Glamorgan" = "South Glamorgan"
    )
  ) 

strandings_df4 %>%
  summarise(counties_before = n_distinct(county))
## # A tibble: 1 × 1
##   counties_before
##             <int>
## 1             146
strandings_df5 %>%
  summarise(counties_after = n_distinct(county))
## # A tibble: 1 × 1
##   counties_after
##            <int>
## 1            109

Whilst count also appears to hold, based on the plot pattern below, species-related information, I’m not going to use it as a predictor as we do not know enough about how it was derived, as reflected in these variable descriptions.

strandings_df5 %>%
  ggplot(aes(species, count, colour = species_uncertainty)) +
  geom_jitter(alpha = 0.5, size = 2) +
  coord_flip() +
  scale_y_log10() +
  scale_colour_manual(values = cols[c(1, 5)]) +
  labs(title = "How 'Count' Relates to Species", 
       x = NULL, y = "Count (log10)", colour = "Species") +
  theme(legend.position = "top")

Modelling

So, I’ll set aside the rows where the species is uncertain (to be used later for new predictions), and I’ll train a model on 75% of known species, and test it on the remaining 25%. I’ll use the following predictors:

  • latitude and longitude
  • Mammal length and sex
  • mass_possible indicating a species history of mass strandings
  • date reported converted into decimal, week, month and year
  • county may be useful, especially where the longitude and latitude are missing
  • fam_genus which narrows the range of likely species

I’d like to also make use of the textrecipes package. I can tokenise the textual information in comment and location to see if these add to the predictive power of the model.

I’ll tune the model using tune_race_anova which quickly discards hyperparameter combinations showing little early promise.

known_species <- strandings_df5 %>%
  filter(species_uncertainty == "Known") %>%
  mutate(across(
    c(
      "species",
      "mass_single",
      "mass_possible",
      "county",
      "location",
      "sex",
      "fam_genus"
    ),
    factor
  ))

set.seed(123)

data_split <-
  known_species %>%
  mutate(species = fct_drop(species)) %>% 
  initial_split(strata = species)

train <- data_split %>% training()

test <- data_split %>% testing()

predictors <-
  c(
    "latitude",
    "longitude",
    "length",
    "mass_single",
    "mass_possible",
    "county",
    "location",
    "comment",
    "sex",
    "fam_genus"
  )

recipe <-
  train %>%
  recipe() %>%
  update_role(species, new_role = "outcome") %>%
  update_role(all_of(predictors), new_role = "predictor") %>%
  update_role(!has_role("outcome") & !has_role("predictor"), 
              new_role = "id") %>%
  step_date(date, features = c("decimal", "week", "month", "year"), 
            label = FALSE) %>%
  step_tokenize(location, comment) %>%
  step_stopwords(location, comment) %>%
  step_tokenfilter(location, comment, max_tokens = tune()) %>% #100
  step_tf(location, comment) %>%
  step_zv(all_predictors()) %>%
  step_dummy(all_nominal_predictors())

xgb_model <-
  boost_tree(trees = tune(), # 440
             mtry = tune(), # 0.6
             learn_rate = 0.02) %>% 
  set_mode("classification") %>%
  set_engine("xgboost", 
             count = FALSE,
             verbosity = 0,
             tree_method = "hist")

xgb_wflow <- workflow() %>%
  add_recipe(recipe) %>%
  add_model(xgb_model)

set.seed(9)

folds <- vfold_cv(train, strata = species)

set.seed(10)

tic()

tune_result <- xgb_wflow %>%
  tune_race_anova(
    resamples = folds,
    grid = crossing(
      trees = seq(200, 520, 40),
      mtry = seq(0.5, 0.7, 0.1),
      max_tokens = seq(80, 120, 20)
      ),
    control = control_race(),
    metrics = metric_set(accuracy)
  )

toc()
## 438.622 sec elapsed
tune_result %>% 
  plot_race() + 
  labs(title = "Early Elimination of Parameter Combinations")

set.seed(123)

xgb_fit <- xgb_wflow %>%
  finalize_workflow(tune_result %>% 
                      select_best(metric = "accuracy")) %>% 
  fit(train)

Having fitted the model with the 3080 records in the training data, I’ll test its accuracy on the 1028 records of known species the model has not yet seen.

Without spending time on alternative models, we’re getting a reasonable result for the “porpoise” of this post, as reflected in both the accuracy metric and confusion matrix.

xgb_results <- xgb_fit %>% 
  augment(new_data = test)

xgb_results %>%
  accuracy(species, .pred_class)
## # A tibble: 1 × 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy multiclass     0.993
xgb_results %>%
  conf_mat(species, .pred_class) %>%
  autoplot(type = "heatmap") +
  scale_fill_gradient2(
    mid = "white",
    high = cols[1],
    midpoint = 0
  ) +
  labs(title = "Confusion Matrix") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

The top variable importance scores include fam_genus, many of the comment tokens, plus length, mass-possible, date_decimal, date_year, and latitude.

vi(xgb_fit %>% extract_fit_parsnip()) %>% 
  arrange(desc(Importance)) %>% 
  mutate(ranking = row_number()) %>% 
  slice_head(n = 40) %>% 
  kbl() %>% 
  kable_material()
Variable Importance ranking
fam_genus_Globicephala 0.1170072 1
fam_genus_Phocoena 0.0944623 2
tf_comment_unidentified 0.0759441 3
fam_genus_Delphinus 0.0711414 4
tf_comment_false 0.0504881 5
fam_genus_Tursiops 0.0465281 6
fam_genus_Lagenorhynchus 0.0443353 7
tf_comment_porpoise 0.0403505 8
tf_comment_sided 0.0285027 9
fam_genus_Grampus 0.0276863 10
tf_comment_lesser 0.0275691 11
length 0.0265941 12
fam_genus_Hyperoodon 0.0232933 13
tf_comment_beaked 0.0230283 14
fam_genus_Orcinus 0.0194158 15
fam_genus_Pseudorca 0.0186163 16
tf_comment_dolphin 0.0182272 17
tf_comment_bottle 0.0169981 18
fam_genus_Physeter 0.0158601 19
tf_comment_rorqual 0.0138746 20
fam_genus_Ziphius 0.0131200 21
mass_possible_S 0.0131200 22
fam_genus_Mesoplodon 0.0125859 23
fam_genus_odontocete 0.0107553 24
tf_comment_fin 0.0106025 25
tf_comment_long 0.0102623 26
date_decimal 0.0100120 27
tf_comment_finned 0.0093989 28
tf_comment_nosed 0.0093516 29
tf_comment_pilot 0.0093481 30
tf_comment_common 0.0089965 31
fam_genus_Stenella 0.0085796 32
tf_comment_whale 0.0084079 33
tf_comment_risso’s 0.0082206 34
date_year 0.0072894 35
tf_comment_killer 0.0070108 36
tf_comment_white 0.0067963 37
fam_genus_cetacean 0.0060907 38
tf_comment_sowerby’s 0.0039810 39
latitude 0.0030011 40

Do the predictions look reasonable?

The class probability is spread across 27 species. I’m going to set a high threshold of 0.9, meaning the predicted species needs to be a pretty confident prediction.

xgb_preds <- xgb_fit %>% 
  augment(new_data = strandings_df5 %>% 
            filter(species_uncertainty == "Uncertain"))

species_levels <- xgb_preds %>% 
  select(starts_with(".pred_"), -.pred_class) %>% 
  names() %>% 
  as.factor()

subset_df <- xgb_preds %>%
  mutate(
    .class_pred = make_class_pred(
      .pred_acutorostrata,
      .pred_acutus,
      .pred_albirostris,
      .pred_ampullatus,
      .pred_bidens,
      .pred_borealis,
      .pred_breviceps,
      .pred_cavirostris,
      .pred_coeruleoalba,
      .pred_crassidens,
      .pred_delphis,
      .pred_electra,
      .pred_europaeus,
      .pred_griseus,
      .pred_leucas,
      .pred_macrocephalus,
      .pred_melaena,
      .pred_melas,
      .pred_mirus,
      .pred_monoceros,
      .pred_musculus,
      .pred_novaeangliae,
      .pred_orca,
      .pred_phocoena,
      .pred_physalus,
      .pred_sp.indet.,
      .pred_truncatus,
      levels = species_levels,
      min_prob = .9
    )
  )

subset_df %>%
  group_by(species, .class_pred) %>% 
  summarise(n = n()) %>% 
  arrange(species, desc(n)) %>% 
  kbl() %>% 
  kable_material()
species .class_pred n
acutorostrata? .pred_acutorostrata 12
acutorostrata/borealis .pred_acutorostrata 1
acutus? .pred_acutus 3
albirostris? .pred_albirostris 9
ampullatus? .pred_ampullatus 3
bidens? .pred_bidens 2
bidens? [EQ] 1
cavirostris? .pred_cavirostris 7
coeruleoalba? .pred_coeruleoalba 1
delphis? .pred_delphis 18
delphis/coeruleoalba [EQ] 48
griseus? .pred_griseus 2
macrocephalus? .pred_macrocephalus 2
melaena? .pred_melaena 20
orca? .pred_orca 4
phocoena? .pred_phocoena 42
physalus? .pred_physalus 4
physalus/acutorostrata [EQ] 1
truncatus? .pred_truncatus 18
truncatus/albirostris [EQ] 5

The majority of the 203 uncertain records are predicted to be as suspected in the original labelling. The remainder are classed as equivocal as they have not met the high bar of a 0.9-or-above probability for a single species.

R Toolbox

Summarising below the packages and functions used in this post enables me to separately create a toolbox visualisation summarising the usage of packages and functions across all posts.

Package Function
base as.character[2]; as.integer[3]; c[9]; comment[3]; conflicts[1]; cumsum[1]; function[2]; is.na[1]; length[2]; log10[1]; mean[2]; min[1]; names[1]; search[1]; seq[3]; set.seed[4]; sum[1]; summary[1]
clock as_iso_year_week_day[1]; date_build[1]; date_parse[1]; get_month[1]; get_week[1]; get_year[1]
doParallel registerDoParallel[1]
dplyr filter[12]; across[2]; arrange[5]; case_when[1]; coalesce[1]; count[4]; desc[5]; group_by[4]; id[1]; if_else[9]; mutate[18]; n[2]; n_distinct[2]; recode[1]; row_number[1]; select[3]; slice_head[2]; summarise[4]; ungroup[2]
finetune control_race[1]; plot_race[1]; tune_race_anova[1]
forcats fct_drop[1]; fct_lump_n[1]; fct_reorder[1]
ggplot2 aes[6]; coord_flip[3]; coord_map[1]; element_text[2]; facet_wrap[2]; geom_jitter[2]; geom_map[1]; geom_violin[2]; ggplot[4]; labs[6]; map_data[1]; scale_colour_manual[2]; scale_fill_gradient2[1]; scale_fill_manual[1]; scale_size_manual[1]; scale_y_log10[1]; theme[3]; theme_bw[1]; theme_set[1]; theme_void[1]
glue glue[2]
janitor clean_names[1]
kableExtra kable_material[5]; kbl[5]
parsnip boost_tree[1]; set_engine[1]; set_mode[1]
probably make_class_pred[1]
purrr map[3]; map2_dfr[1]; possibly[1]
readr parse_number[1]; read_csv[1]; read_lines[1]
recipes all_nominal_predictors[1]; all_predictors[1]; has_role[2]; step_date[1]; step_dummy[1]; step_zv[1]; update_role[3]
rsample initial_split[1]; testing[1]; training[1]; vfold_cv[1]
stats var[3]
stringr str_c[6]; str_count[1]; str_detect[5]; str_remove[2]; str_remove_all[2]; str_starts[1]; str_to_title[2]
textrecipes step_stopwords[1]; step_tf[1]; step_tokenfilter[1]; step_tokenize[1]
tibble as_tibble[1]; tibble[2]; enframe[1]
tictoc tic[1]; toc[1]
tidyr crossing[1]; unnest[1]
tune finalize_workflow[1]; select_best[1]; tune[3]
wesanderson wes_palette[1]
workflows add_model[1]; add_recipe[1]; workflow[1]
yardstick accuracy[2]; conf_mat[1]; metric_set[1]

Citation

Natural History Museum (2019). Data Portal Query on “UK cetacean strandings 1913-1989” created at 2019-08-10 16:41:12.475340 PID https://doi.org/10.5519/qd.iwg63595. Subset of “Historical UK cetacean strandings dataset (1913-1989)” (dataset) PID https://doi.org/10.5519/0028204.

Carl Goodwin
Carl Goodwin
Data Scientist
comments powered by Disqus

Related