Sea Monsters that Lost their Way

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(tictoc)
library(doParallel)

registerDoParallel(cores = 6)

theme_set(theme_bw())

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

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)

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()

strandings_df2 |> 
  select(date, year_val, month_val, day_val) |> 
  summary()

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)

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()

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)")

strandings_df4 |> 
  count(county) |> 
  filter(str_detect(county, "Shet|Northumberland")) |> 
  kbl(col.names = c("County", "Count"))

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))

strandings_df5 |>
  summarise(counties_after = n_distinct(county))

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")

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()

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)

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

xgb_results |>
  accuracy(species, .pred_class)

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))

vi(xgb_fit |> extract_fit_parsnip()) |> 
  arrange(desc(Importance)) |> 
  mutate(ranking = row_number()) |> 
  slice_head(n = 40) |> 
  kbl()

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(col.names = c("Actual Label", "Predicted Label", "Count"))