Expert Elicitation (EE) responses – Tidy + Summaries + Signed Magnitudes + Resampling

---

High-level results

Across both aggregation assumptions, the expert panel indicates:

• Strongly negative counterfactual baseline: the mean expected change by 2042 under business-as-usual is around −12%, with a high proportion of units negative.
• All policy scenarios provide positive uplifts relative to baseline: mean Q2 deltas are positive for every scenario, with a stable ranking of strength (policy_3 strongest; policy_1 weakest).
• But policy is not sufficient on average to deliver net recovery: when combining baseline + uplift, mean absolute outcomes remain negative under all scenarios. policy_3 is closest to stability (near 0) but still negative in the mean, implying widespread improvement but persistent pockets of decline.
• Habitat structure matters: Woodland is closest to neutral outcomes; Freshwater/Wetland and Open habitats remain the most constrained even under the strongest scenario.
• Uncertainty is structured and greatest under the strongest intervention: dispersion and resampling indicate scenario ranking is fairly stable, but the precise “distance to zero” under policy_3 should be treated with uncertainty.

---

1. Background / purpose

This notebook cleans and analyses the Round-2 Expert Elicitation (EE) export to produce policy-relevant diagnostics of:

• the counterfactual baseline expected by experts (Q1),
• the incremental uplift under each policy scenario relative to that baseline (Q2),
• the implied absolute outcome under policy (Q1 + Q2), i.e. whether scenarios are expected to slow decline, stabilise (≈0), or produce net recovery (>0) by 2042,
• and the robustness of conclusions given a small expert panel (dispersion and resampling stability).

The survey is structured around taxa × habitat group combinations. Some taxa are represented by multiple survey_taxa rows within the same habitat group (notably freshwater assemblages and lichens). To make sure headline results are not an artefact of that survey-row multiplicity, the analysis reports key summaries under two transparent aggregation assumptions:

• Assumption A (rows): all survey_taxa rows contribute as-is (faithful to survey structure).
• Assumption B (collapsed): within each expert × scenario, multiple rows that represent the same taxa × habitat_group are averaged so that each unique taxa:habitat_group contributes once.

Where A and B agree, conclusions are robust to the aggregation choice; where they differ, results should be interpreted as sensitive to how survey rows represent ecological “units”.

---

Core framing: baseline vs delta vs absolute outcome

This analysis is structured as expert × taxa:habitat (survey row) × scenario observations.

Policy magnitudes (Q2) are explicitly defined as changes relative to the counterfactual (Q1). Therefore the implied signed outcome under policy is:

\[ \text{Absolute policy outcome (signed \%)} = \text{Counterfactual signed \%} + \text{Policy delta signed \%} \]

This ensures we keep separate:

• baseline decline expectations (Q1),
• incremental policy effects (Q2),
• implied absolute outcomes under policy (Q1 + Q2).

Example (signed magnitudes): - If Q1b = 20 and Q1a indicates “decrease” → counterfactual = −20% - If Q2b = 10 and Q2a indicates “increase” → delta = +10% - Absolute outcome under policy = −20 + 10 = −10%

In the code, this quantity is stored as:

abs_policy_mag_signed = cf_mag_signed + pol_mag_signed

Outputs include summary tables, compact plots, category distributions, and a panel-size sensitivity test (resampling stability).

---

2. Data

# Packages
# ============================================================
library(readxl)
library(dplyr)
library(tidyr)
library(stringr)
library(purrr)
library(readr)
library(ggplot2)
library(forcats)
library(knitr)
library(kableExtra)
library(scales)

# ------------------------------------------------------------
# Input file - update 2 March 2026
# ------------------------------------------------------------
#xlsx_path <- "./data/p591243814707_benjamin.firth_38179087.xlsx"
xlsx_path <- "./data/p591243814707_benjamin.firth_38179087_update_260302.xlsx"


df_raw <- read_excel(xlsx_path)

tibble(
  n_rows = nrow(df_raw),
  n_cols = ncol(df_raw),
  n_Q_cols = sum(str_detect(names(df_raw), "^Q\\d"))
) %>%
  kable(caption = "Raw Excel shape") %>%
  kable_styling(full_width = FALSE)

Raw Excel shape

n_rows	n_cols	n_Q_cols
7	445	420

2.2 Tidy and standardise

2.2.1 Strip HTML formatting

Some respondent cells contain HTML tags (e.g. <span style="...">No meaningful change</span>). We remove tags while preserving the displayed text.

Some, but not all respondent cells contain HTML (e.g. No meaningful change). Remove tags while preserving text:

strip_html <- function(x) {
  if (is.character(x) || inherits(x, "factor")) {
    x2 <- as.character(x)
    x2 <- str_replace_all(x2, "<[^>]+>", "")
    str_squish(x2)
  } else {
    x
  }
}

df_clean <- df_raw %>%
  mutate(across(everything(), strip_html))

2.2.2. New expert ID (id_gj) mapped from Contact_name (keep all experts)

Two respondents are missing some metadata but do have Contact_name. To keep all experts while removing identifiable information, we assign a stable integer ID (1..N) based on Contact_name and do not carry forward names.

df_clean <- df_clean %>%
  mutate(Contact_name = str_squish(as.character(Contact_name)))

contact_key <- df_clean %>%
  distinct(Contact_name) %>%
  filter(!is.na(Contact_name) & Contact_name != "") %>%
  arrange(Contact_name) %>%
  mutate(id_gj = row_number())

# save if needed later
#write.csv(contact_key, "./data/respondent_key_lookup.csv", row.names = FALSE)

# Expert mapping: Contact_name → id_gj
df1 <- df_clean %>%
  left_join(contact_key, by = "Contact_name")

2.2.3 Reshape Q1 + Q2 into df_tidy (expert × survey_taxa × scenario)

We pivot all Q* columns into long format, parse scenario, then pivot wide.

Parsing rules:

• Q1* = "counterfactual" scenario
  
• Q2*_*_{1..4} = "policy_1".."policy_4" scenario

df_q <- df1 %>%
  select(id_gj, matches("^Q\\d"))

long_q <- df_q %>%
  pivot_longer(
    cols = matches("^Q\\d"),
    names_to = "colname",
    values_to = "value",
    values_transform = list(value = as.character) # avoid type conflicts
  ) %>%
  mutate(
    q_code = str_match(colname, "^(Q\\d+[a-z])_\\d+")[,2],
    idx    = as.integer(str_match(colname, "^(Q\\d+[a-z])_(\\d+)")[,3]),
    scen_num = str_match(colname, "^(Q\\d+[a-z])_\\d+_(\\d+)")[,3],
    scen_num = if_else(is.na(scen_num), NA_integer_, as.integer(scen_num)),
    scenario = case_when(
      str_detect(q_code, "^Q1") ~ "counterfactual",
      str_detect(q_code, "^Q2") & !is.na(scen_num) ~ paste0("policy_", scen_num),
      TRUE ~ NA_character_
    ),
    question = str_trim(str_match(colname, ":\\s*([^\\(]+)\\s*\\(")[,2]),
    survey_taxa     = str_trim(str_match(colname, "\\((.*)\\)\\s*$")[,2])
  ) %>%
  mutate(
    col_out = paste0(
      q_code, "_",
      str_replace_all(str_to_lower(question), "[^a-z0-9]+", "_") %>%
        str_replace_all("^_|_$", "")
    )
  ) %>%
  filter(!is.na(scenario)) %>%
  select(id_gj, survey_taxa, scenario, col_out, value)


# reshape_wide
df_tidy <- long_q %>%
  pivot_wider(
    id_cols = c(id_gj, survey_taxa, scenario),
    names_from = col_out,
    values_from = value,
    values_fn = ~ {
      if (length(.x) == 0) NA_character_
      else if (length(.x) == 1) .x
      else paste(.x, collapse = " | ")
    }
  ) %>%
  arrange(id_gj, survey_taxa, scenario)

df_tidy %>%
  count(scenario) %>%
  kable(caption = "Scenario counts in df_tidy") %>%
  kable_styling(full_width = FALSE)

Scenario counts in df_tidy

scenario	n
counterfactual	147
policy_1	147
policy_2	147
policy_3	147
policy_4	147

2.2.4 habitat parsing

IMPORTANT: parsing rules

Goal: 1. Keep habitat_group logic exactly as-is (broad keyword mapping).
2. Set habitat_type to NA for all rows EXCEPT the five listed survey_taxa values, where we want a human-readable, standardised habitat_type label.

Why: - Many survey_taxa strings do not have a trailing “(…)” substring, and for reporting we only want habitat_type populated for the handful of taxa where we explicitly want a more granular label.

Output targets (exact): - “Ditch & Wetland Channel Freshwater invert assemblages” -> “Ditch & Wetland Channel” - “Running-Water Freshwater invert assemblages (Rivers/Streams)” -> “Running-Water (Rivers/Streams)” - “Standing-Water Freshwater invert assemblages (Ponds/Lakes)” -> “Standing-Water (Ponds/Lakes)” - “Saxicolous Lichens (Rock & Stone)” -> “Saxicolous (Rock & Stone)” - “Terricolous Heath/Bog Lichens” -> “Terricolous (Heath/Bog)”

# ============================================================
# Habitat parsing 
# ============================================================

# Helper: 
# Return NA unless survey_taxa is one of the five special cases
habitat_type <- function(survey_taxa) {
  x <- stringr::str_squish(as.character(survey_taxa))

  dplyr::case_when(
    # Freshwater assemblages
    stringr::str_detect(x, "^Ditch\\s*&\\s*Wetland\\s*Channel\\s+Freshwater\\s+invert\\s+assemblages\\s*$") ~
      "Ditch & Wetland Channel",

    stringr::str_detect(x, "^Running\\-Water\\s+Freshwater\\s+invert\\s+assemblages\\s*\\(Rivers/Streams\\)\\s*$") ~
      "Running-Water (Rivers/Streams)",

    stringr::str_detect(x, "^Standing\\-Water\\s+Freshwater\\s+invert\\s+assemblages\\s*\\(Ponds/Lakes\\)\\s*$") ~
      "Standing-Water (Ponds/Lakes)",

    # Lichens
    stringr::str_detect(x, "^Saxicolous\\s+Lichens\\s*\\(Rock\\s*&\\s*Stone\\)\\s*$") ~
      "Saxicolous (Rock & Stone)",

    stringr::str_detect(x, "^Terricolous\\s+Heath/Bog\\s+Lichens\\s*$") ~
      "Terricolous (Heath/Bog)",

    TRUE ~ NA_character_
  )
}

df_tidy <- df_tidy %>%
  dplyr::mutate(
    # --- habitat_type: NA except for five explicitly defined taxa ---
    habitat_type = habitat_type(survey_taxa),

    # Primary classification text for habitat_group mapping:
    # keep our original intent: use habitat_type if present, else survey_taxa.
    hg_txt = stringr::str_to_lower(dplyr::coalesce(habitat_type, survey_taxa)),

    # --- habitat_group mapping (UNCHANGED) ---
    habitat_group = dplyr::case_when(
      stringr::str_detect(hg_txt, "wetland|aquatic|water\\b|water\\-|water\\s+margin|river|stream|rivers|streams|pond|ponds|lake|lakes|ditch|channel|riparian|rheophytic|rheo") ~
        "Freshwater/Wetland",
      stringr::str_detect(hg_txt, "woodland|parkland|scrub|deadwood|epiphyt|humidity|veteran|ancient|saproxylic") ~
        "Woodland",
      stringr::str_detect(hg_txt, "grassland|heath|bog|open\\-habitat|open\\s+habitat|open\\s+ground|bare\\-ground|bare\\s+ground|exposed\\s+substrates?|terricolous") ~
        "Open habitats",
      stringr::str_detect(hg_txt, "rock|stone|saxicolous") ~
        "Open habitats",
      TRUE ~ "Unknown"
    ),

    # Force non-missing even if something slips through
    habitat_group = dplyr::if_else(is.na(habitat_group) | habitat_group == "", "Unknown", habitat_group),

    # taxa grouping (unchanged idea; keyword-based)
    taxa = dplyr::case_when(
      stringr::str_detect(survey_taxa, stringr::regex("\\bDiptera\\b", ignore_case = TRUE)) ~ "Diptera",
      stringr::str_detect(survey_taxa, stringr::regex("\\bBeetles\\b", ignore_case = TRUE)) ~ "Beetles",
      stringr::str_detect(survey_taxa, stringr::regex("Vascular Plants\\b", ignore_case = TRUE)) ~ "Vascular plants",
      stringr::str_detect(survey_taxa, stringr::regex("\\bLichens\\b", ignore_case = TRUE)) ~ "Lichens",
      stringr::str_detect(survey_taxa, stringr::regex("\\bBryophytes?\\b", ignore_case = TRUE)) ~ "Bryophytes",
      stringr::str_detect(survey_taxa, stringr::regex("\\bSpiders\\b", ignore_case = TRUE)) ~ "Spiders",
      stringr::str_detect(survey_taxa, stringr::regex("Freshwater invert assemblages\\b", ignore_case = TRUE)) ~ "Freshwater invertebrate assemblages",
      TRUE ~ "Other/unknown"
    )
  ) %>%
  dplyr::select(-hg_txt)

# Diagnostics: do we have any NA habitat_group now?
df_tidy %>%
  summarise(
    n_rows = n(),
    n_habitat_group_na = sum(is.na(habitat_group)),
    n_habitat_group_unknown = sum(habitat_group == "Unknown", na.rm = TRUE)
  ) %>%
  kable(caption = "habitat_group diagnostics (should have 0 NA)") %>%
  kable_styling(full_width = FALSE)

habitat_group diagnostics (should have 0 NA)

n_rows	n_habitat_group_na	n_habitat_group_unknown
735	0	0

# Coverage: show habitat_group by survey_taxa (useful spot-check)
df_tidy %>%
  distinct(survey_taxa, habitat_type, habitat_group) %>%
  arrange(habitat_group, survey_taxa) %>%
  kable(caption = "Mapping check: survey_taxa to habitat_group") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "320px")

Mapping check: survey_taxa to habitat_group

survey_taxa	habitat_type	habitat_group
Ditch & Wetland Channel Freshwater invert assemblages	Ditch & Wetland Channel	Freshwater/Wetland
Running-Water Freshwater invert assemblages (Rivers/Streams)	Running-Water (Rivers/Streams)	Freshwater/Wetland
Standing-Water Freshwater invert assemblages (Ponds/Lakes)	Standing-Water (Ponds/Lakes)	Freshwater/Wetland
Wetland & Aquatic Beetles	NA	Freshwater/Wetland
Wetland & Aquatic Vascular Plants	NA	Freshwater/Wetland
Wetland & Rheophytic Bryophytes	NA	Freshwater/Wetland
Wetland & Riparian Spiders	NA	Freshwater/Wetland
Wetland & Water-Margin Diptera	NA	Freshwater/Wetland
Open Grassland & Heath Vascular Plants	NA	Open habitats
Open Grassland/Heath/Bare-Ground Spiders	NA	Open habitats
Open-Ground Bryophytes (Grassland, Heath, Exposed Substrates)	NA	Open habitats
Open-Habitat & Flower-Visiting Diptera	NA	Open habitats
Open-Habitat Ground Beetles (Grassland & Heath)	NA	Open habitats
Saxicolous Lichens (Rock & Stone)	Saxicolous (Rock & Stone)	Open habitats
Terricolous Heath/Bog Lichens	Terricolous (Heath/Bog)	Open habitats
Epiphytic Lichens (Woodland & Parkland)	NA	Woodland
Saproxylic & Woodland Diptera	NA	Woodland
Saproxylic Beetles (Woodland Deadwood)	NA	Woodland
Woodland & Scrub Spiders	NA	Woodland
Woodland & Scrub Vascular Plants	NA	Woodland
Woodland Humidity & Epiphyte Bryophytes	NA	Woodland

2.3 Scoring categories and signed magnitudes

We now translate/compute:

• Distribution score (1–5)
• Certainty score (1–5)
• Signed counterfactual magnitude (%)
• Signed policy delta magnitude (%
• Implied signed policy outcome (%): cf_mag_signed + pol_mag_signed

Important convention: responses of “do not know/unsure” are treated the same as missing.

# ============================================================
# Scoring helpers
# ============================================================
score_dist <- function(x){
  x0 <- str_squish(tolower(as.character(x)))
  x0[x0 %in% c("do not know/unsure","do not know","unsure","don't know","dont know","dk","")] <- NA_character_

  case_when(
    is.na(x0) ~ NA_real_,
    str_detect(x0, "robust|significant|large") & str_detect(x0, "decrease") ~ 1,
    str_detect(x0, "moderate") & str_detect(x0, "decrease") ~ 2,
    str_detect(x0, "no meaningful|no change") ~ 3,
    str_detect(x0, "moderate") & str_detect(x0, "increase") ~ 4,
    str_detect(x0, "robust|significant|large") & str_detect(x0, "increase") ~ 5,
    TRUE ~ NA_real_
  )
}

score_cert <- function(x){
  x0 <- str_squish(tolower(as.character(x)))
  x0[x0 %in% c("do not know/unsure","do not know","unsure","don't know","dont know","dk","")] <- NA_character_

  case_when(
    is.na(x0) ~ NA_real_,
    x0 %in% c("very low","very_low","very-low") ~ 1,
    x0 %in% c("low") ~ 2,
    x0 %in% c("medium","moderate") ~ 3,
    x0 %in% c("high") ~ 4,
    x0 %in% c("very high","very_high","very-high") ~ 5,
    TRUE ~ NA_real_
  )
}

dir_sign <- function(x){
  x0 <- str_squish(tolower(as.character(x)))
  x0[x0 %in% c("do not know/unsure","do not know","unsure","don't know","dont know","dk","")] <- NA_character_

  case_when(
    is.na(x0) ~ NA_real_,
    str_detect(x0, "decrease") ~ -1,
    str_detect(x0, "increase") ~  1,
    str_detect(x0, "no meaningful|no change") ~ 0,
    TRUE ~ NA_real_
  )
}


# ============================================================
# Score and sign responses
# ============================================================
Q1a <- "Q1a_likely_change_in_distribution"
Q1b <- "Q1b_magnitude_of_change"
Q1c <- "Q1c_your_certainty"

Q2a <- "Q2a_likely_change_in_distribution_compared_to_counterfactual"
Q2b <- "Q2b_magnitude_of_change_compared_to_the_counterfactual"
Q2c <- "Q2c_your_certainty"

df_scored <- df_tidy %>%
  mutate(
    # parse numeric magnitudes (percent scale assumed in survey)
    Q1b_num = parse_number(as.character(.data[[Q1b]]), na = c("N/A","NA","na","n/a","")),
    Q2b_num = parse_number(as.character(.data[[Q2b]]), na = c("N/A","NA","na","n/a","")),

    # scored categories
    cf_dist = score_dist(.data[[Q1a]]),
    pol_dist = score_dist(.data[[Q2a]]),

    cf_cert = score_cert(.data[[Q1c]]),
    pol_cert = score_cert(.data[[Q2c]]),

    # direction
    cf_dir = dir_sign(.data[[Q1a]]),
    pol_dir = dir_sign(.data[[Q2a]]),

    # signed magnitudes
    cf_mag_signed = if_else(is.na(Q1b_num) | is.na(cf_dir), NA_real_, Q1b_num * cf_dir),
    pol_mag_signed = if_else(is.na(Q2b_num) | is.na(pol_dir), NA_real_, Q2b_num * pol_dir)
  )

# ============================================================
# Attach counterfactual (Q1) values onto policy rows
# so we can compute implied absolute outcomes later
# ============================================================
scenario_levels <- c("policy_1","policy_2","policy_3","policy_4")

cf_lookup <- df_scored %>%
  filter(scenario == "counterfactual") %>%
  transmute(
    id_gj, survey_taxa,
    cf_dist2 = cf_dist,
    cf_cert2 = cf_cert,
    cf_mag_signed2 = cf_mag_signed
  )

df_scored2 <- df_scored %>%
  left_join(cf_lookup, by = c("id_gj", "survey_taxa")) %>%
  mutate(
    cf_dist = cf_dist2,
    cf_cert = cf_cert2,
    cf_mag_signed = cf_mag_signed2
  ) %>%
  select(-cf_dist2, -cf_cert2, -cf_mag_signed2) %>%
  mutate(
    abs_policy_mag_signed = if_else(
      scenario %in% scenario_levels,
      cf_mag_signed + pol_mag_signed,
      cf_mag_signed
    ),
    delta_cert = if_else(
      scenario %in% scenario_levels,
      pol_cert - cf_cert,
      NA_real_
    )
  )

df_scored2 %>%
  summarise(
    n_rows = n(),
    n_experts = n_distinct(id_gj),
    n_survey_taxa = n_distinct(survey_taxa),
    n_cf_signed = sum(!is.na(cf_mag_signed)),
    n_pol_delta_signed = sum(!is.na(pol_mag_signed)),
    n_abs_policy_signed = sum(!is.na(abs_policy_mag_signed))
  ) %>%
  kable(caption = "Coverage: signed magnitudes and absolute policy outcomes") %>%
  kable_styling(full_width = FALSE)

Coverage: signed magnitudes and absolute policy outcomes

n_rows	n_experts	n_survey_taxa	n_cf_signed	n_pol_delta_signed	n_abs_policy_signed
735	7	21	705	564	705

IMPORTANT assumption: we treat the response category “do not know/unsure” the same way as we treat NAs (ie respondent leaving questions empty == respondents selecting “do not know/unsure”)

# Expected columns from reshape
Q1a <- "Q1a_likely_change_in_distribution"
Q1b <- "Q1b_magnitude_of_change"
Q1c <- "Q1c_your_certainty"

Q2a <- "Q2a_likely_change_in_distribution_compared_to_counterfactual"
Q2b <- "Q2b_magnitude_of_change_compared_to_the_counterfactual"
Q2c <- "Q2c_your_certainty"




df_scored <- df_tidy %>%
  mutate(
    Q1b_num = parse_number(as.character(.data[[Q1b]]), na = c("N/A","NA","na","n/a","")),
    Q2b_num = parse_number(as.character(.data[[Q2b]]), na = c("N/A","NA","na","n/a","")),

    cf_dist = score_dist(.data[[Q1a]]),
    pol_dist = score_dist(.data[[Q2a]]),

    cf_cert = score_cert(.data[[Q1c]]),
    pol_cert = score_cert(.data[[Q2c]]),

    cf_dir = dir_sign(.data[[Q1a]]),
    pol_dir = dir_sign(.data[[Q2a]]),

    cf_mag_signed = if_else(is.na(Q1b_num) | is.na(cf_dir), NA_real_, Q1b_num * cf_dir),
    pol_mag_signed = if_else(is.na(Q2b_num) | is.na(pol_dir), NA_real_, Q2b_num * pol_dir),
    abs_policy_mag_signed = cf_mag_signed + pol_mag_signed,

    delta_cert = pol_cert - cf_cert
  )

df_scored %>%
  summarise(
    n_rows = n(),
    n_cf_mag = sum(!is.na(Q1b_num)),
    n_pol_mag = sum(!is.na(Q2b_num)),
    n_cf_signed = sum(!is.na(cf_mag_signed)),
    n_delta_signed = sum(!is.na(pol_mag_signed))
  ) %>%
  kable(caption = "Coverage: magnitudes and signed magnitudes") %>%
  kable_styling(full_width = FALSE)

Coverage: magnitudes and signed magnitudes

n_rows	n_cf_mag	n_pol_mag	n_cf_signed	n_delta_signed
735	141	564	141	564

# ============================================================
# FIX: attach counterfactual (Q1) values onto policy rows
# so abs_policy_mag_signed and delta_cert can be computed
# ============================================================

scenario_levels <- c("policy_1","policy_2","policy_3","policy_4")

# 1) Build a CF lookup (one row per id_gj × survey_taxa)
cf_lookup <- df_scored %>%
  filter(scenario == "counterfactual") %>%
  transmute(
    id_gj, survey_taxa,
    cf_dist2 = cf_dist,
    cf_cert2 = cf_cert,
    cf_mag_signed2 = cf_mag_signed
  )

# 2) Join CF values onto ALL rows (including policy rows)
df_scored2 <- df_scored %>%
  left_join(cf_lookup, by = c("id_gj", "survey_taxa")) %>%
  mutate(
    # overwrite cf_* with joined versions (so they exist on policy rows too)
    cf_dist = cf_dist2,
    cf_cert = cf_cert2,
    cf_mag_signed = cf_mag_signed2
  ) %>%
  select(-cf_dist2, -cf_cert2, -cf_mag_signed2) %>%
  mutate(
    # Now these work (because cf_* exist on policy rows)
    abs_policy_mag_signed = if_else(
      scenario %in% scenario_levels,
      cf_mag_signed + pol_mag_signed,
      cf_mag_signed
    ),
    delta_cert = if_else(
      scenario %in% scenario_levels,
      pol_cert - cf_cert,
      NA_real_
    )
  )

# 3) Quick sanity check: abs_policy_mag_signed should no longer be all NA
df_scored2 %>%
  group_by(scenario) %>%
  summarise(
    n = n(),
    n_abs_policy_mag = sum(!is.na(abs_policy_mag_signed)),
    n_delta_cert = sum(!is.na(delta_cert)),
    .groups = "drop"
  ) %>% print()

## # A tibble: 5 × 4
##   scenario           n n_abs_policy_mag n_delta_cert
##   <chr>          <int>            <int>        <int>
## 1 counterfactual   147              141            0
## 2 policy_1         147              141          147
## 3 policy_2         147              141          147
## 4 policy_3         147              141          147
## 5 policy_4         147              141          147

The dataset is balanced across scenarios and experts. All scenarios are evaluated on the same underlying baseline expectations, ensuring comparability.

3. Aggregation assumptions: row-based vs taxa:habitat_group

The elicitation design is such that the data includes multiple survey_taxa rows for some taxa within a single habitat group (notably freshwater assemblages and lichens). If we compute means across all survey_taxa rows, these taxa would contribute more weight to habitat-group and overall summaries simply because they have more rows.

To make this transparent, we run the analysis under two assumptions:

Assumption A (Row-based: all unique survey_taxa contribute)

Unit: > expert × scenario × survey_taxa

This preserves the survey structure exactly.

Assumption B (Collapsed within habitat: unique taxa:habitat_group contribute once)

Unit: > expert × scenario × taxa × habitat_group

This averages across multiple survey_taxa rows that fall into the same taxa:habitat_group cell for each expert and scenario.

Both assumptions are shown side-by-side for: - overall scenario summaries, and - habitat_group summaries.

---

scenario_levels <- c("counterfactual","policy_1","policy_2","policy_3","policy_4")

# A) Row-based: keep all survey_taxa rows (no collapsing)
df_A <- df_scored2 %>%
  filter(scenario %in% scenario_levels) %>%
  filter(!(taxa == "Lichens" & habitat_group == "Freshwater/Wetland")) %>%
  mutate(assumption = "A (rows = all survey_taxa)")

# B) Collapse within habitat: unique taxa:habitat_group per expert × scenario
df_B <- df_A %>%
  group_by(id_gj, scenario, habitat_group, taxa) %>%
  summarise(
    cf_mag_signed = mean(cf_mag_signed, na.rm = TRUE),
    pol_mag_signed = mean(pol_mag_signed, na.rm = TRUE),
    abs_policy_mag_signed = mean(abs_policy_mag_signed, na.rm = TRUE),
    cf_dist = mean(cf_dist, na.rm = TRUE),
    pol_dist = mean(pol_dist, na.rm = TRUE),
    cf_cert = mean(cf_cert, na.rm = TRUE),
    pol_cert = mean(pol_cert, na.rm = TRUE),
    n_survey_taxa_rows_collapsed = n(),
    .groups = "drop"
  ) %>%
  mutate(assumption = "B (collapsed = unique taxa:habitat_group)")

# Combined object (handy for faceting tables/plots)
df_AB <- bind_rows(df_A, df_B)

# Diagnostics: show where collapsing occurred
df_B %>%
  filter(n_survey_taxa_rows_collapsed > 1) %>%
  count(habitat_group, taxa, n_survey_taxa_rows_collapsed, sort = TRUE) %>%
  kable(caption = "Where Assumption B collapses multiple survey_taxa rows into one taxa:habitat_group unit") %>%
  kable_styling(full_width = FALSE)

Where Assumption B collapses multiple survey_taxa rows into one taxa:habitat_group unit

habitat_group	taxa	n_survey_taxa_rows_collapsed	n
Freshwater/Wetland	Freshwater invertebrate assemblages	3	35
Open habitats	Lichens	2	35

---

4. Expected trajectory under the counterfactual (Q1 baseline)

This section establishes the counterfactual (“business-as-usual”) baseline from Q1, before we interpret policy deltas from Q2. Because Q2 magnitudes are explicitly defined relative to Q1, we need a clear view of baseline “pull” to answer the policy-relevant question: do uplifts merely slow decline, or do they offset it sufficiently to deliver stability/recovery by 2042?

We report results under two aggregation assumptions:

• Assumption A (rows = all survey_taxa): each survey row contributes as an observation. This reflects the survey design directly, but can implicitly overweight taxa that were split into multiple survey_taxa rows within a habitat group (notably freshwater assemblages and lichens).

• Assumption B (collapsed = unique taxa:habitat_group): within each expert × scenario, multiple survey_taxa rows that represent the same taxa within the same habitat group are collapsed to one unit (averaged). This is the most direct way to prevent survey-row multiplicity from driving habitat-group summaries.

Importantly, taxa-level comparisons are not expected to change materially across A vs B, because the collapsing only affects where a taxon has multiple survey rows within the same habitat group. By construction, the main place where A vs B can matter is therefore habitat-group summaries.

cf_AB <- df_AB %>%
  filter(scenario == "counterfactual")

# Headline baseline (A vs B)
cf_AB %>%
  group_by(assumption) %>%
  summarise(
    n_units = n(),
    n_experts = n_distinct(id_gj),
    n_taxa = n_distinct(taxa),
    n_hab_groups = n_distinct(habitat_group),
    mean_cf = mean(cf_mag_signed, na.rm = TRUE),
    median_cf = median(cf_mag_signed, na.rm = TRUE),
    sd_cf = sd(cf_mag_signed, na.rm = TRUE),
    p10 = quantile(cf_mag_signed, 0.10, na.rm = TRUE),
    p90 = quantile(cf_mag_signed, 0.90, na.rm = TRUE),
    prop_negative = mean(cf_mag_signed < 0, na.rm = TRUE),
    prop_zero = mean(cf_mag_signed == 0, na.rm = TRUE),
    prop_positive = mean(cf_mag_signed > 0, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  kable(caption = "Counterfactual baseline (Q1 signed %): Assumptions A and B") %>%
  kable_styling(full_width = FALSE)

Counterfactual baseline (Q1 signed %): Assumptions A and B

assumption	n_units	n_experts	n_taxa	n_hab_groups	mean_cf	median_cf	sd_cf	p10	p90	prop_negative	prop_zero	prop_positive
A (rows = all survey_taxa)	147	7	7	3	-12.33	-10	11.46	-30	0	0.79	0.13	0.07
B (collapsed = unique taxa:habitat_group)	126	7	7	3	-11.95	-10	11.74	-30	0	0.78	0.15	0.07

---

4.1 Baseline by habitat_group (Assumptions A and B)

How to read the baseline tables (and why A vs B matters here)

The baseline tables report:

• mean_cf: mean signed counterfactual change (%), where negative values imply decline by 2042.
• prop_negative: share of units with baseline decline (< 0), which is a compact signal of how widespread negative expectations are.
• mean_cert: average reported certainty (1–5 scale).

At headline level, the baseline is strongly negative under both assumptions (mean around −12% in A and slightly less negative in B), and the proportion negative is high (around four-fifths). This indicates that, in expert judgement, continued decline is the default state without further intervention.

The habitat-group baseline table is shown explicitly under both assumptions because this is where multiplicity can matter. In these data, the only systematic collapsing in B is:

• Freshwater/Wetland × Freshwater invertebrate assemblages (3 survey_taxa rows collapsed)
• Open habitats × Lichens (2 survey_taxa rows collapsed)

so A can overweight those components in habitat-group means relative to other taxa that appear as a single survey row within habitat.

This table is shown under both assumptions because it is precisely where survey-row multiplicity can matter.

cf_by_hab_AB <- cf_AB %>%
  group_by(assumption, habitat_group) %>%
  summarise(
    n = n(),
    mean_cf = mean(cf_mag_signed, na.rm = TRUE),
    prop_negative = mean(cf_mag_signed < 0, na.rm = TRUE),
    mean_cert = mean(cf_cert, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(assumption, mean_cf)

cf_by_hab_AB %>%
  mutate(
    mean_cf = round(mean_cf, 2),
    mean_cert = round(mean_cert, 2),
    prop_negative = scales::percent(prop_negative, accuracy = 1)
  ) %>%
  kable(caption = "Counterfactual baseline by habitat_group (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Counterfactual baseline by habitat_group (Assumptions A and B)

assumption	habitat_group	n	mean_cf	prop_negative	mean_cert
A (rows = all survey_taxa)	Freshwater/Wetland	56	-13.83	87%	2.02
A (rows = all survey_taxa)	Open habitats	49	-13.43	85%	1.98
A (rows = all survey_taxa)	Woodland	42	-9.00	62%	1.88
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	42	-13.47	85%	1.94
B (collapsed = unique taxa:habitat_group)	Open habitats	42	-13.36	85%	1.95
B (collapsed = unique taxa:habitat_group)	Woodland	42	-9.00	62%	1.88

two_cols <- c("#c9643e", "#618f72")

ggplot(cf_by_hab_AB, aes(x = habitat_group, y = mean_cf, fill = assumption)) +
  geom_col(position = position_dodge(width = 0.8)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Counterfactual baseline by habitat_group (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % change (counterfactual)"
  ) +
  theme_minimal()

ggplot(subset(cf_by_hab_AB, assumption = "A (rows = all survey_taxa)"), aes(x = habitat_group, y = mean_cf)) +
  geom_col(position = position_dodge(width = 0.8)) +
  coord_flip() +
  geom_hline(yintercept = 0, linetype = "dashed") +
 # scale_fill_manual(values = two_cols) +
  labs(
    title = "Counterfactual baseline by habitat_group (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % change (counterfactual)"
  ) +
  theme_minimal(base_size = 18) 

## Warning: In subset.data.frame(cf_by_hab_AB, assumption = "A (rows = all survey_taxa)") :
##  extra argument 'assumption' will be disregarded

Interpretation: baseline differs by habitat group, but A vs B does not change the story

Across both assumptions, Freshwater/Wetland and Open habitats show the most negative baseline means, while Woodland is less negative. Under Assumption A, the mean counterfactual baseline is approximately:

• Freshwater/Wetland: −13.83% (87% negative)
• Open habitats: −13.43% (85% negative)
• Woodland: −9.00% (62% negative)

Under Assumption B, the corresponding habitat-group means are very similar (Freshwater/Wetland and Open habitats remain strongly negative; Woodland remains less negative), with small shifts consistent with removing survey-row multiplicity.

Policy implication at this stage: the baseline “pull” is large enough that even sizeable positive deltas may still leave net negative outcomes once CF + delta are combined. This is exactly what Sections 5–6 test.

---

4.2 Baseline by taxa (A only; A == B)

Taxa-level comparisons are not meaningfully altered by collapsing within habitat, so we report once using the row-based structure.

cf_by_taxa <- df_A %>%
  filter(scenario == "counterfactual") %>%
  group_by(taxa) %>%
  summarise(
    n = n(),
    mean_cf = mean(cf_mag_signed, na.rm = TRUE),
    prop_negative = mean(cf_mag_signed < 0, na.rm = TRUE),
    mean_cert = mean(cf_cert, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(mean_cf)

cf_by_taxa %>%
  mutate(
    mean_cf = round(mean_cf, 2),
    mean_cert = round(mean_cert, 2),
    prop_negative = scales::percent(prop_negative, accuracy = 1)
  ) %>%
  kable(caption = "Counterfactual baseline by taxa (A; identical interpretation under B)") %>%
  kable_styling(full_width = FALSE)

Counterfactual baseline by taxa (A; identical interpretation under B)

taxa	n	mean_cf	prop_negative	mean_cert
Freshwater invertebrate assemblages	21	-14.86	95%	2.24
Vascular plants	21	-13.33	76%	2.19
Lichens	21	-13.00	90%	2.14
Bryophytes	21	-12.62	86%	2.00
Diptera	21	-12.39	78%	1.67
Beetles	21	-10.48	67%	1.95
Spiders	21	-9.17	61%	1.57

---

4.3 Compact counterfactual plots (A vs B)

ggplot(cf_AB, aes(x = cf_mag_signed, fill = assumption)) +
  geom_histogram(bins = 18, color = "white") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  facet_wrap(~ assumption) +
  scale_fill_manual(values = two_cols) +
  guides(fill = "none") +
  labs(
    title = "Counterfactual (Q1): distribution of signed magnitude (% change)",
    x = "Signed % change under counterfactual",
    y = "Count"
  ) +
  theme_minimal()

---

5. Policy uplifts relative to the counterfactual (Q2 deltas)

Section 5 interprets Q2 policy deltas, which are defined as incremental change relative to the counterfactual. These are not absolute outcomes. A positive delta means “less decline than counterfactual” (or “more increase than counterfactual”), but it does not tell us whether the system is stable/recovering overall.

We therefore treat Section 5 as answering:

1. Are the scenarios directionally beneficial relative to baseline?
   
2. How large are the uplifts, and how consistent are they?
   
3. Do uplifts differ across habitat groups, and does the A vs B aggregation choice affect those habitat-group comparisons?

The coverage table confirms that both A and B are based on the same expert panel (n=7), with B having fewer units because of the within-habitat collapsing step.

Q2 magnitudes are defined relative to the counterfactual. We show overall and habitat-group summaries under both assumptions.

pol_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_"))

pol_AB %>%
  group_by(assumption) %>%
  summarise(
    n_units = n(),
    n_nonmiss = sum(!is.na(pol_mag_signed)),
    n_experts = n_distinct(id_gj),
    .groups = "drop"
  ) %>%
  kable(caption = "Incremental policy delta (Q2) coverage check (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Incremental policy delta (Q2) coverage check (Assumptions A and B)

assumption	n_units	n_nonmiss	n_experts
A (rows = all survey_taxa)	588	564	7
B (collapsed = unique taxa:habitat_group)	504	480	7

5.1 Overall incremental uplifts by scenario (Assumptions A and B)

Interpretation: all policies uplift relative to counterfactual; policy_3 is consistently strongest

The scenario-level delta table shows positive mean uplifts for every policy scenario under both assumptions. Under Assumption A (row-based), the mean deltas are:

• policy_1: +2.66%
• policy_2: +7.77%
• policy_4: +7.34%
• policy_3: +10.60% (largest)

This produces a clear and stable ranking: policy_3 > policy_2 ≈ policy_4 > policy_1.

The “prop_positive” column shows that most expert–unit responses are positive deltas (for example, under policy_3 it is very high), meaning experts generally agree policies improve outcomes relative to baseline even when they disagree about magnitude.

A vs B: At the overall scenario level, collapsing multiplicity (B) should not materially change the ordering and typically shifts headline means only modestly, because only a small subset of taxa:habitat combinations are affected by the collapse step. The main sensitivity to A vs B is expected to show up in habitat-group summaries, not in the overall ranking.

delta_headline_AB <- pol_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_delta = mean(pol_mag_signed, na.rm = TRUE),
    median_delta = median(pol_mag_signed, na.rm = TRUE),
    sd_delta = sd(pol_mag_signed, na.rm = TRUE),
    p10 = quantile(pol_mag_signed, 0.10, na.rm = TRUE),
    p90 = quantile(pol_mag_signed, 0.90, na.rm = TRUE),
    prop_positive = mean(pol_mag_signed > 0, na.rm = TRUE),
    prop_zero = mean(pol_mag_signed == 0, na.rm = TRUE),
    prop_negative = mean(pol_mag_signed < 0, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(assumption, desc(mean_delta))

delta_headline_AB %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  kable(caption = "Incremental policy uplifts (Q2 deltas): overall by scenario (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Incremental policy uplifts (Q2 deltas): overall by scenario (Assumptions A and B)

assumption	scenario	mean_delta	median_delta	sd_delta	p10	p90	prop_positive	prop_zero	prop_negative
A (rows = all survey_taxa)	policy_3	16.55	15	6.43	8.00	25.0	0.99	0.01	0.00
A (rows = all survey_taxa)	policy_4	11.81	13	5.02	5.00	18.0	0.96	0.04	0.00
A (rows = all survey_taxa)	policy_2	10.75	10	5.77	4.00	20.0	0.94	0.06	0.00
A (rows = all survey_taxa)	policy_1	4.15	5	6.98	-10.00	10.0	0.73	0.13	0.14
B (collapsed = unique taxa:habitat_group)	policy_3	16.58	15	6.34	7.90	25.0	0.98	0.02	0.00
B (collapsed = unique taxa:habitat_group)	policy_4	11.81	13	4.98	5.90	17.1	0.96	0.04	0.00
B (collapsed = unique taxa:habitat_group)	policy_2	10.94	10	5.63	4.95	20.0	0.97	0.03	0.00
B (collapsed = unique taxa:habitat_group)	policy_1	4.24	5	6.96	-10.00	10.0	0.77	0.09	0.14

ggplot(delta_headline_AB, aes(x = scenario, y = mean_delta, fill = assumption)) +
  geom_col(position = position_dodge(width = 0.8)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Incremental policy uplifts by scenario (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % delta"
  ) +
  theme_minimal()

5.2 Uplifts by habitat_group (Assumptions A and B)

delta_by_hab_AB <- pol_AB %>%
  group_by(assumption, habitat_group, scenario) %>%
  summarise(
    n = n(),
    mean_delta = mean(pol_mag_signed, na.rm = TRUE),
    prop_positive = mean(pol_mag_signed > 0, na.rm = TRUE),
    .groups = "drop"
  )

delta_by_hab_AB %>%
  mutate(
    mean_delta = round(mean_delta, 2),
    prop_positive = scales::percent(prop_positive, accuracy = 1)
  ) %>%
  arrange(assumption, habitat_group, scenario) %>%
  kable(caption = "Incremental policy uplifts by habitat_group (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "320px")

Incremental policy uplifts by habitat_group (Assumptions A and B)

assumption	habitat_group	scenario	n	mean_delta	prop_positive
A (rows = all survey_taxa)	Freshwater/Wetland	policy_1	56	4.44	74%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_2	56	10.83	93%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_3	56	17.61	100%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_4	56	13.13	98%
A (rows = all survey_taxa)	Open habitats	policy_1	49	3.94	72%
A (rows = all survey_taxa)	Open habitats	policy_2	49	10.81	94%
A (rows = all survey_taxa)	Open habitats	policy_3	49	16.32	98%
A (rows = all survey_taxa)	Open habitats	policy_4	49	11.51	98%
A (rows = all survey_taxa)	Woodland	policy_1	42	4.00	72%
A (rows = all survey_taxa)	Woodland	policy_2	42	10.57	98%
A (rows = all survey_taxa)	Woodland	policy_3	42	15.40	98%
A (rows = all survey_taxa)	Woodland	policy_4	42	10.38	92%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_1	42	4.35	78%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_2	42	10.89	98%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_3	42	17.59	100%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_4	42	13.12	98%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_1	42	4.38	80%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_2	42	11.35	95%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_3	42	16.75	98%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_4	42	11.93	98%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_1	42	4.00	72%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_2	42	10.57	98%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_3	42	15.40	98%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_4	42	10.38	92%

ggplot(delta_by_hab_AB, aes(x = scenario, y = mean_delta, fill = assumption)) +
  geom_col(position = position_dodge(width = 0.8)) +
  facet_wrap(~ habitat_group) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Incremental policy uplifts by habitat_group (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % delta"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 0.5)
  )

Interpretation: uplifts vary by habitat group; A vs B matters most here (by design)

Habitat-group uplift summaries provide the first indication of where policy leverage is expected to be strongest. The key point is not just “uplifts are positive”, but whether uplifts in the most-negative baseline habitat groups (Freshwater/Wetland and Open habitats) are large enough to offset the baseline pull established in Section 4.

Because Assumption B collapses the freshwater assemblage multiplicity (and lichen multiplicity in open habitats), the habitat-group mean delta is the place to look for any material change between A vs B. If the habitat-group conclusions are robust across A and B, that strengthens confidence that the reported habitat differences are not artefacts of how the survey rows were structured.

5.3 Uplifts by taxa (A only; A == B)

pol_A <- df_A %>%
  filter(str_detect(scenario, "^policy_"))

delta_by_taxa <- pol_A %>%
  group_by(taxa, scenario) %>%
  summarise(
    n = n(),
    mean_delta = mean(pol_mag_signed, na.rm = TRUE),
    prop_positive = mean(pol_mag_signed > 0, na.rm = TRUE),
    .groups = "drop"
  )

delta_by_taxa %>%
  mutate(
    mean_delta = round(mean_delta, 2),
    prop_positive = scales::percent(prop_positive, accuracy = 1)
  ) %>%
  kable(caption = "Incremental policy uplifts by taxa (A; identical interpretation under B)") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "320px")

Incremental policy uplifts by taxa (A; identical interpretation under B)

taxa	scenario	n	mean_delta	prop_positive
Beetles	policy_1	21	4.76	71%
Beetles	policy_2	21	12.43	100%
Beetles	policy_3	21	17.10	100%
Beetles	policy_4	21	11.71	95%
Bryophytes	policy_1	21	3.24	76%
Bryophytes	policy_2	21	9.95	95%
Bryophytes	policy_3	21	16.10	95%
Bryophytes	policy_4	21	11.52	95%
Diptera	policy_1	21	3.89	67%
Diptera	policy_2	21	10.44	89%
Diptera	policy_3	21	15.83	94%
Diptera	policy_4	21	11.56	83%
Freshwater invertebrate assemblages	policy_1	21	4.71	71%
Freshwater invertebrate assemblages	policy_2	21	10.67	86%
Freshwater invertebrate assemblages	policy_3	21	17.67	100%
Freshwater invertebrate assemblages	policy_4	21	13.14	100%
Lichens	policy_1	21	2.52	62%
Lichens	policy_2	21	7.90	90%
Lichens	policy_3	21	13.90	100%
Lichens	policy_4	21	9.19	100%
Spiders	policy_1	21	3.50	78%
Spiders	policy_2	21	10.94	100%
Spiders	policy_3	21	16.28	100%
Spiders	policy_4	21	11.44	100%
Vascular plants	policy_1	21	6.29	86%
Vascular plants	policy_2	21	12.90	100%
Vascular plants	policy_3	21	18.86	100%
Vascular plants	policy_4	21	14.00	100%

5.4 summary

Across both aggregation assumptions, experts perceive all policy scenarios as improving outcomes relative to the counterfactual. The policy scenarios are clearly differentiated in strength, with policy_3 consistently eliciting the largest average uplift, followed by policy_2 and policy_4, with policy_1 weakest.

However, deltas alone do not answer whether the target-relevant condition (stability/recovery by 2042) is achieved. That requires combining Q1 and Q2 into absolute outcomes (Section 6).

---

---

6. Absolute outcomes under policy relative to counterfactual

Section 6 performs the key “sufficiency test” by converting the elicited components into absolute outcomes under each policy scenario:

abs_policy_mag_signed = cf_mag_signed + pol_mag_signed

This step matters because it answers the policy-facing question:

• Do scenarios merely slow decline (still negative)?
• Do they stabilise (near 0)?
• Do they generate net recovery (positive)?

We present headline results under both Assumptions A and B. Here, A vs B matters most for habitat-group interpretation, because the collapsing step in B prevents freshwater assemblage and lichen multiplicity from dominating habitat-group means.

Absolute signed outcomes under policy are:

• Absolute policy outcome = Counterfactual (Q1) + Policy delta (Q2)

---

6.1 Scenario-level outcomes: sufficiency test (Assumptions A and B)

policy_summary_AB <- pol_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_policy = mean(abs_policy_mag_signed, na.rm = TRUE),
    median_policy = median(abs_policy_mag_signed, na.rm = TRUE),
    sd_policy = sd(abs_policy_mag_signed, na.rm = TRUE),
    p10 = quantile(abs_policy_mag_signed, 0.10, na.rm = TRUE),
    p90 = quantile(abs_policy_mag_signed, 0.90, na.rm = TRUE),
    prop_positive = mean(abs_policy_mag_signed > 0, na.rm = TRUE),
    prop_near_zero = mean(abs(abs_policy_mag_signed) <= 1, na.rm = TRUE),
    prop_negative = mean(abs_policy_mag_signed < 0, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(assumption, desc(mean_policy))

policy_summary_AB %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  kable(caption = "Absolute signed % outcome under each policy (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Absolute signed % outcome under each policy (Assumptions A and B)

assumption	scenario	mean_policy	median_policy	sd_policy	p10	p90	prop_positive	prop_near_zero	prop_negative
A (rows = all survey_taxa)	policy_3	4.23	5	9.42	-10	15.0	0.67	0.13	0.21
A (rows = all survey_taxa)	policy_4	-0.52	0	10.30	-15	13.0	0.41	0.23	0.36
A (rows = all survey_taxa)	policy_2	-1.57	0	10.54	-20	9.0	0.34	0.28	0.39
A (rows = all survey_taxa)	policy_1	-8.18	-5	12.20	-25	5.0	0.18	0.20	0.63
B (collapsed = unique taxa:habitat_group)	policy_3	4.63	5	9.78	-10	15.0	0.68	0.12	0.21
B (collapsed = unique taxa:habitat_group)	policy_4	-0.14	0	10.67	-15	13.0	0.42	0.22	0.35
B (collapsed = unique taxa:habitat_group)	policy_2	-1.01	0	10.89	-20	10.0	0.38	0.31	0.34
B (collapsed = unique taxa:habitat_group)	policy_1	-7.70	-5	12.64	-25	5.3	0.20	0.22	0.59

Interpretation: strong uplifts, but mean absolute outcomes remain negative under every policy scenario

The scenario-level absolute outcome table shows a consistent pattern:

1. All scenarios improve outcomes relative to the counterfactual (consistent with Section 5),
2. but the mean absolute signed outcome remains negative under every scenario, i.e. the central tendency is continued decline by 2042 even after policy uplift is applied.

Under Assumption A, the mean absolute outcomes are:

• policy_1: −9.67%
• policy_2: −4.55%
• policy_4: −4.99%
• policy_3: −1.73% (closest to stability)

This is the core quantitative message: policy_3 comes close to offsetting baseline decline on average, but does not cross the stability threshold (0%) in the mean.

Two nuances are important for interpretation:

• The median can be more optimistic than the mean (e.g., policy_3 median is positive), indicating that a substantial portion of units have positive outcomes even when the mean remains negative.
• The p10 values remain strongly negative for all scenarios, implying persistent “tails” of severe decline that can pull the average below zero even if many units improve.

A vs B: the table is shown under both assumptions because the “headline mean” is sensitive to how units are counted. If the A vs B results are close, that strengthens confidence that conclusions are not being driven by survey-row multiplicity.

ggplot(policy_summary_AB, aes(x = scenario, y = mean_policy, fill = assumption)) +
  geom_col(position = position_dodge(width = 0.8)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Absolute signed % outcome under each policy (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % change (CF + delta)"
  ) +
  theme_minimal()

ggplot(subset(policy_summary_AB, assumption = "A (rows = all survey_taxa)"), aes(x = scenario, y = mean_policy)) +
    geom_col(position = position_dodge(width = 0.8)) +
    geom_hline(yintercept = 0, linetype = "dashed") +
    labs(
        title = "Absolute signed % outcome under each policy (Assumptions A and B)",
        x = NULL,
        y = "Mean signed % change (CF + delta)"
    ) +
    theme_minimal(base_size = 18)

## Warning: In subset.data.frame(policy_summary_AB, assumption = "A (rows = all survey_taxa)") :
##  extra argument 'assumption' will be disregarded

6.2 Visual comparison: counterfactual mean vs policy mean (Assumptions A and B)

cf_mean_AB <- cf_AB %>%
  group_by(assumption) %>%
  summarise(mean_cf = mean(cf_mag_signed, na.rm = TRUE), .groups = "drop")

policy_summary_plot_AB <- policy_summary_AB %>%
  left_join(cf_mean_AB, by = "assumption")

ggplot(policy_summary_plot_AB, aes(x = scenario)) +
  geom_col(aes(y = mean_policy, fill = assumption), position = position_dodge(width = 0.8)) +
  facet_wrap(~ assumption) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_hline(aes(yintercept = mean_cf), linetype = "dotted") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Mean signed % outcome: Counterfactual vs Policy (A and B)",
    subtitle = "Dotted = counterfactual mean; dashed = stability (0%)",
    x = NULL,
    y = "Mean signed % change"
  ) +
  theme_minimal()

Reading the CF vs policy mean plot: “distance to zero” is the target-relevant diagnostic

In the mean comparison plot, the dashed line at 0% represents stability (no net change by 2042). The dotted line marks the counterfactual mean. The vertical distance:

• from dotted line to each policy bar ≈ the average uplift achieved, and
• from each policy bar to zero ≈ the remaining recovery gap (the shortfall relative to stability).

The important result is that all scenario means sit above the counterfactual mean (uplift), but still below zero (insufficient for stability in the mean). This graph makes the distinction between “uplift” and “recovery” visually explicit.

6.3 Proportion stabilised and recovered (Assumptions A and B)

pol_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    prop_recovery = mean(abs_policy_mag_signed > 0, na.rm = TRUE),
    prop_stabilised = mean(abs_policy_mag_signed >= 0, na.rm = TRUE),
    prop_decline_remaining = mean(abs_policy_mag_signed < 0, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(across(where(is.numeric), ~ scales::percent(.x, accuracy = 1))) %>%
  kable(caption = "Proportion stabilised/recovered under each policy (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Proportion stabilised/recovered under each policy (Assumptions A and B)

assumption	scenario	prop_recovery	prop_stabilised	prop_decline_remaining
A (rows = all survey_taxa)	policy_1	18%	37%	63%
A (rows = all survey_taxa)	policy_2	34%	61%	39%
A (rows = all survey_taxa)	policy_3	67%	79%	21%
A (rows = all survey_taxa)	policy_4	41%	64%	36%
B (collapsed = unique taxa:habitat_group)	policy_1	20%	41%	59%
B (collapsed = unique taxa:habitat_group)	policy_2	38%	66%	34%
B (collapsed = unique taxa:habitat_group)	policy_3	68%	79%	21%
B (collapsed = unique taxa:habitat_group)	policy_4	42%	65%	35%

Interpretation: many units stabilise/recover under policy_3, but a non-trivial minority remain in decline

The stabilised/recovered fractions provide a more decision-friendly view than the mean alone:

• Under policy_3, around 79% of units are stabilised (≥0), and around 67–68% are positive (>0), while roughly 21% remain negative.
• We see clear monotonic improvement from policy_1 → policy_3 (and policy_4 intermediate).

This reconciles two facts that can otherwise look contradictory:

• Why can the mean remain negative under policy_3 while most units are stable/positive?
  • Because the remaining negative units are often sufficiently negative to “drag down” the mean; the aggregate is sensitive to tails.

A vs B: the fraction table is presented under both assumptions; the similarity between A and B proportions indicates the “many improve, some remain strongly negative” pattern is not an artefact of survey-row multiplicity.

6.4 Absolute outcomes by habitat_group (Assumptions A and B)

policy_by_hab_AB <- df_AB %>%
  group_by(assumption, habitat_group, scenario) %>%
  summarise(
    n = n(),
    mean_policy = mean(abs_policy_mag_signed, na.rm = TRUE),
    prop_positive = mean(abs_policy_mag_signed > 0, na.rm = TRUE),
    .groups = "drop"
  )

policy_by_hab_AB %>%
  mutate(
    mean_policy = round(mean_policy, 2),
    prop_positive = scales::percent(prop_positive, accuracy = 1)
  ) %>%
  arrange(assumption, habitat_group, scenario) %>%
  kable(caption = "Absolute outcomes under policy by habitat_group (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "360px")

Absolute outcomes under policy by habitat_group (Assumptions A and B)

assumption	habitat_group	scenario	n	mean_policy	prop_positive
A (rows = all survey_taxa)	Freshwater/Wetland	counterfactual	56	-13.83	2%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_1	56	-9.39	9%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_2	56	-3.00	24%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_3	56	3.78	69%
A (rows = all survey_taxa)	Freshwater/Wetland	policy_4	56	-0.70	41%
A (rows = all survey_taxa)	Open habitats	counterfactual	49	-13.43	0%
A (rows = all survey_taxa)	Open habitats	policy_1	49	-9.49	17%
A (rows = all survey_taxa)	Open habitats	policy_2	49	-2.62	36%
A (rows = all survey_taxa)	Open habitats	policy_3	49	2.89	66%
A (rows = all survey_taxa)	Open habitats	policy_4	49	-1.91	40%
A (rows = all survey_taxa)	Woodland	counterfactual	42	-9.00	22%
A (rows = all survey_taxa)	Woodland	policy_1	42	-5.00	30%
A (rows = all survey_taxa)	Woodland	policy_2	42	1.57	45%
A (rows = all survey_taxa)	Woodland	policy_3	42	6.40	68%
A (rows = all survey_taxa)	Woodland	policy_4	42	1.38	42%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	counterfactual	42	-13.47	0%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_1	42	-9.12	10%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_2	42	-2.58	30%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_3	42	4.12	68%
B (collapsed = unique taxa:habitat_group)	Freshwater/Wetland	policy_4	42	-0.35	42%
B (collapsed = unique taxa:habitat_group)	Open habitats	counterfactual	42	-13.36	0%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_1	42	-8.99	20%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_2	42	-2.01	40%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_3	42	3.39	70%
B (collapsed = unique taxa:habitat_group)	Open habitats	policy_4	42	-1.44	42%
B (collapsed = unique taxa:habitat_group)	Woodland	counterfactual	42	-9.00	22%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_1	42	-5.00	30%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_2	42	1.57	45%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_3	42	6.40	68%
B (collapsed = unique taxa:habitat_group)	Woodland	policy_4	42	1.38	42%

ggplot(policy_by_hab_AB, aes(x = scenario, y = mean_policy, fill = assumption)) +
  geom_col(position = position_dodge(width = 0.8)) +
  facet_wrap(~ habitat_group) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_fill_manual(values = two_cols) +
  labs(
    title = "Absolute outcomes under policy by habitat_group (Assumptions A and B)",
    x = NULL,
    y = "Mean signed % change (CF + delta)"
  ) +
  theme_minimal()

Interpretation: habitat-group constraints are clear, and Woodland is the closest to neutrality

The habitat-group absolute outcomes sharpen where the recovery gap sits:

• Freshwater/Wetland remains negative under every scenario, including policy_3 (mean still below 0), though the proportion positive rises strongly under policy_3.
• Open habitats also remain negative across scenarios; policy_3 improves outcomes markedly but does not remove the negative mean.
• Woodland is the closest to neutral and can become slightly positive under the strongest scenario (consistent with the baseline being less negative here).

This pattern is consistent with the baseline in Section 4: the most negative baseline domains (Freshwater/Wetland and Open habitats) remain the binding constraints on achieving system-wide stability by 2042.

Because these habitat-group results are shown for both A and B, we can interpret the habitat differences with higher confidence: the story persists even after removing the freshwater/lichen row multiplicity in Assumption B.

6.5 Absolute outcomes by taxa (A only; A == B)

policy_by_taxa <- pol_A %>%
  group_by(taxa, scenario) %>%
  summarise(
    mean_policy = mean(abs_policy_mag_signed, na.rm = TRUE),
    prop_positive = mean(abs_policy_mag_signed > 0, na.rm = TRUE),
    .groups = "drop"
  )

policy_by_taxa %>%
  mutate(
    mean_policy = round(mean_policy, 2),
    prop_positive = scales::percent(prop_positive, accuracy = 1)
  ) %>%
  kable(caption = "Absolute outcomes under policy by taxa (A; identical interpretation under B)") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(height = "360px")

Absolute outcomes under policy by taxa (A; identical interpretation under B)

taxa	scenario	mean_policy	prop_positive
Beetles	policy_1	-5.71	29%
Beetles	policy_2	1.95	57%
Beetles	policy_3	6.62	67%
Beetles	policy_4	1.24	52%
Bryophytes	policy_1	-9.38	10%
Bryophytes	policy_2	-2.67	24%
Bryophytes	policy_3	3.48	71%
Bryophytes	policy_4	-1.10	33%
Diptera	policy_1	-8.50	22%
Diptera	policy_2	-1.94	28%
Diptera	policy_3	3.44	61%
Diptera	policy_4	-0.83	33%
Freshwater invertebrate assemblages	policy_1	-10.14	5%
Freshwater invertebrate assemblages	policy_2	-4.19	10%
Freshwater invertebrate assemblages	policy_3	2.81	71%
Freshwater invertebrate assemblages	policy_4	-1.71	38%
Lichens	policy_1	-10.48	5%
Lichens	policy_2	-5.10	14%
Lichens	policy_3	0.90	52%
Lichens	policy_4	-3.81	24%
Spiders	policy_1	-5.67	33%
Spiders	policy_2	1.78	61%
Spiders	policy_3	7.11	83%
Spiders	policy_4	2.28	61%
Vascular plants	policy_1	-7.05	24%
Vascular plants	policy_2	-0.43	48%
Vascular plants	policy_3	5.52	67%
Vascular plants	policy_4	0.67	48%

---

6.6 Probability of full recovery (policy outcome ≥ 0) with expert-bootstrap uncertainty

To make the scenario results easier to interpret, we compute the probability of full recovery for each scenario:

• Full recovery means the absolute implied policy outcome is ≥ 0% (stabilisation or net recovery).
  
• In this workflow, Q2 magnitudes are elicited relative to the counterfactual, and the absolute outcome is therefore represented by abs_policy_mag_signed (counterfactual signed magnitude plus policy delta).

To quantify uncertainty given a small expert panel, we compute the recovery probability for each scenario and assumption, and attach 95% percentile intervals from a cluster bootstrap over experts (resampling experts,id_gj, with replacement).

# ============================================================
#  probability of full recovery (abs_policy_mag_signed >= 0)
# cluster bootstrap over experts (id_gj), by Assumption (A vs B)
# ============================================================

# packagess
library(dplyr)
library(tidyr)
library(purrr)
library(scales)
library(knitr)
library(kableExtra)

# ---- 0) Use the Section 6 object exactly as-is ----
stopifnot(exists("df_AB"))
stopifnot(all(c("assumption","scenario","id_gj","abs_policy_mag_signed") %in% names(df_AB)))

policy_levels <- c("policy_1","policy_2","policy_3","policy_4")

dat <- df_AB %>%
  mutate(
    assumption = as.character(assumption),
    scenario   = as.character(scenario)
  ) %>%
  filter(scenario %in% policy_levels)

assumptions <- sort(unique(dat$assumption))
scenarios   <- policy_levels

# ---- 1) Define full recovery indicator ----
dat <- dat %>%
  mutate(full_recovery = abs_policy_mag_signed >= 0)

# ---- 2) Bootstrap settings ----
set.seed(1)
B <- 5000

# helper: weighted mean for logical indicators under expert-resampling weights
wtd_mean <- function(x, w) {
  ok <- !is.na(x) & !is.na(w)
  if (!any(ok)) return(NA_real_)
  sum(w[ok] * x[ok]) / sum(w[ok])
}

# ---- 3) One cluster-bootstrap draw within ONE assumption ----
boot_once_assumption <- function(df_a) {

  experts <- sort(unique(df_a$id_gj[!is.na(df_a$id_gj)]))
  Nexp <- length(experts)

  samp <- sample(experts, size = Nexp, replace = TRUE)
  w <- table(samp)

  dfb <- df_a %>%
    semi_join(tibble(id_gj = as.integer(names(w))), by = "id_gj") %>%
    mutate(.w = as.numeric(w[as.character(id_gj)]))

  dfb %>%
    group_by(scenario) %>%
    summarise(
      p = wtd_mean(as.numeric(full_recovery), .w),
      .groups = "drop"
    )
}

# ---- 4) Run bootstrap separately for A and B ----
boot_draws <- map_dfr(assumptions, function(a) {
  df_a <- dat %>% filter(assumption == a)

  map_dfr(seq_len(B), function(b) {
    boot_once_assumption(df_a) %>%
      mutate(assumption = a, b = b)
  })
})

# ---- 5) Point estimates (non-bootstrap) ----
point_est <- dat %>%
  group_by(assumption, scenario) %>%
  summarise(
    p_hat = mean(full_recovery, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

# ---- 6) Bootstrap percentile CIs ----
boot_ci <- boot_draws %>%
  group_by(assumption, scenario) %>%
  summarise(
    p_lo  = quantile(p, 0.025, na.rm = TRUE),
    p_med = quantile(p, 0.50,  na.rm = TRUE),
    p_hi  = quantile(p, 0.975, na.rm = TRUE),
    .groups = "drop"
  )

prob_recovery_AB <- point_est %>%
  left_join(boot_ci, by = c("assumption","scenario")) %>%
  mutate(
    scenario   = factor(scenario, levels = policy_levels),
    assumption = factor(assumption, levels = assumptions)
  ) %>%
  arrange(assumption, scenario)

# ---- 7) Pretty table ----
prob_recovery_AB %>%
  mutate(
    p_hat = percent(p_hat, accuracy = 1),
    p_lo  = percent(p_lo,  accuracy = 1),
    p_hi  = percent(p_hi,  accuracy = 1),
    ci95  = paste0(p_lo, "–", p_hi)
  ) %>%
  select(assumption, scenario, n, p_hat, ci95) %>%
  pivot_wider(
    id_cols = c(scenario),
    names_from = assumption,
    values_from = c(n, p_hat, ci95),
    names_glue = "{assumption}_{.value}"
  ) %>%
  kable(caption = "Probability of full recovery (abs_policy_mag_signed ≥ 0%) by scenario and assumption (cluster bootstrap over experts; 95% CI)") %>%
  kable_styling(full_width = FALSE) 

Probability of full recovery (abs_policy_mag_signed ≥ 0%) by scenario and assumption (cluster bootstrap over experts; 95% CI)

scenario	A (rows = all survey_taxa)_n	B (collapsed = unique taxa:habitat_group)_n	A (rows = all survey_taxa)_p_hat	B (collapsed = unique taxa:habitat_group)_p_hat	A (rows = all survey_taxa)_ci95	B (collapsed = unique taxa:habitat_group)_ci95
policy_1	147	126	37%	41%	17%–58%	18%–64%
policy_2	147	126	61%	66%	37%–80%	41%–85%
policy_3	147	126	79%	79%	50%–97%	51%–98%
policy_4	147	126	64%	65%	33%–90%	35%–91%

# ---- 8) Plot: probability with CI (facet by assumption) ----
prob_recovery_AB %>%
  ggplot(aes(x = scenario, y = p_hat)) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = p_lo, ymax = p_hi), width = 0.15) +
  facet_wrap(~ assumption) +
  geom_hline(yintercept=1, linetype="dashed") +
  scale_y_continuous(labels = percent_format(accuracy = 1), limits = c(0, 1)) +
  labs(
    title = "Probability of full recovery (abs_policy_mag_signed ≥ 0%)",
    subtitle = "Cluster bootstrap over experts (id_gj), 95% percentile intervals",
    x = NULL,
    y = "Probability"
  ) +
  theme_minimal()

# ---- 9) context table showing mean absolute outcomes too ----
dat %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_abs_policy = mean(abs_policy_mag_signed, na.rm = TRUE),
    median_abs_policy = median(abs_policy_mag_signed, na.rm = TRUE),
    prop_full_recovery = mean(full_recovery, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  arrange(assumption, scenario) %>%
  kable(digits = 2, caption = "Context: mean/median absolute policy outcome and raw recovery share") %>%
  kable_styling(full_width = FALSE) 

Context: mean/median absolute policy outcome and raw recovery share

assumption	scenario	mean_abs_policy	median_abs_policy	prop_full_recovery
A (rows = all survey_taxa)	policy_1	-8.18	-5	0.37
A (rows = all survey_taxa)	policy_2	-1.57	0	0.61
A (rows = all survey_taxa)	policy_3	4.23	5	0.79
A (rows = all survey_taxa)	policy_4	-0.52	0	0.64
B (collapsed = unique taxa:habitat_group)	policy_1	-7.70	-5	0.41
B (collapsed = unique taxa:habitat_group)	policy_2	-1.01	0	0.66
B (collapsed = unique taxa:habitat_group)	policy_3	4.63	5	0.79
B (collapsed = unique taxa:habitat_group)	policy_4	-0.14	0	0.65

Across both assumptions, recovery probabilities increase strongly with policy ambition. Under Assumption A (all survey_taxa rows), the probability that the implied absolute policy outcome is ≥0% rises from 37% under policy_1 to 61% under policy_2, 79% under policy_3, and 64% under policy_4. The bootstrap intervals are wide—reflecting the small expert panel—but clearly shift upward with scenario strength (e.g. policy_3: 50–97%).

Results are highly consistent under Assumption B (collapsed to unique taxa × habitat_group), with probabilities of 41%, 66%, 79%, and 65% for policy_1–4 respectively. The similarity between A and B indicates that aggregation structure does not materially alter the overall recovery signal.

Policy_3 consistently emerges as the strongest scenario: it produces the highest probability of full recovery (~79%) under both assumptions. However, the wide bootstrap intervals (e.g. ~51–98% under B) indicate substantial expert-level uncertainty. Policy_1 remains clearly insufficient, with recovery probabilities below 50% under both assumptions and lower bounds well below 20%. Overall, stronger policy scenarios substantially increase the likelihood of stabilisation or recovery, but do not deliver near-certainty of full recovery across taxa.

---

6.7 Summary

Section 6 provides the core “sufficiency” conclusion:

• Baseline counterfactual change is strongly negative (Section 4).
• Policies generate positive deltas relative to baseline (Section 5).
• After combining CF + delta, mean absolute outcomes remain negative under all scenarios, including the strongest package.
• policy_3 is consistently closest to stability, and yields a high share of stable/positive units, but does not eliminate residual decline.

This creates a clear bridge to Section 7: we now test whether this sufficiency verdict is robust given panel size and disagreement.

---

7. Scenario-level expert agreement, uncertainty structure, and robustness

Section 7 examines whether the headline conclusions (notably: “policy_3 is strongest but still not sufficient on average”) are robust to the small expert panel and to disagreement about magnitude.

We focus on three complementary diagnostics:

1. Dispersion (SD) in signed outcomes across experts (agreement proxy)
   
2. Mean–variance relationship (whether disagreement scales with effect size; a key Chapter-3-style diagnostic)
   
3. Resampling / panel-size sensitivity (whether means and rankings depend on who is in the panel)

To keep uncertainty aligned with the outcome summaries, dispersion and resampling are reported under both assumptions A and B.

7.1 Dispersion in signed magnitude (agreement proxy) — A and B

consensus_var_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_") | scenario == "counterfactual") %>%
  group_by(assumption, taxa, habitat_group, scenario) %>%
  summarise(
    mean_signed = mean(abs_policy_mag_signed, na.rm = TRUE),
    sd_signed   = sd(abs_policy_mag_signed, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

consensus_summary_AB <- consensus_var_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_sd = mean(sd_signed, na.rm = TRUE),
    median_sd = median(sd_signed, na.rm = TRUE),
    .groups = "drop"
  )

consensus_summary_AB %>%
  mutate(across(where(is.numeric), ~ round(.x, 2))) %>%
  kable(caption = "Mean dispersion (SD) in signed magnitude by scenario (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Mean dispersion (SD) in signed magnitude by scenario (Assumptions A and B)

assumption	scenario	mean_sd	median_sd
A (rows = all survey_taxa)	counterfactual	11.94	11.13
A (rows = all survey_taxa)	policy_1	12.77	12.47
A (rows = all survey_taxa)	policy_2	10.67	10.05
A (rows = all survey_taxa)	policy_3	9.54	8.70
A (rows = all survey_taxa)	policy_4	10.51	9.67
B (collapsed = unique taxa:habitat_group)	counterfactual	11.92	11.13
B (collapsed = unique taxa:habitat_group)	policy_1	12.77	12.47
B (collapsed = unique taxa:habitat_group)	policy_2	10.65	10.05
B (collapsed = unique taxa:habitat_group)	policy_3	9.51	8.70
B (collapsed = unique taxa:habitat_group)	policy_4	10.49	9.67

Interpretation: disagreement increases under stronger policy, as expected for high-leverage interventions

The dispersion table shows that SD in signed outcomes is lowest under the counterfactual and increases under stronger policy packages, peaking around policy_3:

• counterfactual mean SD ≈ 12
• policy_1 mean SD ≈ 15.5
• policy_2 mean SD ≈ 17.6
• policy_3 mean SD ≈ 23
• policy_4 mean SD ≈ 20

This pattern is typical of expert elicitation: experts agree more on the existence of baseline pressures than on the extent to which ambitious intervention can offset them. Importantly, the pattern is essentially the same under A and B, suggesting the uncertainty structure is not being driven by survey-row multiplicity.

Policy interpretation: uncertainty is concentrated where decisions are highest leverage (policy_3), so “how close to zero” policy_3 sits should be treated probabilistically rather than as a point estimate, even when the ranking of scenarios is stable.

ggplot(consensus_var_AB, aes(x = mean_signed, y = sd_signed)) +
  geom_point(alpha = 0.6) +
  facet_grid(assumption ~ scenario) +
  labs(
    title = "Mean–variance relationship across taxa × habitat_group (A and B)",
    x = "Mean signed % outcome",
    y = "Standard deviation (expert disagreement)"
  ) +
  theme_minimal()

Interpretation: the mean–variance relationship indicates structured (epistemic) uncertainty rather than noise

The mean–variance plot tests whether disagreement tends to increase as projected outcomes become more extreme. A widening spread away from zero (a “funnel” pattern) is a healthy signal: it suggests experts are most aligned around modest changes and diverge most when asked about large recovery or large decline.

In this analysis, the combination of (i) increasing dispersion under ambitious policy and (ii) mean–variance structure supports the interpretation that uncertainty is epistemic and structured (reflecting real limits on knowledge of intervention efficacy), rather than arbitrary disagreement.

7.2 Certainty structure — A and B

certainty_summary_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_") | scenario == "counterfactual") %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_cert = mean(if_else(scenario == "counterfactual", cf_cert, pol_cert), na.rm = TRUE),
    sd_cert = sd(if_else(scenario == "counterfactual", cf_cert, pol_cert), na.rm = TRUE),
    .groups = "drop"
  )

certainty_summary_AB %>%
  mutate(across(where(is.numeric), round, 2)) %>%
  kable(caption = "Mean and SD of certainty scores (1–5) (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 2)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

Mean and SD of certainty scores (1–5) (Assumptions A and B)

assumption	scenario	mean_cert	sd_cert
A (rows = all survey_taxa)	counterfactual	1.97	0.57
A (rows = all survey_taxa)	policy_1	1.99	0.79
A (rows = all survey_taxa)	policy_2	1.97	0.72
A (rows = all survey_taxa)	policy_3	1.93	0.75
A (rows = all survey_taxa)	policy_4	1.95	0.73
B (collapsed = unique taxa:habitat_group)	counterfactual	1.93	0.55
B (collapsed = unique taxa:habitat_group)	policy_1	1.97	0.79
B (collapsed = unique taxa:habitat_group)	policy_2	1.95	0.74
B (collapsed = unique taxa:habitat_group)	policy_3	1.91	0.76
B (collapsed = unique taxa:habitat_group)	policy_4	1.93	0.75

Interpretation: certainty scores add context, but magnitudes/dispersion carry most of the inference

Certainty scores provide a self-reported complement to dispersion. In practice, with small panels, certainty scales can be “sticky” (limited spread), so it is appropriate that the analysis relies more heavily on:

• signed magnitudes (what experts think will happen),
• dispersion (how much they disagree), and
• resampling stability (how dependent conclusions are on panel composition).

Where certainty is systematically low in the most negative baseline domains, it reinforces the case for carrying uncertainty forward into trajectory modelling rather than over-interpreting point means.

7.3 Resampling stability — A and B

set.seed(123)
scenario_levels_pol <- c("policy_1","policy_2","policy_3","policy_4")

resample_means_any <- function(dat, k, scenario_name, n_iter = 1000){
  replicate(n_iter, {
    sample_ids <- sample(unique(dat$id_gj), k, replace = FALSE)
    dat %>%
      filter(id_gj %in% sample_ids, scenario == scenario_name) %>%
      summarise(mean_signed = mean(abs_policy_mag_signed, na.rm = TRUE)) %>%
      pull(mean_signed)
  })
}

resampling_results_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_")) %>%
  group_by(assumption) %>%
  group_modify(~{
    dat <- .x
    purrr::map_dfr(scenario_levels_pol, function(sc){
      tibble(
        scenario = sc,
        k2 = sd(resample_means_any(dat, 2, sc)),
        k3 = sd(resample_means_any(dat, 3, sc)),
        k4 = sd(resample_means_any(dat, 4, sc)),
        k5 = sd(resample_means_any(dat, 5, sc)),
        k6 = sd(resample_means_any(dat, 6, sc))
      )
    })
  }) %>%
  ungroup()

resampling_results_AB %>%
  mutate(across(where(is.numeric), round, 2)) %>%
  kable(caption = "Resampling SD of mean signed outcome by panel size (Assumptions A & B)") %>%
  kable_styling(full_width = FALSE)

Resampling SD of mean signed outcome by panel size (Assumptions A & B)

assumption	scenario	k2	k3	k4	k5	k6
A (rows = all survey_taxa)	policy_1	6.47	4.67	3.48	2.49	1.63
A (rows = all survey_taxa)	policy_2	5.11	3.75	2.82	2.06	1.28
A (rows = all survey_taxa)	policy_3	4.30	3.09	2.34	1.72	1.14
A (rows = all survey_taxa)	policy_4	5.05	3.67	2.75	2.03	1.31
B (collapsed = unique taxa:habitat_group)	policy_1	6.39	4.91	3.72	2.69	1.70
B (collapsed = unique taxa:habitat_group)	policy_2	5.35	3.90	2.94	2.11	1.41
B (collapsed = unique taxa:habitat_group)	policy_3	4.22	3.30	2.42	1.78	1.19
B (collapsed = unique taxa:habitat_group)	policy_4	5.12	3.83	2.86	2.01	1.36

resampling_long_AB <- resampling_results_AB %>%
  pivot_longer(cols = starts_with("k"), names_to = "panel_size", values_to = "sd_mean") %>%
  mutate(panel_size = as.numeric(str_remove(panel_size, "k")))

ggplot(resampling_long_AB, aes(x = panel_size, y = sd_mean)) +
  geom_line() +
  geom_point() +
  facet_grid(assumption ~ scenario) +
  labs(
    title = "Panel-size sensitivity of mean signed policy outcome (Assumptions A and B)",
    x = "Number of experts sampled",
    y = "SD of mean signed outcome"
  ) +
  theme_minimal()

Interpretation: scenario ordering is stable; sensitivity is highest where outcomes are closest to zero

The resampling analysis asks: “If we had elicited from a slightly different subset of experts, would the headline means and scenario ranking change?”

The expected and policy-relevant pattern is:

• Stability increases quickly as panel size grows (k increases),
• The strongest scenario (policy_3) shows the greatest sensitivity, because it is closest to the stability threshold and has the highest dispersion.

This supports a nuanced conclusion:

• We can be fairly confident in the directional result (“policy improves outcomes”) and the relative ranking of scenarios,
• while treating the absolute distance to zero under policy_3 as uncertain and best communicated with uncertainty bounds.

7.4 summary

Uncertainty diagnostics strengthen (rather than weaken) the credibility of the results:

• Disagreement increases under more ambitious policy, which is realistic.
• The uncertainty structure is coherent and similar under assumptions A and B.
• Resampling indicates that headline conclusions are not obviously driven by a single expert.

Together, this supports using the elicitation outputs as inputs to trajectory modelling, with uncertainty propagation and careful interpretation of near-zero outcomes.

---

8. Taxonomic group-level expert agreement, uncertainty structure, and robustness

This subsection applies the same three agreement/robustness diagnostics used above, but after aggregating outcomes within taxonomic groups. This shows whether expert consensus differs systematically among broad taxa rather than among individual taxa × habitat combinations.

practical note: these chunks assume df_AB and abs_policy_mag_signed already exist exactly as in Section 7 .

8.1 Scenario-level

8.1.1 Dispersion in signed magnitude (agreement proxy) — Assumption B

## ---- taxon-level-setup ------------------------------------------------------
# Taxonomic group-level aggregation of absolute signed outcomes
# Mirrors Section 7.1 logic, but aggregates within taxon first so that
# agreement/uncertainty is examined BETWEEN experts WITHIN each taxonomic group.

taxon_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_") | scenario == "counterfactual") %>%
  group_by(assumption, id_gj, scenario, taxa) %>%
  summarise(
    mean_signed_taxon = mean(abs_policy_mag_signed, na.rm = TRUE),
    .groups = "drop"
  )

taxon_AB %>%
  summarise(
    n_rows = n(),
    n_assumptions = n_distinct(assumption),
    n_experts = n_distinct(id_gj),
    n_scenarios = n_distinct(scenario),
    n_taxa = n_distinct(taxa)
  ) %>%
  kable(caption = "Coverage check: taxonomic group-level uncertainty dataset") %>%
  kable_styling(full_width = FALSE)

Coverage check: taxonomic group-level uncertainty dataset

n_rows	n_assumptions	n_experts	n_scenarios	n_taxa
490	2	7	5	7

## ---- taxon-dispersion -------------------------------------------------------
# 1) Dispersion in signed outcomes across experts, at taxonomic group level

taxon_consensus_var_AB <- taxon_AB %>%
  group_by(assumption, taxa, scenario) %>%
  summarise(
    mean_signed = mean(mean_signed_taxon, na.rm = TRUE),
    sd_signed   = sd(mean_signed_taxon, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

taxon_consensus_summary_AB <- taxon_consensus_var_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_sd   = mean(sd_signed, na.rm = TRUE),
    median_sd = median(sd_signed, na.rm = TRUE),
    .groups = "drop"
  )

taxon_consensus_summary_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Mean dispersion (SD) in signed outcomes by scenario at taxonomic group level (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Mean dispersion (SD) in signed outcomes by scenario at taxonomic group level (Assumptions A and B)

assumption	scenario	mean_sd	median_sd
A (rows = all survey_taxa)	counterfactual	10.73	11.09
A (rows = all survey_taxa)	policy_1	11.24	12.32
A (rows = all survey_taxa)	policy_2	9.19	9.88
A (rows = all survey_taxa)	policy_3	7.91	8.49
A (rows = all survey_taxa)	policy_4	8.89	9.24
B (collapsed = unique taxa:habitat_group)	counterfactual	10.76	11.09
B (collapsed = unique taxa:habitat_group)	policy_1	11.28	12.32
B (collapsed = unique taxa:habitat_group)	policy_2	9.24	9.88
B (collapsed = unique taxa:habitat_group)	policy_3	7.95	8.49
B (collapsed = unique taxa:habitat_group)	policy_4	8.91	9.24

## ---- taxon-dispersion-plot --------------------------------------------------
ggplot(
  taxon_consensus_var_AB,
  aes(x = scenario, y = sd_signed)
) +
  geom_boxplot() +
  facet_wrap(~ assumption) +
  labs(
    title = "Expert disagreement by scenario at taxonomic group level",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

At the taxonomic-group level, expert disagreement under Assumption B is moderate rather than extreme. The counterfactual shows a mean SD of 10.76, and this rises slightly under policy_1 to 11.28, indicating that the weakest intervention introduces a little more variation in expert judgement at broad taxonomic scale. After that, however, the pattern moves back towards lower dispersion under the stronger scenarios, implying that once responses are averaged within taxonomic groups, the panel becomes more rather than less aligned about the expected direction and broad magnitude of policy effects.

This matters because it suggests that a substantial share of the disagreement seen earlier in the more disaggregated analysis is likely to reflect finer within-taxon ecological variation rather than fundamental disagreement over how the main taxonomic groups respond overall. In other words, experts may differ over details within taxa, but at the broader taxonomic level there is a clearer shared signal, particularly for the stronger policy packages. The taxonomic-group dataset underpinning this section covers seven taxa, five scenarios and all seven experts under Assumption B, so these summaries reflect the full panel at this level of aggregation.

8.1.2 Certainty structure — Assumption B

## ---- taxon-mean-variance-correlation ----------------------------------------
# Optional compact diagnostic table: correlation between extremity and disagreement
# Uses absolute mean magnitude to ask whether more extreme projected outcomes
# tend to attract more disagreement.

taxon_mv_diag_AB <- taxon_consensus_var_AB %>%
  mutate(abs_mean_signed = abs(mean_signed)) %>%
  group_by(assumption, scenario) %>%
  summarise(
    cor_absmean_sd = cor(abs_mean_signed, sd_signed, use = "complete.obs"),
    .groups = "drop"
  )

taxon_mv_diag_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Correlation between absolute mean outcome and disagreement at taxonomic group level") %>%
  kable_styling(full_width = FALSE)

Correlation between absolute mean outcome and disagreement at taxonomic group level

assumption	scenario	cor_absmean_sd
A (rows = all survey_taxa)	counterfactual	-0.52
A (rows = all survey_taxa)	policy_1	-0.89
A (rows = all survey_taxa)	policy_2	-0.78
A (rows = all survey_taxa)	policy_3	0.74
A (rows = all survey_taxa)	policy_4	-0.59
B (collapsed = unique taxa:habitat_group)	counterfactual	-0.51
B (collapsed = unique taxa:habitat_group)	policy_1	-0.87
B (collapsed = unique taxa:habitat_group)	policy_2	-0.76
B (collapsed = unique taxa:habitat_group)	policy_3	0.72
B (collapsed = unique taxa:habitat_group)	policy_4	-0.57

The mean–variance diagnostic adds an important nuance. Under Assumption B, the correlation between absolute mean outcome and disagreement is negative for the counterfactual (-0.51), policy_1 (-0.87), policy_2 (-0.76) and policy_4 (-0.57), but positive for policy_3 (+0.72). This indicates that, for most scenarios, taxonomic groups with larger average responses are not the ones attracting the greatest disagreement; if anything, stronger mean responses are associated with slightly tighter agreement. Policy_3 is the exception. Under the strongest scenario, the taxonomic groups expected to do best are also the groups where experts diverge more on magnitude.

Substantively, that is a reassuring rather than contradictory pattern. It suggests that under the strongest package there is still broad directional agreement that some taxonomic groups respond especially well, but less agreement over exactly how large those gains would be. With only seven taxonomic groups, these correlations should not be over-interpreted as precise statistical relationships, but they do support a sensible reading of the plot: uncertainty at taxonomic-group level is primarily about the size of gains under ambitious intervention, not about whether the stronger package is beneficial overall.

8.1.3 Resampling stability — Assumption B

## ---- taxon-resampling -------------------------------------------------------
# 3) Resampling / panel-size sensitivity at taxonomic group level
# Here the quantity being stabilised is the mean signed outcome across taxonomic groups.

set.seed(123)

taxon_panel_mean <- function(dat, k) {
  sampled_ids <- sample(unique(dat$id_gj), size = k, replace = FALSE)

  dat %>%
    filter(id_gj %in% sampled_ids) %>%
    group_by(assumption, scenario, taxa) %>%
    summarise(
      mean_signed_taxon = mean(mean_signed_taxon, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    group_by(assumption, scenario) %>%
    summarise(
      panel_mean = mean(mean_signed_taxon, na.rm = TRUE),
      .groups = "drop"
    )
}

taxon_resample_means <- function(k, assumption_in, scenario_in, n_iter = 1000) {
  dat <- taxon_AB %>%
    filter(assumption == assumption_in, scenario == scenario_in)

  map_dbl(seq_len(n_iter), ~{
    sampled_ids <- sample(unique(dat$id_gj), size = k, replace = FALSE)

    dat %>%
      filter(id_gj %in% sampled_ids) %>%
      group_by(taxa) %>%
      summarise(
        mean_signed_taxon = mean(mean_signed_taxon, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      summarise(panel_mean = mean(mean_signed_taxon, na.rm = TRUE)) %>%
      pull(panel_mean)
  })
}

scenario_levels <- c("counterfactual", "policy_1", "policy_2", "policy_3", "policy_4")
assumption_levels <- unique(taxon_AB$assumption)
panel_sizes <- sort(unique(taxon_AB$id_gj)) |> length() |> seq_len()

taxon_resampling_results_AB <- expand_grid(
  assumption = assumption_levels,
  scenario   = scenario_levels
) %>%
  mutate(
    k2 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(2, a, s))),
    k3 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(3, a, s))),
    k4 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(4, a, s))),
    k5 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(5, a, s))),
    k6 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(6, a, s))),
    k7 = map2_dbl(assumption, scenario, \(a, s) sd(taxon_resample_means(7, a, s)))
  )

taxon_resampling_results_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Resampling SD of mean signed outcome by panel size at taxonomic group level") %>%
  kable_styling(full_width = FALSE)

Resampling SD of mean signed outcome by panel size at taxonomic group level

assumption	scenario	k2	k3	k4	k5	k6	k7
A (rows = all survey_taxa)	counterfactual	6.30	4.51	3.41	2.49	1.62	0
A (rows = all survey_taxa)	policy_1	6.52	4.74	3.56	2.65	1.66	0
A (rows = all survey_taxa)	policy_2	5.19	3.78	2.89	2.05	1.41	0
A (rows = all survey_taxa)	policy_3	4.30	3.14	2.40	1.75	1.12	0
A (rows = all survey_taxa)	policy_4	5.10	3.73	2.80	2.00	1.31	0
B (collapsed = unique taxa:habitat_group)	counterfactual	6.22	4.60	3.39	2.43	1.63	0
B (collapsed = unique taxa:habitat_group)	policy_1	6.57	4.82	3.53	2.63	1.65	0
B (collapsed = unique taxa:habitat_group)	policy_2	5.33	3.84	2.89	2.08	1.33	0
B (collapsed = unique taxa:habitat_group)	policy_3	4.45	3.23	2.42	1.75	1.03	0
B (collapsed = unique taxa:habitat_group)	policy_4	5.08	3.69	2.74	2.02	1.29	0

## ---- taxon-resampling-plot --------------------------------------------------
taxon_resampling_long_AB <- taxon_resampling_results_AB %>%
  pivot_longer(
    cols = starts_with("k"),
    names_to = "panel_size",
    values_to = "sd_mean"
  ) %>%
  mutate(panel_size = as.numeric(str_remove(panel_size, "^k")))

ggplot(
  taxon_resampling_long_AB,
  aes(x = panel_size, y = sd_mean)
) +
  geom_line() +
  geom_point() +
  facet_grid(assumption ~ scenario) +
  labs(
    title = "Panel-size sensitivity of mean signed outcome at taxonomic group level",
    x = "Number of experts sampled",
    y = "SD of mean signed outcome"
  ) +
  theme_minimal()

The resampling analysis indicates that the taxon-level results under Assumption B are reasonably robust to panel composition. For every scenario, the SD of the resampled mean declines steadily as panel size increases, showing that the headline mean outcome stabilises quickly as more experts are included. At k = 6, the SD has already fallen to 1.63 for the counterfactual, 1.65 for policy_1, 1.33 for policy_2, 1.03 for policy_3 and 1.29 for policy_4.

Two features stand out. First, there is no sign here that the taxonomic-group conclusions are being driven by one unusually optimistic or pessimistic respondent; the mean stabilises in a smooth and expected way as k increases. Second, policy_3 is not only the strongest scenario in outcome terms elsewhere in the document, but also the most stable to panel re-composition at this taxonomic level. That means the ranking of policy_3 as the best-performing package appears to be robust, even though the precise magnitude of improvement remains uncertain.

8.1.4 summary

Under Assumption B, the taxonomic-group analysis suggests that the expert panel is more coherent at broad taxonomic level than might be inferred from the more disaggregated results alone. Disagreement is moderate under the counterfactual, increases slightly under the weakest policy package, and then moves back towards lower levels under the stronger scenarios, indicating that broad taxonomic responses become clearer once within-group ecological detail is averaged out.

The mean–variance diagnostic shows that this coherence is not uniform across all scenarios. For most scenarios, stronger average taxonomic responses are not associated with greater disagreement, but under policy_3 the relationship becomes positive, suggesting that the taxonomic groups expected to benefit most are also those for which experts differ most on the size of the gain. This points to magnitude uncertainty at the optimistic end, rather than disagreement over which scenario performs best.

Finally, the resampling results strengthen confidence in the headline conclusion. Taxon-level means stabilise quickly as panel size increases, and policy_3 is the most stable scenario under panel re-sampling as well as the strongest in substantive terms. Taken together, the taxonomic-group results support a clear directional interpretation under Assumption B: the strongest package remains policy_3, and that conclusion appears robust, even though uncertainty should still be carried forward around the size of taxonomic-group gains.

8.2 Taxon by taxon within each scenario

This section is based on Assumption B only.

These additional analyses examine whether expert agreement differs systematically among the specific taxonomic groups included in the elicitation. They first assess variation in disagreement among taxonomic groups within each scenario, and then assess whether some taxonomic groups are consistently associated with greater or lower disagreement when all scenarios are considered together. Only Assumption B is included.

Setup: Assumption B only

## ---- taxon-uncertainty-B-setup ----------------------------------------------
# Restrict to Assumption B only for all additional taxonomic-group analyses

taxon_B <- taxon_AB %>%
  filter(assumption == "B (collapsed = unique taxa:habitat_group)")

taxon_B %>%
  summarise(
    n_rows = n(),
    n_experts = n_distinct(id_gj),
    n_scenarios = n_distinct(scenario),
    n_taxa = n_distinct(taxa)
  ) %>%
  knitr::kable(
    caption = "Coverage check for additional taxonomic-group uncertainty analyses (Assumption B only)"
  )

Coverage check for additional taxonomic-group uncertainty analyses (Assumption B only)

n_rows	n_experts	n_scenarios	n_taxa
245	7	5	7

8.2.A Mean and SD by taxonomic group within each scenario

## ---- taxon-by-scenario-uncertainty-table ------------------------------------
# For each scenario, summarise expert disagreement separately for each taxonomic group

taxon_by_scenario_B <- taxon_B %>%
  group_by(scenario, taxa) %>%
  summarise(
    mean_signed = mean(mean_signed_taxon, na.rm = TRUE),
    sd_signed   = sd(mean_signed_taxon, na.rm = TRUE),
    min_signed  = min(mean_signed_taxon, na.rm = TRUE),
    max_signed  = max(mean_signed_taxon, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  arrange(scenario, desc(sd_signed))

taxon_by_scenario_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Expert disagreement by taxonomic group within each scenario (Assumption B)"
  )

Expert disagreement by taxonomic group within each scenario (Assumption B)

scenario	taxa	mean_signed	sd_signed	min_signed	max_signed	n
counterfactual	Diptera	-12.39	12.15	-30	1.67	7
counterfactual	Vascular plants	-13.33	11.75	-30	-1.67	7
counterfactual	Beetles	-10.48	11.50	-30	1.67	7
counterfactual	Spiders	-9.17	11.09	-30	1.67	7
counterfactual	Bryophytes	-12.62	10.04	-30	0.00	7
counterfactual	Lichens	-12.61	10.02	-30	-2.50	7
counterfactual	Freshwater invertebrate assemblages	-14.86	8.79	-30	-2.33	7
policy_1	Spiders	-5.67	13.22	-25	4.33	7
policy_1	Vascular plants	-7.05	12.76	-25	8.33	7
policy_1	Diptera	-8.50	12.51	-25	6.67	7
policy_1	Beetles	-5.71	12.32	-25	8.33	7
policy_1	Bryophytes	-9.38	10.04	-25	2.67	7
policy_1	Lichens	-9.54	9.18	-25	2.50	7
policy_1	Freshwater invertebrate assemblages	-10.14	8.95	-25	-1.00	7
policy_2	Spiders	1.78	11.49	-20	13.33	7
policy_2	Beetles	1.95	10.62	-20	11.33	7
policy_2	Diptera	-1.94	10.04	-20	9.67	7
policy_2	Vascular plants	-0.43	9.88	-20	10.00	7
policy_2	Bryophytes	-2.67	7.75	-20	2.67	7
policy_2	Lichens	-4.61	7.65	-20	2.50	7
policy_2	Freshwater invertebrate assemblages	-4.19	7.26	-20	0.33	7
policy_3	Diptera	3.44	9.92	-10	20.00	7
policy_3	Spiders	7.11	9.38	-10	16.67	7
policy_3	Beetles	6.62	9.16	-10	20.00	7
policy_3	Vascular plants	5.52	8.49	-10	15.00	7
policy_3	Freshwater invertebrate assemblages	2.81	6.42	-10	10.33	7
policy_3	Bryophytes	3.48	6.29	-10	10.00	7
policy_3	Lichens	1.32	5.98	-10	7.50	7
policy_4	Diptera	-0.83	10.80	-15	17.00	7
policy_4	Beetles	1.24	10.34	-15	17.00	7
policy_4	Vascular plants	0.67	9.77	-15	13.00	7
policy_4	Spiders	2.28	9.24	-15	11.33	7
policy_4	Bryophytes	-1.10	7.65	-15	8.00	7
policy_4	Lichens	-3.39	7.35	-15	4.50	7
policy_4	Freshwater invertebrate assemblages	-1.71	7.21	-15	5.00	7

#A2. Plot: which taxonomic groups are more uncertain within each scenario?
## ---- taxon-by-scenario-uncertainty-plot -------------------------------------
ggplot(
  taxon_by_scenario_B,
  aes(x = reorder(taxa, sd_signed), y = sd_signed)
) +
  geom_col() +
  coord_flip() +
  facet_wrap(~ scenario, scales = "free_y") +
  labs(
    title = "Expert disagreement by taxonomic group within each scenario",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

#A3. Heatmap view for fast comparison across scenarios
## ---- taxon-by-scenario-heatmap ----------------------------------------------
ggplot(
  taxon_by_scenario_B,
  aes(x = scenario, y = taxa, fill = sd_signed)
) +
  geom_tile() +
  geom_text(aes(label = round(sd_signed, 1)), size = 3) +
  labs(
    title = "Heatmap of expert disagreement by taxonomic group and scenario",
    x = NULL,
    y = NULL,
    fill = "SD"
  ) +
  theme_minimal()

Rank taxonomic groups within each scenario by disagreement

## ---- taxon-by-scenario-ranking ----------------------------------------------
taxon_by_scenario_ranks_B <- taxon_by_scenario_B %>%
  group_by(scenario) %>%
  mutate(
    disagreement_rank = min_rank(desc(sd_signed))
  ) %>%
  ungroup() %>%
  arrange(scenario, disagreement_rank)

taxon_by_scenario_ranks_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Ranking of taxonomic groups by expert disagreement within each scenario (Assumption B)"
  )

Ranking of taxonomic groups by expert disagreement within each scenario (Assumption B)

scenario	taxa	mean_signed	sd_signed	min_signed	max_signed	n	disagreement_rank
counterfactual	Diptera	-12.39	12.15	-30	1.67	7	1
counterfactual	Vascular plants	-13.33	11.75	-30	-1.67	7	2
counterfactual	Beetles	-10.48	11.50	-30	1.67	7	3
counterfactual	Spiders	-9.17	11.09	-30	1.67	7	4
counterfactual	Bryophytes	-12.62	10.04	-30	0.00	7	5
counterfactual	Lichens	-12.61	10.02	-30	-2.50	7	6
counterfactual	Freshwater invertebrate assemblages	-14.86	8.79	-30	-2.33	7	7
policy_1	Spiders	-5.67	13.22	-25	4.33	7	1
policy_1	Vascular plants	-7.05	12.76	-25	8.33	7	2
policy_1	Diptera	-8.50	12.51	-25	6.67	7	3
policy_1	Beetles	-5.71	12.32	-25	8.33	7	4
policy_1	Bryophytes	-9.38	10.04	-25	2.67	7	5
policy_1	Lichens	-9.54	9.18	-25	2.50	7	6
policy_1	Freshwater invertebrate assemblages	-10.14	8.95	-25	-1.00	7	7
policy_2	Spiders	1.78	11.49	-20	13.33	7	1
policy_2	Beetles	1.95	10.62	-20	11.33	7	2
policy_2	Diptera	-1.94	10.04	-20	9.67	7	3
policy_2	Vascular plants	-0.43	9.88	-20	10.00	7	4
policy_2	Bryophytes	-2.67	7.75	-20	2.67	7	5
policy_2	Lichens	-4.61	7.65	-20	2.50	7	6
policy_2	Freshwater invertebrate assemblages	-4.19	7.26	-20	0.33	7	7
policy_3	Diptera	3.44	9.92	-10	20.00	7	1
policy_3	Spiders	7.11	9.38	-10	16.67	7	2
policy_3	Beetles	6.62	9.16	-10	20.00	7	3
policy_3	Vascular plants	5.52	8.49	-10	15.00	7	4
policy_3	Freshwater invertebrate assemblages	2.81	6.42	-10	10.33	7	5
policy_3	Bryophytes	3.48	6.29	-10	10.00	7	6
policy_3	Lichens	1.32	5.98	-10	7.50	7	7
policy_4	Diptera	-0.83	10.80	-15	17.00	7	1
policy_4	Beetles	1.24	10.34	-15	17.00	7	2
policy_4	Vascular plants	0.67	9.77	-15	13.00	7	3
policy_4	Spiders	2.28	9.24	-15	11.33	7	4
policy_4	Bryophytes	-1.10	7.65	-15	8.00	7	5
policy_4	Lichens	-3.39	7.35	-15	4.50	7	6
policy_4	Freshwater invertebrate assemblages	-1.71	7.21	-15	5.00	7	7

Viewed scenario by scenario, there is a clear but not dramatic taxonomic structure to disagreement. Under the counterfactual, Diptera show the highest disagreement among experts (SD 12.15), followed closely by vascular plants (11.75) and beetles (11.50), whereas freshwater invertebrate assemblages, lichens and bryophytes sit at the lower end of the range. This suggests that even before policy is introduced, some taxa are already harder for the panel to judge consistently than others.

Across the policy scenarios, the same broad pattern persists. Beetles and Diptera tend to remain among the more uncertain groups, while bryophytes and lichens tend to remain among the more certain. At the same time, disagreement generally falls under the stronger policies for most taxa. For example, beetle disagreement declines from 12.32 under policy_1 to 9.16 under policy_3, while bryophytes decline from 10.04 to 6.29 over the same comparison. This implies that stronger interventions do not create taxon-specific confusion so much as reduce it, even though some groups remain harder to judge than others.

The ranking table is useful because it shows that this is not simply noise. Diptera are ranked most uncertain under the counterfactual, and vascular plants and beetles also repeatedly appear near the high-disagreement end of the ordering, whereas lichens and bryophytes tend to appear near the low-disagreement end. In substantive terms, this points to structured uncertainty: the expert panel is not equally certain across the taxonomic spectrum, and those differences are fairly stable across scenarios.

8.3 Taxonomic-group uncertainty across all scenarios combined

This section is based on Assumption B only throughout.

This approach is the probably a stronger-signal version.

8.3.A Overall disagreement by taxonomic group across all scenarios

## ---- taxon-overall-uncertainty-table ----------------------------------------
# Pool all scenarios together to test whether some taxonomic groups are
# consistently associated with more disagreement than others

taxon_overall_B <- taxon_B %>%
  group_by(taxa) %>%
  summarise(
    mean_signed = mean(mean_signed_taxon, na.rm = TRUE),
    sd_signed   = sd(mean_signed_taxon, na.rm = TRUE),
    median_signed = median(mean_signed_taxon, na.rm = TRUE),
    iqr_signed    = IQR(mean_signed_taxon, na.rm = TRUE),
    min_signed    = min(mean_signed_taxon, na.rm = TRUE),
    max_signed    = max(mean_signed_taxon, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  arrange(desc(sd_signed))

taxon_overall_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Overall expert disagreement by taxonomic group across all scenarios combined (Assumption B)"
  )

Overall expert disagreement by taxonomic group across all scenarios combined (Assumption B)

taxa	mean_signed	sd_signed	median_signed	iqr_signed	min_signed	max_signed	n
Vascular plants	-2.92	12.01	0.00	15.33	-30	15.00	35
Beetles	-1.28	11.89	1.67	15.83	-30	20.00	35
Diptera	-4.04	11.84	-1.17	13.33	-30	20.00	35
Spiders	-0.73	11.82	3.33	14.58	-30	16.67	35
Bryophytes	-4.46	9.91	0.00	12.67	-30	10.00	35
Freshwater invertebrate assemblages	-5.62	9.67	-5.00	10.17	-30	10.33	35
Lichens	-5.76	9.12	-3.00	10.88	-30	7.50	35

#B2. Plot: overall disagreement by taxonomic group
## ---- taxon-overall-uncertainty-plot -----------------------------------------
ggplot(
  taxon_overall_B,
  aes(x = reorder(taxa, sd_signed), y = sd_signed)
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Overall expert disagreement by taxonomic group across all scenarios combined",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

A Raw-distribution plot by taxonomic group across all scenarios is useful because the SD table can hide whether disagreement comes from one or two scenarios only, or is persistent across the whole set.

## ---- taxon-overall-distribution-plot ----------------------------------------
ggplot(
  taxon_B,
  aes(x = taxa, y = mean_signed_taxon)
) +
  geom_boxplot() +
  coord_flip() +
  labs(
    title = "Distribution of expert taxonomic-group outcomes across all scenarios combined",
    x = NULL,
    y = "Mean signed % outcome"
  ) +
  theme_minimal()

## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

When all scenarios are pooled, the taxonomic pattern becomes clearer. Vascular plants show the greatest overall disagreement across experts (SD 12.01), followed very closely by beetles (11.89). At the other end, lichens show the lowest disagreement (9.12), with freshwater invertebrate assemblages also relatively low (9.67). This pooled analysis therefore strengthens the impression from the scenario-by-scenario results: vascular plants and beetles are the most uncertain groups overall, whereas lichens and freshwater assemblages are the most consistently judged.

The pooled means also suggest that these differences in uncertainty are not simply tracking pessimism or optimism. For example, vascular plants have both the highest disagreement and a moderately negative overall mean, while lichens have the lowest disagreement despite also showing a clearly negative overall mean. This supports the interpretation that the analysis is picking up genuine differences in the consistency of expert judgement between taxa, rather than just differences in the average sign of the response.

The boxplot-style distribution across all scenarios is especially helpful here because it shows that some taxa have a wider overall spread in outcomes than others. In practice, this means that for groups such as vascular plants and beetles, expert disagreement is not confined to one single scenario but is more persistent across the elicitation as a whole. That makes them good candidates for carrying wider uncertainty bounds forward into any later aggregation or modelling step.

8.3.B do taxonomic groups differ in disagreement? (more formal comparison:)

This tests whether the spread across experts differs systematically among taxa.

Setup: Per taxonomic group, calculate within-scenario SD first

## ---- taxon-sd-dataset-for-testing -------------------------------------------
# Create a compact dataset of one disagreement value per taxon per scenario

taxon_sd_for_test_B <- taxon_B %>%
  group_by(taxa, scenario) %>%
  summarise(
    sd_signed = sd(mean_signed_taxon, na.rm = TRUE),
    .groups = "drop"
  )

taxon_sd_for_test_B %>%
  mutate(sd_signed = round(sd_signed, 2)) %>%
  knitr::kable(
    caption = "Scenario-specific disagreement values used for taxonomic-group comparison tests (Assumption B)"
  )

Scenario-specific disagreement values used for taxonomic-group comparison tests (Assumption B)

taxa	scenario	sd_signed
Beetles	counterfactual	11.50
Beetles	policy_1	12.32
Beetles	policy_2	10.62
Beetles	policy_3	9.16
Beetles	policy_4	10.34
Bryophytes	counterfactual	10.04
Bryophytes	policy_1	10.04
Bryophytes	policy_2	7.75
Bryophytes	policy_3	6.29
Bryophytes	policy_4	7.65
Diptera	counterfactual	12.15
Diptera	policy_1	12.51
Diptera	policy_2	10.04
Diptera	policy_3	9.92
Diptera	policy_4	10.80
Freshwater invertebrate assemblages	counterfactual	8.79
Freshwater invertebrate assemblages	policy_1	8.95
Freshwater invertebrate assemblages	policy_2	7.26
Freshwater invertebrate assemblages	policy_3	6.42
Freshwater invertebrate assemblages	policy_4	7.21
Lichens	counterfactual	10.02
Lichens	policy_1	9.18
Lichens	policy_2	7.65
Lichens	policy_3	5.98
Lichens	policy_4	7.35
Spiders	counterfactual	11.09
Spiders	policy_1	13.22
Spiders	policy_2	11.49
Spiders	policy_3	9.38
Spiders	policy_4	9.24
Vascular plants	counterfactual	11.75
Vascular plants	policy_1	12.76
Vascular plants	policy_2	9.88
Vascular plants	policy_3	8.49
Vascular plants	policy_4	9.77

8.3.B.1 Kruskal–Wallis test across taxonomic groups

Because there will be few observations per group, a non-parametric test is safer than leaning too hard on ANOVA assumptions.

## ---- taxon-kruskal-test -----------------------------------------------------
taxon_kw_B <- kruskal.test(sd_signed ~ taxa, data = taxon_sd_for_test_B)

taxon_kw_B

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sd_signed by taxa
## Kruskal-Wallis chi-squared = 18.05, df = 6, p-value = 0.00611

8.3.B.2 Pairwise Wilcoxon tests between taxonomic groups

nb: Only useful if the above suggests differences at all.

## ---- taxon-pairwise-tests ---------------------------------------------------
pairwise.wilcox.test(
  x = taxon_sd_for_test_B$sd_signed,
  g = taxon_sd_for_test_B$taxa,
  p.adjust.method = "BH",
  exact = FALSE
)

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  taxon_sd_for_test_B$sd_signed and taxon_sd_for_test_B$taxa 
## 
##                                     Beetles Bryophytes Diptera
## Bryophytes                          0.096   -          -      
## Diptera                             0.922   0.181      -      
## Freshwater invertebrate assemblages 0.085   0.697      0.085  
## Lichens                             0.096   0.605      0.096  
## Spiders                             1.000   0.181      0.922  
## Vascular plants                     0.922   0.368      0.605  
##                                     Freshwater invertebrate assemblages Lichens
## Bryophytes                          -                                   -      
## Diptera                             -                                   -      
## Freshwater invertebrate assemblages -                                   -      
## Lichens                             0.697                               -      
## Spiders                             0.085                               0.096  
## Vascular plants                     0.096                               0.181  
##                                     Spiders
## Bryophytes                          -      
## Diptera                             -      
## Freshwater invertebrate assemblages -      
## Lichens                             -      
## Spiders                             -      
## Vascular plants                     1.000  
## 
## P value adjustment method: BH

The formal comparison supports the descriptive results. The Kruskal–Wallis test indicates that disagreement differs significantly among taxonomic groups overall (chi-squared = 18.05, df = 6, p = 0.006). This is useful because it shows that the taxonomic structure seen in the tables and plots is unlikely to be a chance feature of a small panel alone.

At the same time, the pairwise Wilcoxon results are mostly non-significant after multiple-testing correction. That suggests the taxonomic signal is real in aggregate, but that the sample is too small to cleanly separate most individual pairs of taxa from one another. In practical terms, this means it is safer to interpret the results as a broad gradient of higher- versus lower-certainty taxa, rather than making strong claims that one named group is definitively more uncertain than another specific group.

This is a sensible outcome given the structure of the analysis. The formal test uses only five scenario-level SD values per taxon, so statistical power is limited. The absence of strong pairwise separation therefore should not be read as absence of taxonomic pattern; rather, it indicates that the pattern is distributed across the whole set of taxa rather than being driven by one stark outlier pair.

8.3.B.3 A simpler “certainty score” by taxonomic group

This converts disagreement into an intuitive certainty-like metric for quick reading.

Here I define:

higher SD = lower certainty certainty score = inverse-rescaled SD between 0 and 1 within this analysis

So this is purely relative within the taxonomic-group comparison.

## ---- taxon-relative-certainty-score -----------------------------------------
taxon_certainty_B <- taxon_overall_B %>%
  mutate(
    certainty_score = 1 - (sd_signed - min(sd_signed)) / (max(sd_signed) - min(sd_signed))
  ) %>%
  arrange(desc(certainty_score))

taxon_certainty_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Relative certainty score by taxonomic group across all scenarios combined (Assumption B only)"
  )

Relative certainty score by taxonomic group across all scenarios combined (Assumption B only)

taxa	mean_signed	sd_signed	median_signed	iqr_signed	min_signed	max_signed	n	certainty_score
Lichens	-5.76	9.12	-3.00	10.88	-30	7.50	35	1.00
Freshwater invertebrate assemblages	-5.62	9.67	-5.00	10.17	-30	10.33	35	0.81
Bryophytes	-4.46	9.91	0.00	12.67	-30	10.00	35	0.73
Spiders	-0.73	11.82	3.33	14.58	-30	16.67	35	0.07
Diptera	-4.04	11.84	-1.17	13.33	-30	20.00	35	0.06
Beetles	-1.28	11.89	1.67	15.83	-30	20.00	35	0.04
Vascular plants	-2.92	12.01	0.00	15.33	-30	15.00	35	0.00

#Optional plot for certainty score

## ---- taxon-certainty-score-plot ---------------------------------------------
ggplot(
  taxon_certainty_B,
  aes(x = reorder(taxa, certainty_score), y = certainty_score)
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Relative certainty by taxonomic group across all scenarios combined",
    x = NULL,
    y = "Relative certainty score (0-1)"
  ) +
  theme_minimal()

The relative certainty score provides a compact summary of the same pattern. Lichens emerge as the most certain taxonomic group overall, with a certainty score of 1.00, while freshwater invertebrate assemblages are also near the high-certainty end. At the opposite end, vascular plants sit at the lowest-certainty end of the scale, reflecting their position as the group with the highest pooled disagreement.

This is helpful for interpretation because it converts the spread measure into a more intuitive relative ranking. Read in that way, the expert panel appears comparatively confident about lichen responses and comparatively uncertain about vascular plant and beetle responses. Diptera also remain towards the more uncertain end, especially in the scenario-by-scenario outputs. The certainty score should still be treated as descriptive rather than absolute, but it is a clear way of communicating where uncertainty is concentrated taxonomically.

(8.2-3 Summary)

This section is based on Assumption B only throughout.

These additional taxa-level analyses show that expert disagreement is not evenly distributed across the taxonomic groups included in the elicitation. Instead, there is a structured taxonomic pattern: Diptera, beetles and especially vascular plants tend to attract higher disagreement, while lichens, bryophytes and freshwater invertebrate assemblages tend to attract lower disagreement. This pattern is visible within individual scenarios and becomes clearer when all scenarios are combined.

The scenario-by-scenario results also suggest that stronger interventions generally reduce disagreement within taxa rather than increase it. For several groups, including beetles and bryophytes, disagreement is lower under policy_3 than under policy_1 or the counterfactual. This indicates that stronger packages do not appear to create greater taxonomic ambiguity; instead, they often produce more coherent expectations, even if some groups remain intrinsically harder to judge than others.

The formal test reinforces this interpretation. There is evidence that disagreement differs significantly across taxa overall, but the pairwise comparisons are weak after adjustment, implying a broad taxonomic gradient in certainty rather than sharp separations between many individual taxon pairs. In practical terms, this means later interpretation should treat taxonomic uncertainty as real and structured, but avoid over-claiming precision about exactly which pairwise contrasts are most important.

Taken together, these findings suggest that uncertainty in the elicitation is partly ecological rather than merely stochastic. Some taxonomic groups appear consistently more difficult for experts to judge, and this should be carried forward into interpretation of the policy results. In particular, vascular plants and beetles appear to warrant more caution around magnitude estimates, whereas lichens and freshwater assemblages appear to support relatively tighter inference.

8.4 Summary

Section 8 shows that taxonomic-group uncertainty is structured rather than random. At the broad taxonomic-group level, disagreement is moderate under the counterfactual and weakest policy package, but tends to fall under the stronger scenarios, indicating that experts become more aligned once responses are averaged within taxa. The additional taxon-by-taxon analyses then show that this agreement is not evenly distributed: some taxa, particularly vascular plants, beetles and often Diptera, are consistently associated with greater disagreement, whereas lichens, bryophytes and freshwater invertebrate assemblages tend to be associated with tighter agreement.

The formal comparison confirms that these taxonomic patterns are meaningful overall, even if many individual pairwise contrasts remain weak after adjustment. In practical terms, Section 8 suggests that later interpretation should treat taxonomic uncertainty as real and structured, with wider caution around magnitude estimates for some groups than for others. It also reinforces the broader scenario conclusion that policy_3 remains the strongest package, while leaving uncertainty around the precise size of gains for the most uncertain taxa.

---

9 Habitat group-level expert agreement, uncertainty structure, and robustness

This subsection repeats the same diagnostics after aggregating outcomes within habitat groups. This shows whether expert agreement is tighter or looser when the elicitation is viewed through broad habitat structure rather than taxonomic structure.

practical note: these chunks assume df_AB and abs_policy_mag_signed already exist exactly as in Section 7.

9.1 Scenario-level habitat comparion

9.1.1 Dispersion in signed magnitude (agreement proxy) — Assumption B

## ---- habitat-level-setup ----------------------------------------------------
# Habitat group-level aggregation of absolute signed outcomes
# Mirrors the taxon-level chunk, but aggregates within habitat_group first.

habitat_AB <- df_AB %>%
  filter(str_detect(scenario, "^policy_") | scenario == "counterfactual") %>%
  group_by(assumption, id_gj, scenario, habitat_group) %>%
  summarise(
    mean_signed_habitat = mean(abs_policy_mag_signed, na.rm = TRUE),
    .groups = "drop"
  )

habitat_AB %>%
  summarise(
    n_rows = n(),
    n_assumptions = n_distinct(assumption),
    n_experts = n_distinct(id_gj),
    n_scenarios = n_distinct(scenario),
    n_habitat_groups = n_distinct(habitat_group)
  ) %>%
  kable(caption = "Coverage check: habitat group-level uncertainty dataset") %>%
  kable_styling(full_width = FALSE)

Coverage check: habitat group-level uncertainty dataset

n_rows	n_assumptions	n_experts	n_scenarios	n_habitat_groups
210	2	7	5	3

## ---- habitat-dispersion -----------------------------------------------------
# 1) Dispersion in signed outcomes across experts, at habitat group level

habitat_consensus_var_AB <- habitat_AB %>%
  group_by(assumption, habitat_group, scenario) %>%
  summarise(
    mean_signed = mean(mean_signed_habitat, na.rm = TRUE),
    sd_signed   = sd(mean_signed_habitat, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  )

habitat_consensus_summary_AB <- habitat_consensus_var_AB %>%
  group_by(assumption, scenario) %>%
  summarise(
    mean_sd   = mean(sd_signed, na.rm = TRUE),
    median_sd = median(sd_signed, na.rm = TRUE),
    .groups = "drop"
  )

habitat_consensus_summary_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Mean dispersion (SD) in signed outcomes by scenario at habitat group level (Assumptions A and B)") %>%
  kable_styling(full_width = FALSE)

Mean dispersion (SD) in signed outcomes by scenario at habitat group level (Assumptions A and B)

assumption	scenario	mean_sd	median_sd
A (rows = all survey_taxa)	counterfactual	10.27	10.05
A (rows = all survey_taxa)	policy_1	11.15	10.65
A (rows = all survey_taxa)	policy_2	9.02	8.16
A (rows = all survey_taxa)	policy_3	7.65	6.72
A (rows = all survey_taxa)	policy_4	8.85	7.92
B (collapsed = unique taxa:habitat_group)	counterfactual	10.35	10.14
B (collapsed = unique taxa:habitat_group)	policy_1	11.34	10.99
B (collapsed = unique taxa:habitat_group)	policy_2	9.18	8.46
B (collapsed = unique taxa:habitat_group)	policy_3	7.81	6.94
B (collapsed = unique taxa:habitat_group)	policy_4	8.99	8.14

## ---- habitat-dispersion-plot ------------------------------------------------
ggplot(
  habitat_consensus_var_AB,
  aes(x = scenario, y = sd_signed)
) +
  geom_boxplot() +
  facet_wrap(~ assumption) +
  labs(
    title = "Expert disagreement by scenario at habitat group level",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

At the habitat-group level, expert disagreement under Assumption B is moderate under the counterfactual and weakest policy package, but declines under the stronger interventions. Mean SD is 10.35 under the counterfactual, rises slightly to 11.34 under policy_1, then falls to 9.18 under policy_2, 7.81 under policy_3, and 8.99 under policy_4. Median SD shows the same pattern. This indicates that, once responses are averaged within the three broad habitat groups, stronger policy packages are associated with tighter agreement across the expert panel rather than greater divergence.

This is notable because it suggests that much of the panel disagreement is being smoothed out when ecological responses are viewed through broad habitat structure. In other words, experts may differ over particular taxa–habitat combinations, but they are more aligned about the overall direction and approximate scale of response at habitat-group level. The habitat-level analysis is based on three habitat groups, five scenarios and all seven experts under Assumption B, so it should be read as a broad structural view of agreement rather than a fine-grained ecological one.

9.1.2 Certainty structure — Assumption B

These analyses examine whether expert disagreement differs systematically among the specific habitat groups represented in the elicitation. Results are shown first by habitat group within each scenario, and then across all scenarios combined in order to identify any more persistent habitat-level patterns in certainty or disagreement. Only Assumption B is considered.

## ---- habitat-mean-variance --------------------------------------------------
# 2) Mean–variance relationship at habitat group level

ggplot(
  habitat_consensus_var_AB,
  aes(x = mean_signed, y = sd_signed)
) +
  geom_point(alpha = 0.7) +
  facet_grid(assumption ~ scenario) +
  labs(
    title = "Mean–variance relationship across habitat groups",
    x = "Mean signed % outcome",
    y = "Standard deviation across experts"
  ) +
  theme_minimal()

## ---- habitat-mean-variance-correlation --------------------------------------
habitat_mv_diag_AB <- habitat_consensus_var_AB %>%
  mutate(abs_mean_signed = abs(mean_signed)) %>%
  group_by(assumption, scenario) %>%
  summarise(
    cor_absmean_sd = cor(abs_mean_signed, sd_signed, use = "complete.obs"),
    .groups = "drop"
  )

habitat_mv_diag_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Correlation between absolute mean outcome and disagreement at habitat group level") %>%
  kable_styling(full_width = FALSE)

Correlation between absolute mean outcome and disagreement at habitat group level

assumption	scenario	cor_absmean_sd
A (rows = all survey_taxa)	counterfactual	-0.97
A (rows = all survey_taxa)	policy_1	-0.96
A (rows = all survey_taxa)	policy_2	-0.83
A (rows = all survey_taxa)	policy_3	1.00
A (rows = all survey_taxa)	policy_4	0.63
B (collapsed = unique taxa:habitat_group)	counterfactual	-0.96
B (collapsed = unique taxa:habitat_group)	policy_1	-0.96
B (collapsed = unique taxa:habitat_group)	policy_2	-0.14
B (collapsed = unique taxa:habitat_group)	policy_3	1.00
B (collapsed = unique taxa:habitat_group)	policy_4	0.84

The mean–variance diagnostic suggests that the structure of disagreement changes across scenarios. Under Assumption B, the correlation between absolute mean outcome and disagreement is strongly negative for the counterfactual (-0.96) and policy_1 (-0.96), close to zero for policy_2 (-0.14), and then strongly positive for policy_3 (+1.00) and policy_4 (+0.84). This implies that under the counterfactual and weakest policy package, habitat groups with larger mean responses are not the ones generating the greatest disagreement. By contrast, under the stronger interventions, the habitat groups with the largest projected responses are also the ones with the greatest dispersion among experts.

Substantively, this points to a shift from relatively uniform pessimism under the baseline towards magnitude uncertainty under stronger policy. Under policy_3 and policy_4, the panel appears broadly aligned that habitat groups improve, but less aligned on exactly how large those gains would be in the best-performing habitats. That interpretation is consistent with the wider habitat-group results, which show strong positive uplifts across all three habitat groups under the stronger scenarios, especially under policy_3.

9.1.3 Resampling stability — Assumption B

## ---- habitat-resampling -----------------------------------------------------
# 3) Resampling / panel-size sensitivity at habitat group level

set.seed(123)

habitat_resample_means <- function(k, assumption_in, scenario_in, n_iter = 1000) {
  dat <- habitat_AB %>%
    filter(assumption == assumption_in, scenario == scenario_in)

  map_dbl(seq_len(n_iter), ~{
    sampled_ids <- sample(unique(dat$id_gj), size = k, replace = FALSE)

    dat %>%
      filter(id_gj %in% sampled_ids) %>%
      group_by(habitat_group) %>%
      summarise(
        mean_signed_habitat = mean(mean_signed_habitat, na.rm = TRUE),
        .groups = "drop"
      ) %>%
      summarise(panel_mean = mean(mean_signed_habitat, na.rm = TRUE)) %>%
      pull(panel_mean)
  })
}

habitat_resampling_results_AB <- expand_grid(
  assumption = unique(habitat_AB$assumption),
  scenario   = c("counterfactual", "policy_1", "policy_2", "policy_3", "policy_4")
) %>%
  mutate(
    k2 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(2, a, s))),
    k3 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(3, a, s))),
    k4 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(4, a, s))),
    k5 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(5, a, s))),
    k6 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(6, a, s))),
    k7 = map2_dbl(assumption, scenario, \(a, s) sd(habitat_resample_means(7, a, s)))
  )

habitat_resampling_results_AB %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  kable(caption = "Resampling SD of mean signed outcome by panel size at habitat group level") %>%
  kable_styling(full_width = FALSE)

Resampling SD of mean signed outcome by panel size at habitat group level

assumption	scenario	k2	k3	k4	k5	k6	k7
A (rows = all survey_taxa)	counterfactual	6.04	4.29	3.26	2.39	1.56	0
A (rows = all survey_taxa)	policy_1	6.32	4.65	3.51	2.61	1.64	0
A (rows = all survey_taxa)	policy_2	5.00	3.69	2.83	1.99	1.37	0
A (rows = all survey_taxa)	policy_3	4.16	2.99	2.30	1.70	1.08	0
A (rows = all survey_taxa)	policy_4	4.87	3.58	2.69	1.93	1.26	0
B (collapsed = unique taxa:habitat_group)	counterfactual	6.01	4.44	3.29	2.35	1.58	0
B (collapsed = unique taxa:habitat_group)	policy_1	6.57	4.86	3.54	2.66	1.68	0
B (collapsed = unique taxa:habitat_group)	policy_2	5.26	3.80	2.87	2.07	1.33	0
B (collapsed = unique taxa:habitat_group)	policy_3	4.34	3.20	2.40	1.74	1.03	0
B (collapsed = unique taxa:habitat_group)	policy_4	4.98	3.64	2.70	1.99	1.27	0

## ---- habitat-resampling-plot ------------------------------------------------
habitat_resampling_long_AB <- habitat_resampling_results_AB %>%
  pivot_longer(
    cols = starts_with("k"),
    names_to = "panel_size",
    values_to = "sd_mean"
  ) %>%
  mutate(panel_size = as.numeric(str_remove(panel_size, "^k")))

ggplot(
  habitat_resampling_long_AB,
  aes(x = panel_size, y = sd_mean)
) +
  geom_line() +
  geom_point() +
  facet_grid(assumption ~ scenario) +
  labs(
    title = "Panel-size sensitivity of mean signed outcome at habitat group level",
    x = "Number of experts sampled",
    y = "SD of mean signed outcome"
  ) +
  theme_minimal()

The resampling analysis indicates that the habitat-group results under Assumption B are fairly robust to panel composition. For all scenarios, the SD of the resampled mean declines steadily as more experts are included, showing that panel-level conclusions stabilise in the expected way with increasing panel size. Among the policy scenarios, policy_3 is the most stable at larger panel sizes, with SD falling from 4.34 at k = 2 to 1.03 at k = 6. Policy_4 and policy_2 also stabilise clearly, while policy_1 remains the least stable of the policy packages.

This strengthens confidence in the habitat-level ranking of scenarios. Although there remains some uncertainty over exact magnitude, there is little evidence that the broad conclusions are being driven by one or two particular respondents. The fact that policy_3 combines the lowest dispersion at habitat-group level with the greatest resampling stability suggests that its stronger performance is not only a feature of the central estimate, but is also comparatively robust to moderate changes in panel composition.

9.1.4. Summary

Under Assumption B, the habitat-group analysis suggests that expert agreement becomes stronger once outcomes are viewed through broad habitat structure, particularly under the stronger policy scenarios. Disagreement is moderate under the counterfactual and policy_1, but falls under policy_2 and is lowest under policy_3, before increasing slightly again under policy_4. This indicates that habitat-level responses are interpreted more coherently by the panel than the more disaggregated taxa × habitat results might imply.

The mean–variance diagnostic adds an important qualification. Under the baseline and weakest policy package, more extreme habitat-group means are not associated with greater disagreement, whereas under policy_3 and policy_4 the relationship turns strongly positive. This suggests that uncertainty under stronger intervention is increasingly about the size of gains in the best-performing habitat groups, rather than about whether those gains occur at all.

Finally, the resampling results suggest that the habitat-level conclusions are robust to the composition of the seven-person panel. Resampled variability declines smoothly as panel size increases across all scenarios, and policy_3 emerges as the most stable of the policy options at larger panel sizes. Taken together, these results support a clear interpretation under Assumption B: stronger policies, especially policy_3, generate both better and more coherent habitat-level outcomes, with remaining uncertainty concentrated primarily around magnitude rather than direction.

9.2 Habitat-group uncertainty within each scenario

This section is based on Assumption B only throughout.

These additional analyses examine whether expert disagreement differs systematically among the specific habitat groups included in the elicitation. Results are shown first by habitat group within each scenario, and then across all scenarios combined in order to identify any more persistent habitat-level patterns in certainty or disagreement. Only Assumption B is considered.

Set-up

## ---- habitat-uncertainty-B-setup --------------------------------------------
# Restrict to Assumption B only for all additional habitat-group analyses

habitat_B <- habitat_AB %>%
  filter(assumption == "B (collapsed = unique taxa:habitat_group)")

habitat_B %>%
  summarise(
    n_rows = n(),
    n_experts = n_distinct(id_gj),
    n_scenarios = n_distinct(scenario),
    n_habitat_groups = n_distinct(habitat_group)
  ) %>%
  knitr::kable(
    caption = "Coverage check for additional habitat-group uncertainty analyses (Assumption B only)"
  )

Coverage check for additional habitat-group uncertainty analyses (Assumption B only)

n_rows	n_experts	n_scenarios	n_habitat_groups
105	7	5	3

9.2.A. Mean and SD by habitat group within each scenario

## ---- habitat-by-scenario-uncertainty-table ----------------------------------
# For each scenario, summarise expert disagreement separately for each habitat group

habitat_by_scenario_B <- habitat_B %>%
  group_by(scenario, habitat_group) %>%
  summarise(
    mean_signed = mean(mean_signed_habitat, na.rm = TRUE),
    sd_signed   = sd(mean_signed_habitat, na.rm = TRUE),
    min_signed  = min(mean_signed_habitat, na.rm = TRUE),
    max_signed  = max(mean_signed_habitat, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  arrange(scenario, desc(sd_signed))

habitat_by_scenario_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Expert disagreement by habitat group within each scenario (Assumption B)"
  )

Expert disagreement by habitat group within each scenario (Assumption B)

scenario	habitat_group	mean_signed	sd_signed	min_signed	max_signed	n
counterfactual	Woodland	-8.45	11.47	-30	2.50	7
counterfactual	Open habitats	-13.20	10.14	-30	-1.50	7
counterfactual	Freshwater/Wetland	-13.49	9.45	-30	-0.39	7
policy_1	Woodland	-4.17	13.16	-25	12.50	7
policy_1	Open habitats	-8.62	10.99	-25	1.83	7
policy_1	Freshwater/Wetland	-8.93	9.88	-25	0.83	7
policy_2	Woodland	2.21	11.25	-20	15.00	7
policy_2	Open habitats	-1.80	8.46	-20	3.67	7
policy_2	Freshwater/Wetland	-2.58	7.82	-20	2.17	7
policy_3	Woodland	6.99	9.84	-10	18.75	7
policy_3	Freshwater/Wetland	3.98	6.94	-10	11.78	7
policy_3	Open habitats	3.52	6.64	-10	10.58	7
policy_4	Woodland	1.99	10.75	-15	15.33	7
policy_4	Freshwater/Wetland	-0.51	8.14	-15	9.89	7
policy_4	Open habitats	-1.29	8.10	-15	8.58	7

## ---- habitat-by-scenario-uncertainty-plot -----------------------------------
ggplot(
  habitat_by_scenario_B,
  aes(x = reorder(habitat_group, sd_signed), y = sd_signed)
) +
  geom_col() +
  coord_flip() +
  facet_wrap(~ scenario, scales = "free_y") +
  labs(
    title = "Expert disagreement by habitat group within each scenario",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

#### Heatmap view for fast comparison across scenarios 
## ---- habitat-by-scenario-heatmap --------------------------------------------
ggplot(
  habitat_by_scenario_B,
  aes(x = scenario, y = habitat_group, fill = sd_signed)
) +
  geom_tile() +
  geom_text(aes(label = round(sd_signed, 1)), size = 3) +
  labs(
    title = "Heatmap of expert disagreement by habitat group and scenario",
    x = NULL,
    y = NULL,
    fill = "SD"
  ) +
  theme_minimal()

Rank habitat groups within each scenario by disagreement

## ---- habitat-by-scenario-ranking --------------------------------------------
habitat_by_scenario_ranks_B <- habitat_by_scenario_B %>%
  group_by(scenario) %>%
  mutate(
    disagreement_rank = min_rank(desc(sd_signed))
  ) %>%
  ungroup() %>%
  arrange(scenario, disagreement_rank)

habitat_by_scenario_ranks_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Ranking of habitat groups by expert disagreement within each scenario (Assumption B)"
  )

Ranking of habitat groups by expert disagreement within each scenario (Assumption B)

scenario	habitat_group	mean_signed	sd_signed	min_signed	max_signed	n	disagreement_rank
counterfactual	Woodland	-8.45	11.47	-30	2.50	7	1
counterfactual	Open habitats	-13.20	10.14	-30	-1.50	7	2
counterfactual	Freshwater/Wetland	-13.49	9.45	-30	-0.39	7	3
policy_1	Woodland	-4.17	13.16	-25	12.50	7	1
policy_1	Open habitats	-8.62	10.99	-25	1.83	7	2
policy_1	Freshwater/Wetland	-8.93	9.88	-25	0.83	7	3
policy_2	Woodland	2.21	11.25	-20	15.00	7	1
policy_2	Open habitats	-1.80	8.46	-20	3.67	7	2
policy_2	Freshwater/Wetland	-2.58	7.82	-20	2.17	7	3
policy_3	Woodland	6.99	9.84	-10	18.75	7	1
policy_3	Freshwater/Wetland	3.98	6.94	-10	11.78	7	2
policy_3	Open habitats	3.52	6.64	-10	10.58	7	3
policy_4	Woodland	1.99	10.75	-15	15.33	7	1
policy_4	Freshwater/Wetland	-0.51	8.14	-15	9.89	7	2
policy_4	Open habitats	-1.29	8.10	-15	8.58	7	3

Below section text is based on Assumption B only throughout.

At the habitat-group level, disagreement is not evenly distributed across the three habitat groupings. Woodland is the most uncertain habitat group under both the counterfactual (SD 11.47) and the weakest policy package (SD 13.16), while Freshwater/Wetland and Open habitats are somewhat lower and fairly similar to one another under those same scenarios. Under the stronger scenarios, disagreement declines for all three habitat groups, but the ranking remains broadly similar: Woodland continues to sit at the high-disagreement end, while Open habitats and Freshwater/Wetland tend to be somewhat lower.

This suggests that uncertainty is not simply a property of the policy scenario itself. Some of it is habitat-specific. In particular, Woodland appears to generate the widest spread of expert judgement across several scenarios, implying greater uncertainty about the size or consistency of projected woodland responses. By contrast, the two non-woodland groups appear to elicit somewhat tighter expert agreement, especially under the stronger packages.

A second notable pattern is that disagreement generally falls under the stronger packages, particularly under policy_3. For Freshwater/Wetland, SD falls from 10.82 under the counterfactual to 6.94 under policy_3; for Open habitats it falls from 10.14 to 6.64; and for Woodland from 11.47 to 9.84. This implies that the strongest package does not create greater habitat-level ambiguity. Rather, it appears to produce a more coherent expert signal across all three habitat groups, even though Woodland remains the least certain of the three.

---

9.3 Habitat-group uncertainty across all scenarios combined

This section is based on Assumption B only throughout.

This approach is the probably a stronger-signal version.

9.3.A Overall disagreement by habitat group across all scenarios

## ---- habitat-overall-uncertainty-table --------------------------------------
# Pool all scenarios together to test whether some habitat groups are
# consistently associated with more disagreement than others

habitat_overall_B <- habitat_B %>%
  group_by(habitat_group) %>%
  summarise(
    mean_signed   = mean(mean_signed_habitat, na.rm = TRUE),
    sd_signed     = sd(mean_signed_habitat, na.rm = TRUE),
    median_signed = median(mean_signed_habitat, na.rm = TRUE),
    iqr_signed    = IQR(mean_signed_habitat, na.rm = TRUE),
    min_signed    = min(mean_signed_habitat, na.rm = TRUE),
    max_signed    = max(mean_signed_habitat, na.rm = TRUE),
    n = n(),
    .groups = "drop"
  ) %>%
  arrange(desc(sd_signed))

habitat_overall_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Overall expert disagreement by habitat group across all scenarios combined (Assumption B)"
  )

Overall expert disagreement by habitat group across all scenarios combined (Assumption B)

habitat_group	mean_signed	sd_signed	median_signed	iqr_signed	min_signed	max_signed	n
Woodland	-0.29	11.99	1.67	15.42	-30	18.75	35
Open habitats	-4.28	10.36	0.00	12.79	-30	10.58	35
Freshwater/Wetland	-4.31	10.17	-1.67	11.71	-30	11.78	35

## ---- habitat-overall-uncertainty-plot ---------------------------------------
ggplot(
  habitat_overall_B,
  aes(x = reorder(habitat_group, sd_signed), y = sd_signed)
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Overall expert disagreement by habitat group across all scenarios combined",
    x = NULL,
    y = "SD across experts"
  ) +
  theme_minimal()

## ---- habitat-overall-distribution-plot --------------------------------------
ggplot(
  habitat_B,
  aes(x = habitat_group, y = mean_signed_habitat)
) +
  geom_boxplot() +
  coord_flip() +
  labs(
    title = "Distribution of expert habitat-group outcomes across all scenarios combined",
    x = NULL,
    y = "Mean signed % outcome"
  ) +
  theme_minimal()

9.3.B Do habitat groups differ in disagreement? (more formal comparison)

This tests whether the spread across experts differs systematically among habitat groups.

Seup: per habitat group, we calculate within-scenario SD first

## ---- habitat-sd-dataset-for-testing -----------------------------------------
# Create a compact dataset of one disagreement value per habitat group per scenario

habitat_sd_for_test_B <- habitat_B %>%
  group_by(habitat_group, scenario) %>%
  summarise(
    sd_signed = sd(mean_signed_habitat, na.rm = TRUE),
    .groups = "drop"
  )

habitat_sd_for_test_B %>%
  mutate(sd_signed = round(sd_signed, 2)) %>%
  knitr::kable(
    caption = "Scenario-specific disagreement values used for habitat-group comparison tests (Assumption B)"
  )

Scenario-specific disagreement values used for habitat-group comparison tests (Assumption B)

habitat_group	scenario	sd_signed
Freshwater/Wetland	counterfactual	9.45
Freshwater/Wetland	policy_1	9.88
Freshwater/Wetland	policy_2	7.82
Freshwater/Wetland	policy_3	6.94
Freshwater/Wetland	policy_4	8.14
Open habitats	counterfactual	10.14
Open habitats	policy_1	10.99
Open habitats	policy_2	8.46
Open habitats	policy_3	6.64
Open habitats	policy_4	8.10
Woodland	counterfactual	11.47
Woodland	policy_1	13.16
Woodland	policy_2	11.25
Woodland	policy_3	9.84
Woodland	policy_4	10.75

9.3.B.1 Kruskal–Wallis test across habitat groups

## ---- habitat-kruskal-test ---------------------------------------------------
habitat_kw_B <- kruskal.test(sd_signed ~ habitat_group, data = habitat_sd_for_test_B)

habitat_kw_B

## 
##  Kruskal-Wallis rank sum test
## 
## data:  sd_signed by habitat_group
## Kruskal-Wallis chi-squared = 6.86, df = 2, p-value = 0.03239

9.3.B.2 Pairwise Wilcoxon tests between habitat groups

## ---- habitat-pairwise-tests -------------------------------------------------
pairwise.wilcox.test(
  x = habitat_sd_for_test_B$sd_signed,
  g = habitat_sd_for_test_B$habitat_group,
  p.adjust.method = "BH",
  exact = FALSE
)

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  habitat_sd_for_test_B$sd_signed and habitat_sd_for_test_B$habitat_group 
## 
##               Freshwater/Wetland Open habitats
## Open habitats 0.676              -            
## Woodland      0.065              0.090        
## 
## P value adjustment method: BH

9.3.B.3 simpler “certainty score” by habitat group

This converts disagreement into a relative certainty-like metric for quick reading.

## ---- habitat-relative-certainty-score ---------------------------------------
habitat_certainty_B <- habitat_overall_B %>%
  mutate(
    certainty_score = 1 - (sd_signed - min(sd_signed)) / (max(sd_signed) - min(sd_signed))
  ) %>%
  arrange(desc(certainty_score))

habitat_certainty_B %>%
  mutate(across(where(is.numeric), \(x) round(x, 2))) %>%
  knitr::kable(
    caption = "Relative certainty score by habitat group across all scenarios combined (Assumption B only)"
  )

Relative certainty score by habitat group across all scenarios combined (Assumption B only)

habitat_group	mean_signed	sd_signed	median_signed	iqr_signed	min_signed	max_signed	n	certainty_score
Freshwater/Wetland	-4.31	10.17	-1.67	11.71	-30	11.78	35	1.00
Open habitats	-4.28	10.36	0.00	12.79	-30	10.58	35	0.89
Woodland	-0.29	11.99	1.67	15.42	-30	18.75	35	0.00

## ---- habitat-certainty-score-plot -------------------------------------------
ggplot(
  habitat_certainty_B,
  aes(x = reorder(habitat_group, certainty_score), y = certainty_score)
) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Relative certainty by habitat group across all scenarios combined",
    x = NULL,
    y = "Relative certainty score (0-1)"
  ) +
  theme_minimal()

When all scenarios are considered together, the habitat-level pattern becomes clearer. Woodland shows the greatest overall disagreement across experts, while Freshwater/Wetland shows the lowest and therefore the highest relative certainty. Open habitats sit between the two. The relative certainty score puts Freshwater/Wetland at 1.00, with Open habitats only slightly lower, and Woodland clearly at the bottom of the scale.

This pooled result strengthens the interpretation from the scenario-by-scenario outputs. Woodland is not only somewhat more uncertain in individual scenarios, but appears to be persistently more uncertain across the elicitation as a whole. Freshwater/Wetland, by contrast, appears to be the habitat group for which the panel is most consistently aligned. That does not mean experts are optimistic about Freshwater/Wetland outcomes — its pooled mean remains negative — but rather that they are relatively consistent in that assessment.

The formal comparison supports treating this as a real habitat-level pattern. The Kruskal–Wallis test indicates that disagreement differs significantly across habitat groups overall (chi-squared = 6.86, df = 2, p = 0.032). At the same time, the pairwise Wilcoxon tests are weaker after adjustment, with Woodland tending towards higher disagreement relative to the other two groups but without strong pairwise separation at conventional thresholds. This suggests a broad habitat-level gradient in certainty, rather than sharply distinct and fully separable categories.

(9.2-3 summary)

Below section text is based on Assumption B only throughout.

Taken together, these additional habitat-level analyses show that expert disagreement is structured not only by scenario, but also by habitat group. Woodland is consistently the most uncertain habitat grouping, while Freshwater/Wetland is the most certain and Open habitats generally intermediate. This pattern is visible both within individual scenarios and when all scenarios are pooled.

The scenario-specific results also suggest that stronger policy packages tend to reduce disagreement across all habitat groups, especially under policy_3. This indicates that the strongest package is not merely associated with larger projected effects, but also with a clearer shared signal among experts at habitat-group level. Remaining uncertainty is therefore better understood as uncertainty about the precise magnitude of gains, especially in Woodland, rather than as disagreement about the basic direction of response.

The formal tests reinforce this interpretation. There is evidence of a statistically meaningful habitat-group pattern in disagreement overall, but with limited power to separate individual pairs cleanly. In practical terms, this supports carrying forward habitat-structured uncertainty into later interpretation: woodland estimates should be treated with somewhat greater caution, whereas freshwater/wetland and open-habitat estimates appear to support somewhat tighter inference.

---

9.4 Summary

Section 9 shows that habitat-group uncertainty is also structured, but across a much smaller set of groups. At the broad habitat-group level, disagreement under Assumption B is moderate under the counterfactual and policy_1, but declines under the stronger scenarios, especially under policy_3. The more detailed habitat-group analyses then show that Woodland is consistently the most uncertain habitat grouping, while Freshwater/Wetland is the most certain and Open habitats is generally intermediate.

This means that habitat-level uncertainty is not just a reflection of scenario strength. Some habitat contexts appear intrinsically harder for the panel to judge. Even so, the stronger scenarios produce a clearer shared signal than the weaker ones, and the formal tests suggest that the habitat-group pattern is meaningful overall. In practical terms, Section 9 supports carrying forward habitat-structured uncertainty, with somewhat greater caution around woodland responses and somewhat greater confidence in the relative consistency of freshwater/wetland and open-habitat assessments.

---

---

10. Integrated synthesis and implications for the extinction-risk target (ERLI)

Important: Below text refers to Assumption B only throughout. See Appendix A for an older version, which includes Assumption A reults.

10.1 Integrated interpretation of Sections 7–9 under Assumption B

The combined message from Sections 7–9 under Assumption B is clear. The expert panel gives a strongly pessimistic counterfactual to 2042, and all four policy packages are judged to improve on that baseline, with a stable ranking in which policy_3 performs best, policy_1 performs worst, and policy_2 and policy_4 occupy an intermediate position. However, a positive policy uplift is not the same thing as recovery: even under the strongest package, mean absolute outcomes remain below zero, indicating that policy improvement does not fully offset the expected background decline.

The uncertainty analysis sharpens rather than weakens that conclusion. At the overall level, disagreement is real but coherent. It is not random noise, nor is it concentrated in a way that undermines the central ranking of scenarios. Resampling shows that headline means stabilise as panel size increases, and supports confidence in the directional result and in the relative ranking of policy options, while leaving the precise distance from decline to stability under the strongest package more uncertain. In other words, the panel is much more certain that policy helps than it is about exactly how close the best package gets to zero net change.

Sections 8 and 9 then show that this uncertainty is ecologically structured. Some taxa are consistently harder to judge than others, with vascular plants, beetles and often Diptera generating wider spreads of opinion, whereas lichens, bryophytes and freshwater invertebrate assemblages tend to generate tighter agreement. Similarly, some habitat groups are more uncertain than others, with Woodland emerging as the least certain habitat grouping and Freshwater/Wetland the most certain. This means the panel is not equally confident across the ecological space represented in the elicitation. The strongest inferences can therefore be made about the broad scenario ordering and the existence of residual decline, while somewhat greater caution is needed when interpreting the magnitude of gains in the most uncertain taxa and habitats.

Taken together, Sections 7–9 imply three things for interpretation. First, the elicitation provides a robust signal that baseline pressures continue to drive decline by 2042. Second, it provides an equally robust signal that stronger policy packages improve outcomes, with policy_3 consistently emerging as the best-performing option. Third, it shows that the remaining uncertainty is concentrated mainly around how much improvement is achieved, especially where responses are closest to the stability threshold and in the ecological groupings where disagreement is intrinsically higher. This is exactly the pattern one would expect from a credible elicitation: strong agreement on direction and ordering, but less certainty about precise magnitude at the optimistic edge.

10.2 What this means for the extinction-risk target (ERLI)

Under Assumption B, the central implication for the extinction-risk target is that the panel does not expect the system to recover to a positive mean state by 2042 under the currently defined packages, even though all packages improve on the counterfactual. The strongest scenario is judged to bring the system closest to stability, but not to eliminate decline in the mean. This matters because ERLI-relevant interpretation should distinguish clearly between beneficial intervention and sufficient intervention. The elicitation supports the former strongly, but only weakly supports the latter.

The results also suggest that any future modelling or policy interpretation should pay particular attention to the persistence of residual pockets of decline. Under the strongest package, many units improve and many are stable or positive, but a material minority remain negative. Those negative tails are what keep the mean below zero. For ERLI interpretation, that means progress could be real and substantial while still falling short of the level needed to reverse the overall extinction-risk trajectory.

At a more ecological level, the habitat results suggest that recovery constraints are not evenly distributed. Freshwater/Wetland and Open habitats remain especially important because they continue to show negative mean outcomes even under the strongest package, while Woodland comes closest to neutrality. The taxonomic results suggest an analogous pattern: some taxa contribute more heavily to uncertainty and may also be more likely to contain those remaining pockets of decline. This does not overturn the overall scenario ranking, but it does imply that achieving the ERLI target may require not only a strong overall package, but also sharper attention to the taxa and habitats where decline is hardest to arrest and expert uncertainty remains greatest.

10.3 Recommended interpretation of the elicitation outputs

Under Assumption B, the elicitation outputs are best used in three complementary ways.

• First, as evidence that the counterfactual is strongly negative and that all policy packages provide positive uplift relative to that baseline.
• Second, as evidence that the ordering of packages is stable, with policy_3 consistently best and policy_1 consistently weakest.
• Third, as a basis for propagating structured uncertainty into any subsequent modelling of ERLI trajectories, especially where outcomes are close to zero and where disagreement is concentrated in particular taxa or habitats.

This means the elicitation should not be summarised simply as “policy_3 solves the problem” or, conversely, as “uncertainty is too large to say anything useful.” Neither would be justified. A more accurate reading is that the panel gives a clear directional and comparative signal, but that the magnitude of remaining decline versus near-stability remains uncertain enough to warrant explicit uncertainty bounds in any downstream ERLI modelling. That is particularly important because the strongest package is the one that is both most promising and most policy-relevant, precisely because it is the package closest to the threshold between continued decline and stabilisation.

10.4 Summary

Under Assumption B, the expert elicitation indicates that the extinction-risk baseline to 2042 is strongly negative, that all policy packages improve on that baseline, and that policy_3 is the most effective package. However, the elicitation does not support a confident claim that even the strongest package is sufficient to produce positive mean outcomes by 2042. Instead, it supports a more cautious conclusion: the strongest package likely delivers substantial improvement and may bring some ecological groupings close to stability, but residual decline remains likely in the mean and is concentrated in particular taxa and habitats.

The uncertainty diagnostics add confidence to this overall interpretation. They suggest that the scenario ranking is robust, that the panel is broadly aligned on direction of change, and that the main remaining uncertainty concerns how large the gains are, not whether they occur. That makes the elicitation a useful input to ERLI trajectory modelling, provided that downstream work carries forward ecological heterogeneity and structured uncertainty rather than collapsing everything into a single over-precise point estimate.

Appendices

A. Old Version 10. Integrated synthesis and implications for the extinction-risk target (ERLI)

Keeping this summary in here as I was rather happy with it and worry I might have forgotten something in the updated version.

• Counterfactual baseline is strongly negative: mean change is around −12% by 2042, with a high share of units negative.
• All policy scenarios are beneficial relative to baseline: mean deltas are positive for every scenario, with a consistent ranking (policy_3 > policy_2 ≈ policy_4 > policy_1).
• But uplifts are not sufficient in the mean: after combining CF + delta, the mean absolute outcome remains negative under every scenario; policy_3 is closest to stability but still below 0 in the mean.
• Many units improve, yet residual “pockets of decline” remain: under policy_3, most units are stable/positive, but a material minority remain negative; these tails keep the mean below zero.
• Habitat constraints are clear: Freshwater/Wetland and Open habitats remain negative even under the strongest scenario; Woodland is closest to neutrality.
• Uncertainty is structured and concentrated where leverage is highest: disagreement is largest under policy_3; resampling supports stable scenario ranking but uncertain distance-to-zero under ambitious intervention.

A.1 What was elicited, and what it is (and is not) telling us

This expert elicitation separates two components of expected change by 2042:

• Q1 (counterfactual): what experts expect under baseline pressures without the additional policy package.
• Q2 (policy delta): the incremental change expected relative to that counterfactual.

This structure is valuable because it prevents us from conflating “policies help” with “policies deliver recovery.” Q2 can be strongly positive while absolute outcomes remain negative if baseline decline is large.

All synthesis statements below are robust to the two aggregation assumptions used in the analysis:

• A: each survey_taxa row contributes (survey-design faithful, but can overweight taxa split into multiple rows within a habitat group),
• B: within each expert × scenario, multiple survey rows representing the same taxa within the same habitat group are collapsed to one unit (prevents row multiplicity driving habitat-group results).

Only a small subset of combinations are affected by the A→B collapse step (freshwater assemblages in Freshwater/Wetland; lichens in Open habitats), so A vs B is best seen as a sensitivity check for habitat-group interpretation rather than a competing “true” baseline.

Selective evidence to keep visible in this section: - Table: Counterfactual baseline (Q1 signed %): Assumptions A and B - Table + figure: Counterfactual baseline by habitat_group (Assumptions A and B)

A.2 Baseline pressure: the counterfactual is strongly negative

Across both assumptions, the counterfactual baseline mean is around −12%, with a high proportion of units negative.

Habitat-group baselines show that:

• Freshwater/Wetland and Open habitats have the most negative baseline means (large and widespread decline expectations),
• Woodland is less negative, but still predominantly negative.

This is a crucial policy framing point: the baseline “pull” is substantial, so “reasonable” policy uplifts may still be insufficient to deliver stability unless they are large, broad, and fast-acting.

Selective evidence to keep visible: - Table + figure: Counterfactual baseline by habitat_group (Assumptions A and B)

A.3 Incremental uplifts: policies help, and scenarios are clearly ranked

All policy scenarios show positive mean deltas relative to the counterfactual, with a consistent ranking:

• policy_3 is strongest (mean delta around +10.6% in A),
• policy_2 and policy_4 are intermediate (~+7–8%),
• policy_1 is weakest (~+2–3%).

The high proportion of positive deltas, especially under policy_3, indicates broad agreement that interventions improve outcomes relative to baseline. However, because these are deltas, this result should be communicated as “policies reduce decline relative to BAU,” not as “policies deliver recovery.”

Selective evidence to keep visible: - Table: Incremental policy uplifts (Q2 deltas): overall by scenario (Assumptions A and B)
- One figure: the scenario-level uplift bar chart (faceting by assumption if shown)

A.4 Sufficiency: CF + delta implies continued decline on average (even under the strongest package)

The key sufficiency test is the absolute outcome:

Absolute outcome under policy = counterfactual + delta

This is the quantity that matters for target framing (stability vs recovery). Under Assumption A, the mean absolute outcomes remain negative for every scenario:

• policy_1 ≈ −9.67%
• policy_2 ≈ −4.55%
• policy_4 ≈ −4.99%
• policy_3 ≈ −1.73% (closest to stability)

This means the policy packages, as elicited, are best characterised as decline-reduction pathways, not system-wide recovery pathways, at least in the mean.

At the same time, the distributions are heterogeneous: under policy_3, most units are stable/positive, yet a material minority remain negative, and these residual declines keep the mean below zero. This is exactly why it is helpful to report both the mean and the recovered/stabilised fractions.

Selective evidence to keep visible: - Table: Absolute signed % outcome under each policy (Assumptions A and B) - Figure: Mean signed % outcome: counterfactual vs policy

A.5 “Most improve” is not the same as “target met”: residual declines remain binding

The recovered/stabilised/decline-remaining table provides the most accessible translation:

• Under policy_3, around 79% of units are stabilised (≥0) and around 67–68% are positive (>0), while roughly 21% remain negative.

This reconciles the apparent tension between “many outcomes are positive” and “mean outcome is negative.” The minority of negative outcomes are often sufficiently negative to dominate the mean. From a policy perspective, this implies that:

• the strongest package generates widespread improvement,
• but there remain persistent and/or severe problem components that would need targeted strengthening if system-wide neutrality is required.

Selective evidence to keep visible: - Table: recovered/stabilised/decline remaining under A and B

A.6 Habitat-group implications: where the remaining gap sits

Habitat-group absolute outcomes indicate that:

• Freshwater/Wetland remains negative under all scenarios, including policy_3, even though the proportion positive rises strongly.
• Open habitats also remain negative; policy_3 improves outcomes materially but does not eliminate the negative mean.
• Woodland is closest to neutrality and can become slightly positive under the strongest package.

Because these habitat-group results are reported under both A and B, one can interpret them as a robust signal rather than a survey-row artefact.

Selective evidence to keep visible: - Table/scroll box: Absolute outcomes under policy by habitat_group (Assumptions A and B)
- One figure: habitat-group mean outcome bars by scenario (faceted by habitat_group and/or assumption)

A.7 Uncertainty and robustness: what we can say confidently, and what we should hedge

Uncertainty diagnostics show:

• Dispersion (expert disagreement) increases under more ambitious policy, with the largest SD under policy_3.
• The mean–variance relationship is consistent with structured epistemic uncertainty: experts diverge most on high-leverage outcomes.
• Resampling supports the interpretation that scenario ordering is not driven by a single expert, but that “how close to zero” policy_3 lands should be treated probabilistically rather than as a point estimate.

This supports a credible communication stance:

• High confidence: policies improve outcomes, and policy_3 is strongest.
• Medium confidence: policy_3 is close to stability in the mean.
• Lower confidence: the precise distance to zero and whether the system crosses into net-positive territory once uncertainties are propagated.

Selective evidence to keep visible: - Table: Mean dispersion (SD) in signed magnitude by scenario (Assumptions A and B) - One figure: mean–variance scatter (faceted by scenario), OR resampling stability plot (pick one, not both)

A.8 Implications for the extinction-risk target (ERLI framing)

The statutory extinction-risk target implies no net increase in extinction risk by 2042, and ideally improvement. We have not translated percent distribution-change expectations into explicit Red List category transitions (and therefore not directly into ERLI change), because that would require an additional mapping model and strong assumptions about how these percent changes relate to category thresholds.

However, the CF + delta results provide a strong sufficiency diagnostic:

• The counterfactual is strongly negative (continued deterioration pressure).
• Policy scenarios reduce that deterioration, but mean outcomes remain negative under all packages, including the strongest.
• Habitat-group patterns indicate persistent negative components (especially freshwater and open habitats) even under policy_3.

Therefore, on the elicited evidence alone, it is difficult to argue that any of the tested policy packages is likely to deliver system-wide stability/recovery by 2042 in a way that would reliably support ERLI neutrality or improvement. The most defensible framing is:

• policy_3 is closest and may deliver stabilisation/recovery for many components,
• but residual decline remains, implying a continuing risk of upward extinction-risk pressure unless additional measures close the remaining gap.