TBBI vs AMBI-TB Concordance

Setup

Load packages, compute scores for both indexes, and join at the station-year level. anlz_ambiscr returns TBAMBI, the adjusted AMBI-TB score on a 0–10 scale where higher values indicate better benthic conditions – matching the direction of TBBI.

Code

devtools::load_all('C:/proj/tbeptools', quiet = TRUE)
library(dplyr)
library(tidyr)
library(ggplot2)
library(patchwork)
library(knitr)

segs <- c('OTB', 'HB', 'MTB', 'LTB', 'BCB', 'TCB', 'MR')

tbbiscr <- anlz_tbbiscr(benthicdata)
ambiscr  <- anlz_ambiscr(benthicdata, type = 'AMBI-TB')

# Station-year joined dataset; exclude azoic stations and empty TBBI samples
dat_stn <- inner_join(
  tbbiscr |>
    filter(!is.na(TBBI), TBBICat != 'empty sample', AreaAbbr %in% segs) |>
    select(yr, StationNumber, AreaAbbr, TBBI, TBBICat) |> 
    mutate(TBBI = pmax(0, TBBI)),
  ambiscr |>
    filter(!is.na(BC), BC != 7, AreaAbbr %in% segs) |>
    select(yr, StationNumber, TBAMBI, BC),
  by = c('yr', 'StationNumber')
) |>
  filter(!is.na(TBAMBI)) |>
  mutate(AreaAbbr = factor(AreaAbbr, levels = segs))

# Segment-year means and standard errors
seg_yr <- dat_stn |>
  group_by(AreaAbbr, yr) |>
  summarise(
    n         = n(),
    TBBI_mean = mean(TBBI),
    TBBI_se   = sd(TBBI) / sqrt(n()),
    AMBI_mean = mean(TBAMBI),
    AMBI_se   = sd(TBAMBI) / sqrt(n()),
    .groups   = 'drop'
  )

cat('Matched station-years:', nrow(dat_stn), '\n')

Matched station-years: 2983

Code

cat('Segment-years:', nrow(seg_yr), '\n')

Segment-years: 222

Time series by bay segment

Mean ± 1 SE of each index by segment and year. Both axes are scaled so that higher values indicate better benthic conditions (TBBI: 0–100; adjusted AMBI-TB: 0–10).

Code

p_tbbi <- ggplot(seg_yr, aes(x = yr, y = TBBI_mean)) +
  geom_ribbon(aes(ymin = TBBI_mean - TBBI_se, ymax = TBBI_mean + TBBI_se),
              alpha = 0.25, fill = 'steelblue') +
  geom_line(colour = 'steelblue', linewidth = 0.7) +
  facet_wrap(~ AreaAbbr, nrow = 1) +
  labs(x = NULL, y = 'Mean TBBI') +
  theme_bw(base_size = 10) +
  theme(strip.text = element_text(size = 9))

p_ambi <- ggplot(seg_yr, aes(x = yr, y = AMBI_mean)) +
  geom_ribbon(aes(ymin = AMBI_mean - AMBI_se, ymax = AMBI_mean + AMBI_se),
              alpha = 0.25, fill = 'coral') +
  geom_line(colour = 'coral3', linewidth = 0.7) +
  facet_wrap(~ AreaAbbr, nrow = 1) +
  labs(x = 'Year', y = 'Mean adj. AMBI-TB') +
  theme_bw(base_size = 10) +
  theme(strip.text = element_text(size = 9))

p_tbbi / p_ambi

Mean TBBI (top) and adjusted AMBI-TB (bottom) by bay segment and year. Ribbons show ±1 SE.

Scatterplots

Station-year level

Each point is one matched station-year. The blue line is a linear fit with a 95% CI. Pearson r is computed from all non-missing pairs within each segment.

Code

r_labels <- dat_stn |>
  group_by(AreaAbbr) |>
  summarise(
    r   = cor(TBAMBI, TBBI),
    lab = paste0('r = ', round(r, 2)),
    .groups = 'drop'
  )

ggplot(dat_stn, aes(x = TBAMBI, y = TBBI)) +
  geom_point(alpha = 0.2, size = 0.8, colour = 'grey40') +
  geom_smooth(method = 'lm', se = TRUE, colour = 'steelblue', linewidth = 0.7) +
  geom_text(data = r_labels, aes(label = lab),
            x = Inf, y = -Inf, hjust = 1.1, vjust = -0.5, size = 3) +
  facet_wrap(~ AreaAbbr, nrow = 2) +
  labs(x = 'Adjusted AMBI-TB (0-10, higher = better)', y = 'TBBI') +
  theme_bw(base_size = 10)

Adjusted AMBI-TB vs TBBI at the station-year level, by bay segment.

Segment-year level

Each point is one segment-year mean value. Segment-year means reduce within-year station noise and highlight longer-term co-movement between the two indexes.

Code

rho <- cor(seg_yr$AMBI_mean, seg_yr$TBBI_mean, method = 'spearman')

ggplot(seg_yr, aes(x = AMBI_mean, y = TBBI_mean, colour = AreaAbbr)) +
  geom_point(size = 2.2, alpha = 0.85) +
  geom_smooth(aes(group = 1), method = 'lm', se = TRUE,
              colour = 'grey30', linewidth = 0.7) +
  annotate('text', x = Inf, y = -Inf,
           label = paste0('Spearman ρ = ', round(rho, 2)),
           hjust = 1.1, vjust = -0.5, size = 4) +
  scale_colour_brewer(palette = 'Dark2', name = 'Segment') +
  labs(x = 'Mean adj. AMBI-TB', y = 'Mean TBBI') +
  theme_bw(base_size = 11)

Segment-year mean adjusted AMBI-TB vs mean TBBI, coloured by bay segment.

Concordance metrics by segment

Pearson r and Spearman rho computed at the station-year level for each bay segment. Both statistics use all matched, non-missing station-year pairs within each segment.

Code

cor_tab <- dat_stn |>
  group_by(Segment = AreaAbbr) |>
  summarise(
    n              = n(),
    Pearson_r      = round(cor(TBAMBI, TBBI), 3),
    Pearson_p      = cor.test(TBAMBI, TBBI)$p.value,
    Spearman_rho   = round(cor(TBAMBI, TBBI, method = 'spearman'), 3),
    Spearman_p     = cor.test(TBAMBI, TBBI, method = 'spearman')$p.value,
    .groups        = 'drop'
  ) |>
  mutate(across(c(Pearson_p, Spearman_p),
                ~ ifelse(.x < 0.001, '<0.001', as.character(round(.x, 3)))))

kable(cor_tab,
      col.names = c('Segment', 'n', 'Pearson r', 'Pearson p', 'Spearman rho', 'Spearman p'),
      caption   = 'Station-year correlations between adjusted AMBI-TB and TBBI by bay segment.')

Station-year correlations between adjusted AMBI-TB and TBBI by bay segment.
Segment	n	Pearson r	Pearson p	Spearman rho	Spearman p
OTB	570	0.485	<0.001	0.557	<0.001
HB	569	0.285	<0.001	0.487	<0.001
MTB	622	0.372	<0.001	0.420	<0.001
LTB	377	0.317	<0.001	0.099	0.055
BCB	434	0.511	<0.001	0.523	<0.001
TCB	168	0.441	<0.001	0.525	<0.001
MR	243	0.135	0.035	0.350	<0.001

Directional agreement

Does AMBI-TB move in the same direction as TBBI from year to year within each segment? Year-over-year changes in segment-mean values are computed for both indexes. Points landing in the upper-right (both improve) or lower-left (both decline) quadrant indicate agreement.

Code

delta <- seg_yr |>
  arrange(AreaAbbr, yr) |>
  group_by(AreaAbbr) |>
  mutate(
    dTBBI = TBBI_mean - lag(TBBI_mean),
    dAMBI = AMBI_mean - lag(AMBI_mean)
  ) |>
  filter(!is.na(dTBBI), !is.na(dAMBI)) |>
  ungroup()

Code

ggplot(delta, aes(x = dAMBI, y = dTBBI, colour = AreaAbbr)) +
  geom_vline(xintercept = 0, linetype = 'dashed', colour = 'grey60') +
  geom_hline(yintercept = 0, linetype = 'dashed', colour = 'grey60') +
  geom_point(size = 2.5, alpha = 0.8) +
  scale_colour_brewer(palette = 'Dark2', name = 'Segment') +
  labs(x = 'Δ adj. AMBI-TB', y = 'Δ TBBI') +
  theme_bw(base_size = 11)

Year-over-year change in segment-mean AMBI-TB vs TBBI. Dashed lines are zero-change references. Points in quadrants I and III indicate directional agreement.

Code

agree_tab <- delta |>
  group_by(Segment = AreaAbbr) |>
  summarise(
    Transitions       = n(),
    Agree             = sum(sign(dTBBI) == sign(dAMBI)),
    `% Concordance`   = round(100 * Agree / Transitions, 1),
    `Kendall's tau`   = round(cor(dAMBI, dTBBI, method = 'kendall'), 3),
    .groups           = 'drop'
  )

kable(agree_tab,
      caption = 'Directional agreement between adjusted AMBI-TB and TBBI for year-over-year changes in segment means.')

Directional agreement between adjusted AMBI-TB and TBBI for year-over-year changes in segment means.
Segment	Transitions	Agree	% Concordance	Kendall’s tau
OTB	31	22	71.0	0.359
HB	31	16	51.6	0.131
MTB	31	18	58.1	0.260
LTB	31	20	64.5	0.282
BCB	29	20	69.0	0.458
TCB	31	22	71.0	0.415
MR	31	20	64.5	0.239