Outlier Exploration

Author

Jeffrey Girard

Published

February 13, 2026

1 Setup

1.1 Packages

library(MASS)
library(pracma)
library(tidyverse)
library(patchwork)
library(gt)
library(furrr)
plan(multisession, workers = availableCores() - 1)

1.2 Data

val_df <- read_rds("valence_tidy.rds")
hol_df <- read_rds("holistic_tidy.rds")

2 Synchrony (Dynamic ISC)

To detect participants who may be disengaged or watching “a different movie” than the rest of the sample, we calculate Inter-Subject Correlation (ISC). For each participant, we calculate the correlation between their dynamic rating time-series and the average time-series of all other participants (leave-one-out average).

  • Metric: Mean Pearson correlation (Fisher Z-transformed).
  • Variable Name: D_ISC (Dynamic ISC).
  • Threshold: \(r < 0.30\) (Flag for review).
# Calculate Leave-One-Out correlations for a single video's data
calc_video_isc <- function(video_data) {
  
  # 1. Pivot to wide format
  wide_mat <- video_data |>
    select(Rater, Second, Rating) |>
    pivot_wider(names_from = Rater, values_from = Rating, values_fill = NA) |>
    select(-Second) |>
    as.matrix()
  
  n_raters <- ncol(wide_mat)
  rater_ids <- colnames(wide_mat)
  correlations <- numeric(n_raters)
  
  for (i in 1:n_raters) {
    # 2. Target rater vector
    target_vec <- wide_mat[, i]
    
    # 3. "Gold Standard" (mean of everyone else)
    others_mean <- rowMeans(wide_mat[, -i], na.rm = TRUE)
    
    # 4. Safety Checks
    # Check if we have enough overlapping data points (>= 3) to correlate
    valid_overlap <- sum(!is.na(target_vec) & !is.na(others_mean))
    
    # Check for flat-lines (SD=0) which break correlation
    # We only check SD on the *overlapping* segments
    is_flat <- function(x) {
       v <- na.omit(x)
       if(length(v) < 2) return(TRUE)
       return(sd(v) == 0)
    }
    
    if (valid_overlap < 3 || is_flat(target_vec) || is_flat(others_mean)) {
      correlations[i] <- NA
    } else {
      # 5. Pearson correlation with pairwise deletion
      cor_val <- suppressWarnings(
        cor(target_vec, others_mean, use = "pairwise.complete.obs")
      )
      correlations[i] <- cor_val
    }
  }
  
  tibble(Rater = as.integer(rater_ids), r_val = correlations)
}
# Process ISC: Nest by video, map function, then average per Rater
rater_isc_summary <- 
  val_df |>
  nest(.by = Abbrev) |>
  mutate(isc_data = map(data, calc_video_isc)) |>
  select(Abbrev, isc_data) |>
  unnest(isc_data) |>
  mutate(z_val = 0.5 * log((1 + r_val) / (1 - r_val))) |> # Fisher Z
  summarize(
    D_ISC_Z = mean(z_val, na.rm = TRUE),
    .by = Rater
  ) |>
  mutate(
    # Transform back for interpretable thresholding
    D_ISC = (exp(2 * D_ISC_Z) - 1) / (exp(2 * D_ISC_Z) + 1)
  )

The histogram below visualizes the distribution of synchrony. Participants flagged in red fall below the 0.30 threshold.

bin_w <- 0.05

outliers <- rater_isc_summary |> 
  filter(D_ISC < 0.3) |> 
  mutate(
    bin_center = floor(D_ISC / bin_w) * bin_w + (bin_w / 2)
  )

ggplot(rater_isc_summary, aes(x = D_ISC)) +
  geom_histogram(
    color = "white",
    fill = "steelblue",
    breaks = seq(0, 1, by = bin_w)
  ) +
  geom_vline(xintercept = 0.3, linetype = "dashed", color = "red") +
  geom_text(
    data = outliers, 
    aes(x = bin_center, y = 1, label = Rater),
    vjust = -0.5,
    color = "red",
    fontface = "bold"
  ) +
  scale_x_continuous(
    limits = c(0, 1), 
    breaks = seq(0, 1, by = bin_w*2)
  ) +
  scale_y_continuous(limits = c(0, 20), expand = c(0, 0)) +
  labs(
    x = "Mean Inter-Subject Correlation (Dynamic)", 
    y = "Count (Participants)"
  ) +
  theme_classic()

3 Complexity (Dynamic ENT)

To detect participants who are either “flat-lining” (too smooth) or “randomly clicking” (too jagged), we calculate the Sample Entropy of their ratings.

We normalize this metric using robust Z-scores (Median/MAD) on the Log-Transformed entropy values to account for the floor effect at zero.

  • Metric: Sample Entropy (Log-Transformed).
  • Variable Name: D_ENT (Dynamic Entropy).
  • Threshold: \(Z < -2.5\) (Flat-liner) or \(Z > 2.5\) (Random Clicker).
calc_entropy_safe <- function(rating_vec) {
  clean_vec <- as.numeric(as.vector(na.omit(rating_vec)))
  if (length(clean_vec) < 3) {
    return(NA_real_)
  }
  if (sd(clean_vec) == 0) {
    return(0)
  }
  tryCatch(
    {
      val <- sample_entropy(clean_vec, edim = 2, r = 0.2 * sd(clean_vec))
      if (is.infinite(val) || is.nan(val)) {
        return(NA_real_)
      }
      return(as.numeric(val))
    },
    error = function(e) {
      warning(paste("Entropy Failed:", e$message))
      return(NA_real_)
    }
  )
}
# Process Entropy (Complexity)
rater_entropy_summary <- val_df |>
  nest(.by = c(Rater, Abbrev)) |>
  mutate(
    Samp_Entropy = future_map_dbl(
      data, 
      ~ calc_entropy_safe(.x$Rating),
      .options = furrr_options(seed = TRUE)
    )
  ) |>
  select(-data) |>
  summarize(
    Mean_Entropy = mean(Samp_Entropy, na.rm = TRUE),
    .by = Rater
  ) |>
  mutate(
    D_ENT = log(Mean_Entropy + 1e-6),
    D_ENT_Z = (D_ENT - median(D_ENT, na.rm = TRUE)) / mad(D_ENT, na.rm = TRUE)
  )

The histogram below shows the distribution of entropy. Outliers on the far left represent unusually static ratings, while outliers on the far right represent unusually chaotic ratings.

bin_w <- 0.25

outliers <- 
  rater_entropy_summary |> 
  filter(abs(D_ENT_Z) > 2.5) |> 
  mutate(bin_center = floor(D_ENT_Z / bin_w) * bin_w + (bin_w / 2)) |> 
  mutate(
    .by = bin_center,
    stack_rank = row_number(),
    y_stacked = 1 + max(stack_rank) + (stack_rank - 1)
  )

ggplot(rater_entropy_summary, aes(x = D_ENT_Z)) +
  geom_histogram(
    color = "white",
    fill = "steelblue",
    breaks = seq(-5, 5, by = bin_w)
  ) +
  geom_vline(xintercept = c(-2.5, 2.5), linetype = "dashed", color = "red") +
  geom_text(
    data = outliers, 
    aes(x = bin_center, y = y_stacked, label = Rater),
    size = 3,
    color = "red",
    fontface = "bold"
  ) +
  labs(
    x = "Standardized Log Entropy (Dynamic)", 
    y = "Count (Participants)"
  ) +
  theme_classic()

4 Intensity (Dynamic VAL)

To detect participants who are rating the videos with unusually positive or negative bias on average across the time-series, we calculate the Mean Dynamic Valence. We normalize this using Robust Z-scores.

  • Metric: Robust Z-Score (Median/MAD) of the mean dynamic rating.
  • Variable Name: D_VAL (Dynamic Valence).
  • Threshold: \(Z > 2.5\) (Extreme Positive Bias) or \(Z < -2.5\) (Extreme Negative Bias).
# Process Dynamic Valence
rater_val_summary <- val_df |>
  summarize(
    D_VAL = mean(Rating, na.rm = TRUE), 
    .by = Rater
  ) |>
  mutate(
    # Robust Z-score
    D_VAL_Z = (D_VAL - median(D_VAL, na.rm = TRUE)) / mad(D_VAL, na.rm = TRUE)
  )

The histogram below shows the distribution of entropy. Outliers on the far left represent unusually static ratings, while outliers on the far right represent unusually chaotic ratings.

bin_w <- 0.5

outliers <- 
  rater_val_summary |> 
  filter(abs(D_VAL_Z) > 2.5) |> 
  mutate(bin_center = floor(D_VAL_Z / bin_w) * bin_w + (bin_w / 2)) |> 
  mutate(
    .by = bin_center,
    stack_rank = row_number(),
    y_stacked = 1 + max(stack_rank) + (stack_rank - 1) * 1.25
  )

ggplot(rater_val_summary, aes(x = D_VAL_Z)) +
  geom_histogram(
    color = "white",
    fill = "steelblue",
    breaks = seq(-5, 5, by = bin_w)
  ) +
  geom_vline(xintercept = c(-2.5, 2.5), linetype = "dashed", color = "red") +
  geom_text(
    data = outliers, 
    aes(x = bin_center, y = y_stacked, label = Rater),
    size = 3,
    color = "red",
    fontface = "bold"
  ) +
  labs(
    x = "Standardized Valence Mean (Dynamic)", 
    y = "Count (Participants)"
  ) +
  theme_classic()

5 Intensity (Holistic POS/NEG)

To detect participants who are consistently rating unnaturally high or low on the holistic scales, we calculate Robust Z-Scores for both Positive and Negative affect separately.

  • Metric: Robust Z-Score (Median/MAD).
  • Variable Names: H_POS, H_NEG.
  • Threshold: \(Z > 2.5\) (Extreme Intensity) or \(Z < -2.5\) (Unusually Low).
# Aggregate holistic data
rater_hol_summary <- 
  hol_df |>
  summarize(
    Mean_Score = mean(Rating, na.rm = TRUE), 
    .by = c(Rater, Scale)
  ) |>
  pivot_wider(names_from = Scale, values_from = Mean_Score) |>
  rename(H_POS = Positive, H_NEG = Negative) |> 
  mutate(
    # Calculate Robust Z-scores for each scale separately
    H_POS_Z = (H_POS - median(H_POS, na.rm = TRUE)) / mad(H_POS, na.rm = TRUE),
    
    H_NEG_Z = (H_NEG - median(H_NEG, na.rm = TRUE)) / mad(H_NEG, na.rm = TRUE)
  )

The plots below show the distribution of Positive and Negative intensity.

# Function to plot univariate Z-scores
plot_intensity <- function(data, column, label) {
  bin_w <- 0.5
  col_sym <- rlang::ensym(column)
  
  break_seq <- seq(-5, 5, by = bin_w)
  
  outliers <- data |> 
    filter(abs(!!col_sym) > 2.5) |> 
    mutate(bin_center = floor(!!col_sym / bin_w) * bin_w + (bin_w / 2)) |> 
    arrange(bin_center, Rater) |> 
    mutate(
      stack_rank = row_number(),
      y_stacked = 1 + max(stack_rank) + (stack_rank - 1) * 2,
      .by = bin_center
    )
    
  ggplot(data, aes(x = !!col_sym)) +
    # Use breaks instead of binwidth
    geom_histogram(color = "white", fill = "steelblue", breaks = break_seq) +
    geom_vline(xintercept = c(-2.5, 2.5), linetype = "dashed", color = "red") +
    geom_text(
      data = outliers, 
      aes(x = bin_center, y = y_stacked, label = Rater),
      size = 3.25, color = "red", fontface = "bold"
    ) +
    labs(x = label, y = "Count") +
    theme_classic()
}

p1 <- plot_intensity(rater_hol_summary, H_POS_Z, "Standardized Positive Affect (Holistic)")
p2 <- plot_intensity(rater_hol_summary, H_NEG_Z, "Standardized Negative Affect (Holistic)")
p1 / p2

6 Final Integration

We combine all metrics into a single quality control table using standardized variable names:

  • D_ISC_Z: Dynamic Inter-Subject Correlation (Fisher Z)
  • D_ENT_Z: Dynamic Entropy (Log-transformed Z)
  • D_VAL_Z: Dynamic Valence (Mean Z)
  • H_POS_Z: Holistic Positive (Mean Z)
  • H_NEG_Z: Holistic Negative (Mean Z)

6.1 Correlation Check

# 1. Join and Standardize
final_qc_prep <- 
  rater_isc_summary |>
  left_join(rater_entropy_summary, by = "Rater") |>
  left_join(rater_val_summary, by = "Rater") |>
  left_join(rater_hol_summary, by = "Rater") |>
  transmute(
    Rater,
    D_ISC_Z,           # Fisher Z
    D_ENT_Z,           # Log Robust Z
    D_VAL_Z,           # Mean Robust Z
    H_POS_Z,           # Robust Z
    H_NEG_Z,           # Robust Z
    D_ISC,
    D_ENT,
    D_VAL,
    H_POS,
    H_NEG
  ) |>
  drop_na()

# 2. Correlation Matrix
cor_matrix <- 
  final_qc_prep |> 
  select(ends_with("_Z")) |> 
  cor()

# Display Correlations
cor_matrix |>
  as.data.frame() |>
  rownames_to_column("Metric") |>
  gt() |>
  fmt_number(columns = -Metric, decimals = 2) |>
  tab_header(title = "Correlation of QC Metrics")
Correlation of QC Metrics
Metric D_ISC_Z D_ENT_Z D_VAL_Z H_POS_Z H_NEG_Z
D_ISC_Z 1.00 −0.06 −0.25 −0.09 0.02
D_ENT_Z −0.06 1.00 0.00 0.13 0.26
D_VAL_Z −0.25 0.00 1.00 0.40 −0.37
H_POS_Z −0.09 0.13 0.40 1.00 0.37
H_NEG_Z 0.02 0.26 −0.37 0.37 1.00

6.2 Global Outlier Detection

We calculate the Global Mahalanobis Distance.

  • Pink cells indicate a potential failure.
  • Green cells indicate passing values.
# Define variables for MCD
mcd_vars <- c("D_ISC_Z", "D_ENT_Z", "D_VAL_Z", "H_POS_Z", "H_NEG_Z")

# 4. Calculate MCD and Distances
mcd_global <- 
  final_qc_prep |>
  select(all_of(mcd_vars)) |>
  cov.mcd()

final_qc <- 
  final_qc_prep |>
  mutate(
    Global_Dist = mahalanobis(
      x = pick(all_of(mcd_vars)),
      center = mcd_global$center,
      cov = mcd_global$cov
    ),
    Global_pval = pchisq(Global_Dist, df = length(mcd_vars), lower.tail = FALSE),
    # Flags
    Global_MD = Global_pval < .001,
    Low_ISC = D_ISC < 0.3,
    Low_ENT = D_ENT_Z < -2.5,
    High_ENT = D_ENT_Z > 2.5,
    Low_VAL = D_VAL_Z < -2.5,
    High_VAL = D_VAL_Z > 2.5,
    Low_H_POS = H_POS_Z < -2.5,
    High_H_POS = H_POS_Z > 2.5,
    Low_H_NEG = H_NEG_Z < -2.5,
    High_H_NEG = H_NEG_Z > 2.5,
    # Count Total Flags
    Flag_Count = rowSums(pick(Global_MD, Low_ISC, Low_ENT, High_ENT, 
                              Low_VAL, High_VAL,
                              Low_H_POS, High_H_POS, Low_H_NEG, High_H_NEG))
  )
# 5. Final Output Table
final_qc |> 
  filter(Flag_Count > 0) |> 
  relocate(Flag_Count, .after = Rater) |> 
  arrange(Rater) |> 
  select(
    Rater, Flag_Count, Global_MD, 
    Low_ISC, Low_ENT, High_ENT, 
    Low_VAL, High_VAL,
    Low_H_POS, High_H_POS, Low_H_NEG, High_H_NEG
  ) |>
  gt() |> 
  data_color(
    columns = where(is.logical),
    fn = function(x) ifelse(x == TRUE, "#ffcccc", "#ccffcc")
  ) |>
  tab_header(
    title = "Quality Control Summary",
    subtitle = "Flags for Synchrony, Complexity, Intensity, and Incoherence"
  )
Quality Control Summary
Flags for Synchrony, Complexity, Intensity, and Incoherence
Rater Flag_Count Global_MD Low_ISC Low_ENT High_ENT Low_VAL High_VAL Low_H_POS High_H_POS Low_H_NEG High_H_NEG
17 3 TRUE FALSE FALSE FALSE FALSE TRUE FALSE TRUE FALSE FALSE
29 1 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
40 1 FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
57 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
58 1 TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
59 1 FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
67 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
80 1 FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
81 4 TRUE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
99 1 FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE

7 Normative Baselines

To establish a reference standard for future data collection, we calculate the normative parameters (Median and MAD) and the multivariate covariance structure of the current sample. Crucially, we calculate these norms using only the “Valid” subset of participants (i.e., excluding those flagged for low synchrony, flat-lining, or data incoherence). This ensures our baseline represents “attentive human performance” without noise contamination.

# 1. Define the "Valid" Reference Set
# We exclude anyone with a QC flag to prevent outliers from contaminating the norms
valid_raters <- final_qc |> 
  filter(Flag_Count == 0) |> 
  select(Rater)

clean_data <- final_qc_prep |> 
  semi_join(valid_raters, by = "Rater")

# 2. Calculate Univariate Norms (on Clean Data Only)
normative_stats <- clean_data |> 
  summarize(
    # Synchrony (Fisher Z)
    D_ISC_Median = median(D_ISC, na.rm = TRUE),
    D_ISC_MAD    = mad(D_ISC, na.rm = TRUE),
    
    # Complexity (Log Entropy)
    D_ENT_Median = median(D_ENT, na.rm = TRUE),
    D_ENT_MAD    = mad(D_ENT, na.rm = TRUE),
    
    # Dynamic Valence
    D_VAL_Median = median(D_VAL, na.rm = TRUE),
    D_VAL_MAD    = mad(D_VAL, na.rm = TRUE),
    
    # Intensity Positive
    H_POS_Median = median(H_POS, na.rm = TRUE),
    H_POS_MAD    = mad(H_POS, na.rm = TRUE),
    
    # Intensity Negative
    H_NEG_Median = median(H_NEG, na.rm = TRUE),
    H_NEG_MAD    = mad(H_NEG, na.rm = TRUE)
  ) |> 
  pivot_longer(
    everything(), 
    names_to = "Temp", 
    values_to = "Value"
  ) |> 
  mutate(
    Stat = str_extract(Temp, "(Median|MAD)$"),
    Metric = str_remove(Temp, "_(Median|MAD)$")
  ) |> 
  select(-Temp) |>
  pivot_wider(names_from = Stat, values_from = Value)

# 3. Save & Display
write_rds(normative_stats, "normative_univariate_stats.rds")

normative_stats |> 
  gt() |> 
  tab_header(
    title = "Normative Parameters (Healthy Reference)",
    subtitle = "Derived from valid participant subset"
  ) |> 
  fmt_number(columns = c(Median, MAD), decimals = 3) |> 
  cols_label(
    Metric = "QC Metric",
    Median = "Reference Median",
    MAD = "Reference MAD"
  )
Normative Parameters (Healthy Reference)
Derived from valid participant subset
QC Metric Reference Median Reference MAD
D_ISC 0.722 0.126
D_ENT −2.419 0.497
D_VAL −0.025 0.437
H_POS 1.181 0.416
H_NEG 0.771 0.421 
# 4. Calculate Multivariate Norms (MCD on Clean Data)
# This provides the reference 'Cloud' for future Mahalanobis Distance checks
mcd_clean <- 
  clean_data |> 
  select(D_ISC, D_ENT, D_VAL, H_POS, H_NEG) |> 
  cov.mcd()

multivariate_norms <- list(
  center = mcd_clean$center,
  cov    = mcd_clean$cov
)

write_rds(multivariate_norms, "normative_multivariate_stats.rds")

# 5. Display Reference Correlation Matrix
# This confirms the relationships in the healthy baseline
cov2cor(mcd_clean$cov) |> 
  as.data.frame() |> 
  rownames_to_column("Metric") |> 
  gt() |> 
  tab_header(
    title = "Reference Correlation Structure",
    subtitle = "Correlations among QC metrics in valid sample"
  ) |> 
  fmt_number(columns = -Metric, decimals = 2)
Reference Correlation Structure
Correlations among QC metrics in valid sample
Metric D_ISC D_ENT D_VAL H_POS H_NEG
D_ISC 1.00 −0.04 −0.23 −0.02 0.19
D_ENT −0.04 1.00 0.26 0.30 0.23
D_VAL −0.23 0.26 1.00 0.35 −0.29
H_POS −0.02 0.30 0.35 1.00 0.48
H_NEG 0.19 0.23 −0.29 0.48 1.00

8 Strategic Selection for Clinical Discriminability

To optimize the detection of negative affective bias (a transdiagnostic marker of emotional disorders), stimulus selection must account for the psychometric properties of the healthy normative sample. Specifically, we aim to avoid ceiling effects that obscure clinical elevations while ensuring discriminant validity between appropriate and inappropriate negative responding. We stratify the stimuli into three functional categories based on their intensity and variance in the healthy sample:

  1. Ambiguous/Moderate Probes (Primary Diagnostic Utility)
    • Definition: Stimuli eliciting mild-to-moderate negative affect (\(M \approx 1.0-2.0\)) with high inter-individual variance (\(SD > 0.8\)).
    • Rationale: These stimuli possess the greatest discriminative power for negative bias. Because healthy controls do not rate these at the scale maximum, they provide sufficient psychometric “headroom” to detect hypersensitivity. Elevated ratings in this category reflect a tendency to interpret ambiguous or moderately negative cues as severe.
  2. High-Intensity Provocation (Reactivity & Regulation)
    • Definition: Stimuli eliciting near-maximal negative affect (\(M > 2.5\)) in the healthy sample.
    • Rationale: While these stimuli may suffer from ceiling effects regarding peak reactivity (limiting their ability to distinguish patients from controls on intensity alone), they remain essential for assessing emotion regulation and recovery. They serve as a “stress test” to observe the duration of the affective response post-stimulus.
  3. Low-Intensity Controls (Specificity)
    • Definition: Stimuli perceived as neutral or minimally negative (\(M < 0.8\)) by healthy controls.
    • Rationale: These stimuli establish the specificity of the negative bias. High negative ratings in this category suggest a decoupling of affect from context (e.g., hostile attribution bias or generalized negative affectivity) rather than a potentiated response to actual threat.

Table Interpretation:

  • Healthy_Mean: The normative baseline for negative affect. Lower values indicate greater capacity to detect clinical elevations.
  • Bias_Headroom: The remaining range on the Likert scale (Max - Healthy Mean). Larger values indicate greater sensitivity to exaggerated responses.
  • Clinical_Role: The proposed functional classification for the stimulus in a clinical battery.
# -------------------------------------------------------------------------
# Clinical Stimulus Selection (Negative Bias Focus)
# -------------------------------------------------------------------------

# 1. Calculate Reliability (Mean ISC per Video)
video_reliability <- val_df |>
  nest(.by = Abbrev) |>
  mutate(isc_data = map(data, calc_video_isc)) |> 
  select(Abbrev, isc_data) |>
  unnest(isc_data) |>
  summarize(
    Reliability_ISC = mean(r_val, na.rm = TRUE), 
    .by = Abbrev
  )

# 2. Calculate Discriminability (SD of Holistic Negative Affect)
video_discriminability <- hol_df |>
  filter(Scale == "Negative") |> 
  summarize(
    Mean_Neg = mean(Rating, na.rm = TRUE),
    SD_Neg   = sd(Rating, na.rm = TRUE), 
    .by = Abbrev
  )

# 3. Create Clinical Selection Table
clinical_selection <- 
  video_reliability |>
  left_join(video_discriminability, by = "Abbrev") |>
  mutate(
    # --- RENAME FOR CLARITY ---
    # Explicitly state this is Holistic Negative Affect
    Mean_Holistic_NA = Mean_Neg,
    SD_Holistic_NA   = SD_Neg,
    
    # 2. Assign Clinical Roles
    Clinical_Role = case_when(
      Reliability_ISC < 0.40 ~ "Exclusion (Low Reliability)",
      Mean_Holistic_NA > 2.0 ~ "High-Intensity Provocation",
      Mean_Holistic_NA >= 0.8 & Mean_Holistic_NA <= 2.0 & SD_Holistic_NA > 0.8 ~ "Ambiguous/Moderate Probe",
      Mean_Holistic_NA < 0.8 ~ "Low-Intensity Control",
      TRUE ~ "General Stimulus"
    ),
    
    # 3. "Headroom" Score
    Bias_Headroom = 4.0 - Mean_Holistic_NA
  ) |>
  select(Abbrev, Clinical_Role, Reliability_ISC, Mean_Holistic_NA, SD_Holistic_NA, Bias_Headroom) |> 
  arrange(Clinical_Role, Abbrev)

# 4. The Clinical Selection Table
clinical_selection |>
  gt() |>
  tab_header(
    title = "Stimulus Portfolio for Clinical Applications",
    subtitle = "Stratification by Intensity and Discriminative Potential"
  ) |>
  fmt_number(columns = where(is.numeric), decimals = 2) |>
  
  # --- IMPROVED COLUMN LABELS ---
  cols_label(
    Reliability_ISC = "Reliability (ISC)",
    Mean_Holistic_NA = "Mean Negative (Holistic)",
    SD_Holistic_NA = "SD Negative (Holistic)",
    Bias_Headroom = "Clinical Headroom"
  ) |>
  
  # Color code the roles
  data_color(
    columns = Clinical_Role,
    fn = function(x) {
      colors <- c(
        "Ambiguous/Moderate Probe" = "#cce5ff",
        "High-Intensity Provocation" = "#ffe5cc",
        "Low-Intensity Control" = "#e5ffcc",
        "Exclusion (Low Reliability)" = "#ffcccc",
        "General Stimulus" = "white"
      )
      colors[as.character(x)]
    }
  ) |>
  # Visualizing the Headroom
  data_color(
    columns = Bias_Headroom,
    fn = scales::col_numeric(palette = c("white", "grey"), domain = NULL)
  )
Stimulus Portfolio for Clinical Applications
Stratification by Intensity and Discriminative Potential
Abbrev Clinical_Role Reliability (ISC) Mean Negative (Holistic) SD Negative (Holistic) Clinical Headroom
Fences Ambiguous/Moderate Probe 0.77 1.79 1.24 2.21
GoodWill Ambiguous/Moderate Probe 0.61 0.97 0.99 3.03
MarriageStory Ambiguous/Moderate Probe 0.69 1.51 1.27 2.49
NoCountry Ambiguous/Moderate Probe 0.55 1.34 1.21 2.66
SilenceLambs Ambiguous/Moderate Probe 0.62 1.59 1.16 2.41
Zodiac Ambiguous/Moderate Probe 0.82 1.89 1.29 2.11
DaysSummer Exclusion (Low Reliability) 0.37 0.43 0.75 3.57
ForrestGump Exclusion (Low Reliability) 0.37 0.60 0.95 3.40
LittleMiss Exclusion (Low Reliability) 0.32 0.57 0.89 3.43
PulpFiction Exclusion (Low Reliability) 0.40 0.86 1.03 3.14
GreenMile High-Intensity Provocation 0.68 2.19 1.39 1.81
AkeelahBee Low-Intensity Control 0.86 0.37 0.76 3.63
BealeStreet Low-Intensity Control 0.50 0.66 0.91 3.34
CatchMe Low-Intensity Control 0.65 0.46 0.78 3.54
KingsSpeech Low-Intensity Control 0.43 0.74 1.02 3.26
LadyBird Low-Intensity Control 0.55 0.41 0.72 3.59
LegallyBlonde Low-Intensity Control 0.57 0.26 0.54 3.74
Miracle Low-Intensity Control 0.88 0.31 0.67 3.69
Moonlight Low-Intensity Control 0.61 0.47 0.74 3.53
ParentTrap Low-Intensity Control 0.56 0.25 0.62 3.75
PursuitHappyness Low-Intensity Control 0.76 0.34 0.67 3.66
SocialNetwork Low-Intensity Control 0.50 0.59 0.90 3.41