OP 7 & 8 Data Package 1: Baseline data collection and model development

Download as PDF

1 Data package overview

DRAFT REPORT V 0.0.1

This report forms contribution to OP7 and OP8. Data package 1 contains raw validation videos, the processed sensor datasets, the Python model code and outputs, and the comparative blade strike estimates. At this stage, no inlet/hub guidance has been incorperated. Data package 1 therefore forms a baseline for development throughout the KTP. The folder structure is:

Data package 1 contents
Folder	Contents
Video library	348 validation passage videos (Tests 1-9, 400 rpm), organised by test and run
Sensor data library	Shortcut to the external sensor data library (raw sensor files)
HIFI_BSM- python code and model datasets	Python model training scripts and the processed model datasets (python_data)
Model1_1_binary	Binary blade strike model outputs: cross-validated predictions, performance metrics, diagnostic figures and the deployable model
Model1_1_multiclass	Multiclass (collapsed concentric region) model outputs: predictions, performance metrics and diagnostic figures
BSM - FS700-N400Q09	CEN (semi-mechanistic) blade strike model code and results for the FS700 pump at 400 rpm, Q0.9
Data-driven predictions_alldata	Data-driven model predictions applied across the full dataset, with per-file and per-treatment summaries

The core dataset (Pumpflow_2024_and_2026_datatset_.csv) combines the 2024 and 2026 PumpFlow sensor deployments. Each row is a time step of the sensor signal at 2000hz with a window size of 200ms. Each sensor file corresponds to one passage through the pump and is annotated using the video ground truth (passage type, leading edge strike type, other collision type and concentric pump region).

2 Data import

All data used in the report are imported here.

# Combined 2024 + 2026 sensor dataset with video ground truth annotations
data <- read_csv("./HIFI_BSM- python code and model datasets/python_data/Pumpflow_2024_and_2026_datatset_.csv") %>%
  mutate(
    treatment = factor(
      treatment,
      levels = c("500 (100%)", "400 (70%)", "400_model", "400 (100%)"),
      labels = c("500 (100%)", "400 (70%)", "400 (80%)", "400 (100%)")
    ),
    passage_type = as.factor(passage_type),
    leading_type = as.factor(leading_type),
    other_type = as.factor(other_type)
  )

# Copy of the global sensor index from the data library
global_sensor_dataset2 <- read_csv("./HIFI_BSM- python code and model datasets/python_data/global_sensor_dataset2.csv")

# Existing video metadata file
uoh_meta <- read_csv("./HIFI_BSM- python code and model datasets/python_data/UOH_PF_KTP_datatset2_meta.csv")

# Binary model out-of-fold (OOF) blade strike estimates by treatment
blade_strike_predictions <- read_csv("./Model1_1_binary/blade_strike_predictions.csv")

# Binary model cross-validated per-file predictions
cv_predictions <- read_csv("./Model1_1_binary/cv_predictions.csv")

# Reviewed misclassified files with correction outcomes
misclass1_done <- read_csv("./HIFI_BSM- python code and model datasets/python_data/misclass1_DONE.csv")

# CEN model blade strike estimates (FS700, 400 rpm, Q0.9)
pf_q09_results <- read_csv("./BSM - FS700-N400Q09/pf_q09_results.csv")

# Multiclass model: predicted vs ground truth concentric region by treatment
class_by_treatment <- read_csv("./Model1_1_multiclass/class_by_treatment.csv")

# Data-driven predictions across the full dataset, summarised by treatment
prediction_summary <- read_csv("./Data-driven predictions_alldata/prediction_summary.csv")

# Binary model performance metrics
bin_metrics <- fromJSON("./Model1_1_binary/performance_metrics.json")

3 Data preparation

All cleaning, labelling and assembly happens here. The sensor data are summarised to one row per file, keeping the annotation columns and the maximum high-g acceleration magnitude (max_accmag) as a proxy for strike severity. Treatments are relabelled consistently across the model output files (400_model = 400 (80%), the model 1.1 dataset).

treatment_levels <- c("500 (100%)", "400 (70%)", "400 (80%)", "400 (100%)")

relabel_treatment <- function(x) {
  factor(recode(x, "400_model" = "400 (80%)"), levels = treatment_levels)
}

drop_cols <- c(
  "time_s","higacc_x_g","higacc_y_g","higacc_z_g","higacc_mag_g",
  "inacc_x_ms","inacc_y_ms","inacc_z_ms","inacc_mag_ms",
  "rot_x_degs","rot_y_degs","rot_z_degs","rot_mag_degs","pressure_kpa"
)

blade_strike_data <- data %>%
  group_by(file) %>%
  summarise(
    across(-higacc_mag_g, first),
    max_accmag = max(higacc_mag_g, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(
    concentric_two = case_when(
      is.na(concentric_pump_region) ~ NA_character_,
      concentric_pump_region == 1 ~ "C1",
      concentric_pump_region %in% c(2, 3) ~ "C2_3",
      concentric_pump_region %in% c(4, 5) ~ "C4_5"
    ),
    treatment2 = case_when(
      treatment == "400 (80%)" ~ "Model 1.1",
      TRUE ~ "Model 1"
    ),
    strike = case_when(
      passage_type != "No contact" ~ 1L,
      other_type == "Impeller surface" ~ 0L,
      TRUE ~ 0L
    )
  ) %>%
  select(-any_of(drop_cols))

# Apply the same treatment relabelling to the model output files
blade_strike_predictions <- blade_strike_predictions %>%
  mutate(treatment = relabel_treatment(treatment))

class_by_treatment <- class_by_treatment %>%
  mutate(treatment = relabel_treatment(treatment))

prediction_summary <- prediction_summary %>%
  mutate(treatment = relabel_treatment(treatment))

3.1 model 1.1 metadata assembly

The model 1.1 (2026) blade strike records are matched to the global sensor index and video metadata, then exported as the final metadataset for the package. This can be expanded to the whole dataset by substituting blade_strike_data for the model 1.1 subset.

# Step 1 keep only files present in both global_sensor_dataset2 and
# blade_strike_data (model 1.1)
bsm2 <- blade_strike_data %>% filter(treatment2 == "Model 1.1")

common_files <- intersect(bsm2$file, global_sensor_dataset2$file)

bsm2 <- bsm2 %>% filter(file %in% common_files)
global_matched <- global_sensor_dataset2 %>% filter(file %in% common_files)

# Step 2 combine sensor columns with video metadata (meta keyed on sens_file)
sensor_meta <- global_matched %>%
  select(file, sensor, `duration.mm.ss.`, deployment_id, pump_turbine,
         type, head, point_bep) %>%
  left_join(
    uoh_meta %>% select(sens_file, date, video_file, pump_n, pump_q),
    by = c("file" = "sens_file")
  )

# Step 3 append to the model 1.1 blade strike data
model2_data <- bsm2 %>%
  left_join(sensor_meta, by = "file") %>%
  select(
    file, sensor, date, `duration.mm.ss.`, video_file, deployment_id,
    pump_turbine, type, head, pump_n, pump_q, point_bep, strike,
    passage_type, leading_type, other_type, concentric_pump_region,
    concentric_two, max_accmag
  )

# Export and save final metadataset
write_csv(model2_data, "./HIFI_BSM- python code and model datasets/UoH_AA_PUMPFLOW_KTP_DATA_PACKAGE1.csv")

The combined dataset holds 565 sensor files (307 in the model 1.1 subset).

4 Ground truth annotation summary

Distribution of the video ground truth annotations across the model 1.1 dataset.

annotation_long <- model2_data %>%
  transmute(
    `Passage type`      = as.character(passage_type),
    `Leading type`      = as.character(leading_type),
    `Other type`        = as.character(other_type),
    `Concentric region` = as.character(concentric_pump_region)
  ) %>%
  pivot_longer(everything(), names_to = "panel", values_to = "level") %>%
  filter(!is.na(level), level != "None") %>%
  count(panel, level, name = "Sensors") %>%
  group_by(panel) %>%
  mutate(`%` = round(100 * Sensors / sum(Sensors), 1)) %>%
  ungroup() %>%
  mutate(panel = factor(panel, levels = c("Passage type", "Leading type",
                                          "Other type", "Concentric region")))

annotation_tbl <- annotation_long %>% arrange(panel, desc(Sensors))

annotation_tbl %>%
  gt(groupname_col = "panel") %>%
  tab_header(title = "Annotation summary (ground truth)") %>%
  gt_style(save_as = "annotation_summary")

Annotation summary (ground truth)
level	Sensors	%
Passage type
No contact	161	52.4
Leading edge strike	107	34.9
Other impeller collision	39	12.7
Leading type
Direct	98	91.6
Indirect	9	8.4
Other type
Impeller hub	29	74.4
Impeller surface	10	25.6
Concentric region
4	40	29.4
2	35	25.7
1	29	21.3
3	19	14.0
5	13	9.6

ggplot(annotation_long, aes(level, Sensors)) +
  geom_col(fill = "white", colour = "black", width = 0.6) +
  geom_text(aes(label = Sensors), vjust = -0.4, size = 3, fontface = "bold") +
  facet_wrap(~ panel, ncol = 2, scales = "free", strip.position = "bottom") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(x = NULL, y = "Sensors (n)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Ground truth annotation counts for the model 1.1 dataset.

5 Blade strike ground truth summary

Video ground truth blade strike counts and rates by treatment, with Wilson 95% confidence intervals. Concentric region counts (C1, C2-3, C4-5) are given.

blade_tbl1 <- blade_strike_data %>%
  group_by(treatment) %>%
  summarise(
    N = n(),
    Strikes = sum(strike, na.rm = TRUE),
    C1 = sum(concentric_pump_region == 1, na.rm = TRUE),
    C2_3 = sum(concentric_pump_region %in% c(2, 3), na.rm = TRUE),
    C4_5 = sum(concentric_pump_region %in% c(4, 5), na.rm = TRUE),
    `No strike` = sum(strike == 0, na.rm = TRUE),
    `Strike rate` = round(100 * Strikes / N, 1),
    .groups = "drop"
  ) %>%
  mutate(
    ci = map2(Strikes, N, ~ binom.confint(.x, .y, method = "wilson")),
    `95% CI` = map_chr(ci, ~ paste0(
      round(100 * .x$lower, 1),
      "–",
      round(100 * .x$upper, 1)
    ))
  ) %>%
  select(-ci)

blade_tbl1 %>%
  gt() %>%
  tab_header(title = "Blade strike summary (video ground truth)") %>%
  gt_style(save_as = "blade_strike_summary")

Blade strike summary (video ground truth)
treatment	N	Strikes	C1	C2_3	C4_5	No strike	Strike rate	95% CI
500 (100%)	76	38	7	12	16	38	50.0	39–61
400 (70%)	80	40	6	18	16	40	50.0	39.3–60.7
400 (80%)	307	146	29	54	53	161	47.6	42–53.1
400 (100%)	102	40	5	12	19	62	39.2	30.3–48.9

plot_dat1 <- blade_strike_data %>%
  group_by(treatment) %>%
  summarise(
    N = n(),
    strikes = sum(strike, na.rm = TRUE),
    .groups = "drop"
  )

plot_dat1 <- plot_dat1 %>%
  bind_cols(
    binom.confint(plot_dat1$strikes, plot_dat1$N, method = "wilson")[, c("lower", "upper")]
  ) %>%
  mutate(rate = 100 * strikes / N)

ggplot(plot_dat1, aes(treatment, rate)) +
  geom_col(fill = "white", colour = "black") +
  geom_errorbar(aes(ymin = lower * 100, ymax = upper * 100), width = 0.2) +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  labs(x = NULL, y = "Strike rate (%)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Video ground truth strike rate by treatment (Wilson 95% CI).

6 Blade strike estimate comparison

This section is how blade strike estimates are produced and compared for a given operating point. The example compares, for the model 1.1 dataset (400 rpm, Q0.9):

• Video ground truth from blade_tbl1 above;
• Model 1.1 (OOF) the out-of-fold blade strike estimate from blade_strike_predictions.csv, matched by treatment;
• CEN the semi-mechanistic estimate from pf_q09_results.csv;
• Custom a user-supplied x strikes out of n fish with a custom label.

Comment or uncomment the blocks below to control which estimates are shown; the figure displays whichever groups it is given.

estimates <- list(

  #Video ground truth
  video = blade_tbl1 %>%
    filter(treatment == "400 (80%)") %>%
    transmute(group = "Video ground truth", x = Strikes, n = N),

  # Model 1.1 out-of-fold estimate (matched by treatment)
  model_oof = blade_strike_predictions %>%
    filter(treatment == "400 (80%)") %>%
    transmute(group = "Model 1.1 (OOF)", x = n_predicted_strike, n = n_fish),

  # Custom estimate (user supplied x strikes out of n fish)
  # custom = tibble(group = "Custom", x = 25, n = 100),

  NULL
)

comparison_dat <- bind_rows(estimates)
ci <- binom.confint(comparison_dat$x, comparison_dat$n, method = "wilson")
comparison_dat <- comparison_dat %>%
  mutate(rate = 100 * x / n, lower = 100 * ci$lower, upper = 100 * ci$upper) %>%
  select(group, rate, lower, upper)

# CEN estimate (no CI available) comment out to drop from the plot
comparison_dat <- bind_rows(
  comparison_dat,
  pf_q09_results %>%
    filter(method == "CEN") %>%
    transmute(group = "CEN", rate = Pco_percent,
              lower = NA_real_, upper = NA_real_)
)

comparison_dat <- comparison_dat %>% mutate(group = fct_inorder(group))

ggplot(comparison_dat, aes(group, rate)) +
  geom_col(fill = "white", colour = "black", width = 0.6) +
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0.2, na.rm = TRUE) +
  geom_text(aes(y = coalesce(upper, rate), label = sprintf("%.1f%%", rate)),
            vjust = -0.4, size = 3, fontface = "bold") +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  labs(x = NULL, y = "Collision probability (%)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Collision probability estimates for the 400 (80%) treatment (Wilson 95% CI where available).

7 Strike severity: maximum acceleration

The maximum high-g acceleration magnitude per passage (max_accmag) is used as a proxy for strike severity.

7.1 By dataset and passage type

The model 1.1 sensors had rubber sleeves. Were strikes lower severity in sleeves?

tbl5 <- blade_strike_data %>%
  group_by(treatment2, passage_type) %>%
  summarise(
    Median = median(max_accmag, na.rm = TRUE),
    Min = min(max_accmag, na.rm = TRUE),
    Max = max(max_accmag, na.rm = TRUE),
    IQR = IQR(max_accmag, na.rm = TRUE),
    .groups = "drop"
  )

tbl5 %>%
  gt() %>%
  tab_header(title = "Maximum acceleration (g) by dataset and passage type") %>%
  fmt_number(columns = c(Median, Min, Max, IQR), decimals = 1) %>%
  gt_style(save_as = "maxaccmag_by_dataset_passage")

Maximum acceleration (g) by dataset and passage type
treatment2	passage_type	Median	Min	Max	IQR
Model 1	Leading edge strike	513.2	195.4	692.4	132.3
Model 1	No contact	18.0	0.0	566.4	22.7
Model 1	Other impeller collision	165.9	13.7	570.4	372.7
Model 1.1	Leading edge strike	420.7	15.5	600.5	233.1
Model 1.1	No contact	15.0	0.0	403.5	13.2
Model 1.1	Other impeller collision	47.2	17.9	566.6	60.0

ggplot(blade_strike_data,
       aes(interaction(treatment2, passage_type), max_accmag)) +
  geom_boxplot(colour = "black") +
  labs(x = NULL, y = "Maximum acceleration (g)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Maximum acceleration by dataset and passage type.

7.2 Contact events only

tbl7 <- blade_strike_data %>%
  filter(passage_type != "No contact") %>%
  group_by(treatment2) %>%
  summarise(
    Median = median(max_accmag, na.rm = TRUE),
    Min = min(max_accmag, na.rm = TRUE),
    Max = max(max_accmag, na.rm = TRUE),
    IQR = IQR(max_accmag, na.rm = TRUE),
    .groups = "drop"
  )

tbl7 %>%
  gt() %>%
  tab_header(title = "Maximum acceleration (g), contact events only") %>%
  fmt_number(columns = c(Median, Min, Max, IQR), decimals = 1) %>%
  gt_style(save_as = "maxaccmag_contact_only")

Maximum acceleration (g), contact events only
treatment2	Median	Min	Max	IQR
Model 1	490.4	13.7	692.4	159.0
Model 1.1	398.1	15.5	600.5	396.2

ggplot(
  filter(blade_strike_data, passage_type == "Leading edge strike"),
  aes(treatment2, max_accmag)
) +
  geom_boxplot(colour = "black") +
  labs(x = NULL, y = "Maximum acceleration (g)")

Maximum acceleration for leading edge strikes, Model 1 vs model 1.1.

A refinement for future iterations would be to isolate the older data by excluding the 500 (100%) treatment from the Model 1 group.

7.3 By concentric region

NA cases are retained as “No contact” passages.

tbl9 <- blade_strike_data %>%
  group_by(concentric_pump_region) %>%
  summarise(
    Median = median(max_accmag, na.rm = TRUE),
    Min = min(max_accmag, na.rm = TRUE),
    Max = max(max_accmag, na.rm = TRUE),
    IQR = IQR(max_accmag, na.rm = TRUE),
    .groups = "drop"
  )

tbl9 %>%
  gt() %>%
  tab_header(title = "Maximum acceleration (g) by concentric region") %>%
  fmt_number(columns = c(Median, Min, Max, IQR), decimals = 1) %>%
  gt_style(save_as = "maxaccmag_by_concentric_region")

Maximum acceleration (g) by concentric region
concentric_pump_region	Median	Min	Max	IQR
1	70.7	15.0	570.4	357.3
2	409.8	33.6	599.6	284.0
3	467.2	42.1	616.9	165.6
4	465.0	15.5	692.4	153.9
5	533.6	316.6	692.4	130.9
NA	16.9	0.0	567.1	18.5

ggplot(
  blade_strike_data %>%
    mutate(
      concentric_pump_region = fct_na_value_to_level(
        as.factor(concentric_pump_region),
        level = "No contact"
      ),
      concentric_pump_region = fct_relevel(
        concentric_pump_region,
        "No contact",
        "1", "2", "3", "4", "5"
      )
    ),
  aes(concentric_pump_region, max_accmag)
) +
  geom_jitter(
    width = 0.15,
    alpha = 0.3,
    size = 1,
    colour = "darkblue"
  ) +
  geom_boxplot(
    colour = "black",
    fill = NA,
    width = 0.45,
    outlier.shape = 4
  ) +
  scale_y_continuous(breaks = seq(0, 700, 100)) +
  labs(
    x = "Concentric region",
    y = "Maximum acceleration-mag (g)"
  )

Maximum acceleration by concentric pump region.

df <- blade_strike_data %>%
  filter(!is.na(concentric_pump_region), !is.na(max_accmag))

res <- cor.test(
  as.numeric(df$concentric_pump_region),
  df$max_accmag,
  method = "spearman",
  exact = FALSE
)
res

## 
##  Spearman's rank correlation rho
## 
## data:  as.numeric(df$concentric_pump_region) and df$max_accmag
## S = 1278726, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.4908507

Maximum acceleration increases with concentric region (Spearman’s rho = 0.49, p < 0.001). This gradient indicates the signal data carry a meaningful structure.

7.4 Collapsed concentric regions

tbl11 <- blade_strike_data %>%
  group_by(concentric_two) %>%
  summarise(
    Median = median(max_accmag, na.rm = TRUE),
    Min = min(max_accmag, na.rm = TRUE),
    Max = max(max_accmag, na.rm = TRUE),
    IQR = IQR(max_accmag, na.rm = TRUE),
    .groups = "drop"
  )

tbl11 %>%
  gt() %>%
  tab_header(title = "Maximum acceleration (g) by collapsed concentric region") %>%
  fmt_number(columns = c(Median, Min, Max, IQR), decimals = 1) %>%
  gt_style(save_as = "maxaccmag_concentric_collapsed")

Maximum acceleration (g) by collapsed concentric region
concentric_two	Median	Min	Max	IQR
C1	70.7	15.0	570.4	357.3
C2_3	427.6	33.6	616.9	275.2
C4_5	481.6	15.5	692.4	149.8
NA	16.9	0.0	567.1	18.5

ggplot(
  blade_strike_data %>%
    mutate(
      concentric_two = fct_na_value_to_level(
        as.factor(concentric_two),
        level = "No contact"
      ),
      concentric_two = fct_relevel(
        concentric_two,
        "No contact",
        "C1", "C2_3", "C4_5"
      )
    ),
  aes(concentric_two, max_accmag)
) +
  geom_boxplot(
    colour = "black",
    width = 0.45,
    outlier.shape = 4
  ) +
  geom_jitter(
    width = 0.15,
    alpha = 0.5,
    size = 1,
    colour = "darkblue"
  ) +
  labs(
    x = "Concentric region",
    y = "Maximum acceleration-mag (g)"
  )

Maximum acceleration by collapsed concentric region.

The separation between the region groups suggests the collapsed model will perform best in the model training scenario. again, this text needs expanding with justification

8 Concentric region distribution by treatment

Concentric pump region

Ground truth distribution of strike location (collapsed concentric regions) across treatments, from class_by_treatment.csv.

region_labels <- c(region_1 = "C1", region_2_3 = "C2_3", region_4_5 = "C4_5")

class_gt_long <- class_by_treatment %>%
  select(treatment, n_true_region_1, n_true_region_2_3, n_true_region_4_5) %>%
  pivot_longer(-treatment, names_prefix = "n_true_",
               names_to = "region", values_to = "n") %>%
  mutate(region = region_labels[region]) %>%
  group_by(treatment) %>%
  mutate(pct = 100 * n / sum(n)) %>%
  ungroup()

ggplot(class_gt_long, aes(treatment, pct, fill = region)) +
  geom_col(colour = "black", width = 0.6) +
  geom_text(aes(label = n), position = position_stack(vjust = 0.5), size = 3) +
  scale_fill_grey(start = 0.95, end = 0.4, name = "Region") +
  scale_y_continuous(limits = c(0, 100.001), expand = c(0, 0)) +
  labs(x = NULL, y = "Strikes (%)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Ground truth strike location composition by treatment.

9 Model predictions vs ground truth (model 1.1)

Left panel: out-of-fold model strike rate against the video ground truth for the model 1.1 dataset. Right panel: predicted vs ground truth strike location composition (collapsed concentric regions).

m2_pred <- blade_strike_predictions %>% filter(treatment == "400 (80%)")

strike_cmp <- tibble(
  source = c("Ground truth", "Model 1.1 (OOF)"),
  rate = 100 * c(m2_pred$n_true_strike / m2_pred$n_fish,
                 m2_pred$predicted_strike_rate),
  panel = "Strike rate"
)

region_cmp <- class_by_treatment %>%
  filter(treatment == "400 (80%)") %>%
  select(starts_with("n_pred_region"), starts_with("n_true_region")) %>%
  pivot_longer(
    everything(),
    names_pattern = "n_(pred|true)_(region_.*)",
    names_to = c("source", "region"),
    values_to = "n"
  ) %>%
  mutate(
    source = ifelse(source == "pred", "Model 1.1 (OOF)", "Ground truth"),
    region = region_labels[region],
    panel = "Concentric region composition"
  ) %>%
  group_by(source) %>%
  mutate(pct = 100 * n / sum(n)) %>%
  ungroup()

ggplot() +
  geom_col(data = strike_cmp, aes(source, rate),
           fill = "white", colour = "black", width = 0.6) +
  geom_text(data = strike_cmp,
            aes(source, rate, label = sprintf("%.1f%%", rate)),
            vjust = -0.4, size = 3, fontface = "bold") +
  geom_col(data = region_cmp, aes(source, pct, fill = region),
           colour = "black", width = 0.6) +
  geom_text(data = region_cmp, aes(source, pct, group = region, label = n),
            position = position_stack(vjust = 0.5), size = 3) +
  facet_wrap(~ factor(panel, levels = c("Strike rate",
                                        "Concentric region composition"))) +
  scale_fill_grey(start = 0.95, end = 0.4, name = "Region") +
  scale_y_continuous(limits = c(0, 100.001), expand = c(0, 0)) +
  labs(x = NULL, y = "Percent (%)")

Model 1.1 out-of-fold predictions vs video ground truth (model 1.1 dataset).

10 Predicted strike rates with region composition

Data-driven model predictions applied across the full dataset (prediction_summary.csv). Each bar is the predicted strike rate for a treatment, subdivided by the predicted strike location; error bars are the Wilson 95% CI on the strike rate.

strike_region_dat <- prediction_summary %>%
  pivot_longer(
    c(n_region_1, n_region_2_3, n_region_4_5),
    names_prefix = "n_",
    names_to = "region", values_to = "n_r"
  ) %>%
  mutate(
    region = region_labels[region],
    segment = 100 * n_r / n
  )

ggplot(strike_region_dat, aes(treatment, segment, fill = region)) +
  geom_col(colour = "black", width = 0.6) +
  geom_errorbar(
    data = prediction_summary,
    aes(treatment, ymin = 100 * ci_lo, ymax = 100 * ci_hi),
    inherit.aes = FALSE, width = 0.2
  ) +
  geom_text(
    data = prediction_summary,
    aes(treatment, y = 100 * ci_hi, label = sprintf("%.1f%%", 100 * strike_rate)),
    inherit.aes = FALSE, vjust = -0.4, size = 3, fontface = "bold"
  ) +
  scale_fill_grey(start = 0.95, end = 0.4, name = "Region") +
  scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10), expand = c(0, 0)) +
  labs(x = NULL, y = "Predicted strike rate (%)") +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Predicted strike rate by treatment with predicted region composition (Wilson 95% CI).

11 Model performance

11.1 Misclassification analysis

To identify misclassified files, build the first iteration of the model and import its cross-validated predictions (cv_predictions.csv). Subsetting to error_type != "correct" gives the files where the model predictions are wrong.

misclassified <- cv_predictions %>%
  filter(error_type != "correct")

misclassified %>%
  select(file, treatment, passage_type, other_type, probability,
         y_pred, y_true, cv_fold, error_type) %>%
  gt() %>%
  tab_header(title = "Misclassified files (binary model, out-of-fold)") %>%
  fmt_number(columns = probability, decimals = 3) %>%
  gt_style(save_as = "misclassified_files")

Misclassified files (binary model, out-of-fold)
file	treatment	passage_type	other_type	probability	y_pred	y_true	cv_fold	error_type
B304-0609143711	400_model	No contact	NA	0.509	1	0	9	false_positive
B304-0611135809	400_model	Other impeller collision	Impeller surface	0.482	0	1	3	false_negative
B304-0611144937	400_model	Other impeller collision	Impeller surface	0.429	0	1	2	false_negative
B307-0605123359	400_model	Other impeller collision	Impeller surface	0.405	0	1	7	false_negative
B312-0608115526	400_model	Other impeller collision	Impeller hub	0.461	0	1	5	false_negative
B312-0609150543	400_model	Other impeller collision	Impeller hub	0.288	0	1	6	false_negative
B315-0609144811	400_model	Leading edge strike	NA	0.283	0	1	2	false_negative
B321-0605124405	400_model	Other impeller collision	Impeller hub	0.479	0	1	0	false_negative
B321-0609150138	400_model	Other impeller collision	Impeller hub	0.400	0	1	8	false_negative
B61-0102042637	400 (70%)	No contact	NA	0.580	1	0	6	false_positive
B63-0717130205	400 (70%)	No contact	NA	0.535	1	0	1	false_positive
B67-0703174220	500 (100%)	Other impeller collision	Impeller surface	0.446	0	1	3	false_negative
B67-0717145407	400 (70%)	Leading edge strike	NA	0.458	0	1	6	false_negative
B67-0806155637	400 (70%)	No contact	NA	0.523	1	0	7	false_positive
B67-0806173643	400 (70%)	No contact	NA	0.509	1	0	9	false_positive
B76-0716155619	400 (100%)	Other impeller collision	Impeller surface	0.437	0	1	7	false_negative
B76-0717114435	400 (100%)	No contact	NA	0.538	1	0	5	false_positive
B76-0806134449	400 (70%)	Leading edge strike	NA	0.396	0	1	2	false_negative
B77-0103143826	400 (100%)	Other impeller collision	Impeller surface	0.352	0	1	8	false_negative

The workflow for each misclassified file is:

1. Visually check the corresponding video file.
2. If the label is wrong, modify the underlying data label in a copy of the model training dataset, then apply the same modification to the original dataset (data).
3. Retrain the model with the corrected files.

For example, to identify, check and correct a single file:

blade_strike_data_corrected <- blade_strike_data %>%
  mutate(
    passage_type = case_when(
      file == "B312-0608115526" ~ "No contact",
      TRUE ~ as.character(passage_type)
    ),
    leading_type = case_when(
      file == "B312-0608115526" ~ NA_character_,
      TRUE ~ as.character(leading_type)
    ),
    other_type = case_when(
      file == "B312-0608115526" ~ NA_character_,
      TRUE ~ as.character(other_type)
    ),
    concentric_pump_region = case_when(
      file == "B312-0608115526" ~ NA_real_,
      TRUE ~ concentric_pump_region
    )
  )

During this process with the first model iteration, the flagged files were reviewed against video (misclass1_DONE.csv): files with corrected = 1 had their labels corrected, while the remainder kept their original labels as the disagreements were genuine model errors. The corrected data were then joined back to the existing model data to retrain the model.

misclass_review <- misclass1_done %>%
  select(file, treatment, passage_type, other_type, probability,
         error_type, corrected)

misclass_review %>%
  gt() %>%
  tab_header(title = "Reviewed misclassified files and correction outcomes") %>%
  fmt_number(columns = probability, decimals = 3) %>%
  gt_style(save_as = "misclassified_review_outcomes")

Reviewed misclassified files and correction outcomes
file	treatment	passage_type	other_type	probability	error_type	corrected
B300-0609142227	400_model	No contact	NA	0.740	false_positive	1
B300-0609150319	400_model	Other impeller collision	Impeller hub	0.373	false_negative	1
B300-0610105254	400_model	Other impeller collision	Impeller surface	0.329	false_negative	1
B300-0611145724	400_model	Other impeller collision	Impeller hub	0.390	false_negative	1
B302-0611150142	400_model	Other impeller collision	Impeller hub	0.471	false_negative	1
B304-0609145359	400_model	No contact	NA	0.609	false_positive	1
B304-0611135337	400_model	Leading edge strike	NA	0.490	false_negative	1
B304-0611135809	400_model	Other impeller collision	Impeller surface	0.268	false_negative	0
B304-0611144937	400_model	No contact	NA	0.632	false_positive	1
B307-0605123359	400_model	No contact	NA	0.646	false_positive	1
B309-0520151950	400_model	No contact	NA	0.574	false_positive	1
B309-0608113825	400_model	Other impeller collision	Impeller hub	0.457	false_negative	1
B309-0608115848	400_model	No contact	NA	0.547	false_positive	0
B309-0608124015	400_model	Leading edge strike	NA	0.500	false_negative	0
B309-0609144529	400_model	No contact	NA	0.822	false_positive	1
B312-0608115526	400_model	Other impeller collision	Impeller surface	0.356	false_negative	1
B312-0609150543	400_model	Other impeller collision	Impeller hub	0.378	false_negative	0
B312-0611153455	400_model	Leading edge strike	NA	0.438	false_negative	0
B315-0519082900	400_model	Other impeller collision	Impeller hub	0.378	false_negative	1
B315-0520124203	400_model	No contact	NA	0.543	false_positive	1
B315-0609144811	400_model	Leading edge strike	NA	0.262	false_negative	0
B321-0609153202	400_model	Other impeller collision	Impeller hub	0.481	false_negative	0

Of the 22 reviewed files, 15 were corrected and 7 retained their original labels.

11.2 Binary model performance

The binary blade strike model (MiniRocket + RidgeClassifierCV (Binary)) was trained on 565 files (264 strikes, 301 no-strikes) using 10 sensor channels, with 5-fold cross-validation.

cv <- bin_metrics$cross_validation

tibble(
  Metric = c("AUC", "Accuracy", "Sensitivity", "Specificity", "Precision", "F1"),
  Mean = c(cv$mean_auc, cv$mean_accuracy, cv$mean_sensitivity,
           cv$mean_specificity, cv$mean_precision, cv$mean_f1),
  SD = c(cv$std_auc, cv$std_accuracy, cv$std_sensitivity,
         cv$std_specificity, cv$std_precision, cv$std_f1)
) %>%
  gt() %>%
  tab_header(title = "Cross-validation performance (5-fold)") %>%
  fmt_number(columns = c(Mean, SD), decimals = 3) %>%
  gt_style(save_as = "binary_cv_performance")

Cross-validation performance (5-fold)
Metric	Mean	SD
AUC	0.992	0.009
Accuracy	0.975	0.016
Sensitivity	0.977	0.025
Specificity	0.973	0.029
Precision	0.971	0.031
F1	0.974	0.017

oof <- bin_metrics$out_of_fold_performance

oof[map_lgl(oof, is.numeric)] %>%
  unlist() %>%
  enframe("Metric", "Value") %>%
  gt() %>%
  tab_header(title = "Out-of-fold performance") %>%
  fmt_number(columns = Value, decimals = 3) %>%
  gt_style(save_as = "binary_oof_performance")

Out-of-fold performance
Metric	Value
roc_auc	0.991
pr_auc	0.992
overall_accuracy	0.966
sensitivity	0.951
specificity	0.980
precision	0.977
f1_score	0.964
mcc	0.933
FNR	0.049
FPR	0.020
optimal_threshold	0.488

with(oof$confusion_matrix, tibble(
  ` ` = c("Actual: no strike", "Actual: strike"),
  `Predicted: no strike` = c(tn, fn),
  `Predicted: strike` = c(fp, tp)
)) %>%
  gt() %>%
  tab_header(title = "Out-of-fold confusion matrix") %>%
  gt_style(save_as = "binary_confusion_matrix")

Out-of-fold confusion matrix
	Predicted: no strike	Predicted: strike
Actual: no strike	295	6
Actual: strike	13	251

bin_metrics$performance_by_strike_type %>%
  map_dfr(as_tibble, .id = "Strike type") %>%
  gt() %>%
  tab_header(title = "Accuracy by strike type") %>%
  fmt_number(columns = accuracy, decimals = 3) %>%
  gt_style(save_as = "binary_accuracy_by_strike_type")

Accuracy by strike type
Strike type	n_files	accuracy
leading_direct	168	0.982
leading_indirect	32	1.000
no_contact	301	0.980
other_impeller_hub	47	0.915
other_impeller_surface	17	0.647

Overall out-of-fold accuracy is 96.6% with an ROC AUC of 0.991. Performance is high across strike types, with the weakest class being impeller surface collisions. Is this consistent with these being lower-energy contacts that are harder to separate from clean passages?

knitr::include_graphics(c(
  "Model1_1_binary/confusion_matrix.png",
  "Model1_1_binary/roc_curve.png",
  "Model1_1_binary/precision_recall_curve.png",
  "Model1_1_binary/probability_distribution.png"
))

Binary model diagnostics: confusion matrix, ROC curve, precision-recall curve and predicted probability distribution.

knitr::include_graphics("Model1_1_binary/accuracy_by_strike_type.png")

Binary model accuracy by strike type.

11.3 Multiclass model performance

Section to be completed for the multiclass (collapsed concentric region) model, using the outputs in Model1_1_multiclass/.