From d58e741f5db07c09e72487fc30eb8bf5cb8115ff Mon Sep 17 00:00:00 2001 From: Orihara Izaya Date: Wed, 30 Sep 2026 10:09:34 +0200 Subject: [PATCH] initial commit --- README.md | 6 + .../archivegpt_reliability/README.md | 234 ++++++++++ .../archivegpt_reliability/config.yaml | 98 +++++ .../archivegpt_reliability/requirements.txt | 9 + .../archivegpt_reliability/run_all.py | 109 +++++ .../archivegpt_reliability/src/__init__.py | 2 + .../src/build_dataset.py | 211 +++++++++ .../archivegpt_reliability/src/config.py | 162 +++++++ .../archivegpt_reliability/src/diagnostics.py | 416 ++++++++++++++++++ .../archivegpt_reliability/src/icc.py | 277 ++++++++++++ .../archivegpt_reliability/src/report.py | 302 +++++++++++++ .../archivegpt_reliability/src/validate.py | 186 ++++++++ 12 files changed, 2012 insertions(+) create mode 100644 README.md create mode 100644 archivegpt_reliability/archivegpt_reliability/README.md create mode 100644 archivegpt_reliability/archivegpt_reliability/config.yaml create mode 100644 archivegpt_reliability/archivegpt_reliability/requirements.txt create mode 100644 archivegpt_reliability/archivegpt_reliability/run_all.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/__init__.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/build_dataset.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/config.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/diagnostics.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/icc.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/report.py create mode 100644 archivegpt_reliability/archivegpt_reliability/src/validate.py diff --git a/README.md b/README.md new file mode 100644 index 0000000..083a201 --- /dev/null +++ b/README.md @@ -0,0 +1,6 @@ +# project_ArchiveGPT_scripts + +This repository contains non-main analysis for the *ArchiveGPT* project. Each folder contains its own analysis. + +The following analyses are contained: +- archivegpt_reliability: Analysis of interrater reliability, due to a question from a researcher. \ No newline at end of file diff --git a/archivegpt_reliability/archivegpt_reliability/README.md b/archivegpt_reliability/archivegpt_reliability/README.md new file mode 100644 index 0000000..6497674 --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/README.md @@ -0,0 +1,234 @@ +# ArchiveGPT inter-rater reliability + +This is an additional reproducible analysis of how far raters agreed with one another about the quality +of the descriptions in the [ArchiveGPT study](https://doi.org/10.1057/s41599-026-08367-6), and of what the reported +single-rater intraclass correlation of ~0.11 does and does not establish. + +It is accompanied by the technical note `output/reliability_report.txt` — every result in that note is +produced by `run_all.py`. + +--- + +## Why this code starts from the raw exports + +The published file `analysis/data/clean_data/correctness_per_trial.csv` assigns +its `trial` column as `rep(1:80, ...)` — serial position (line 230 of +`01_pre-processing.R`). Presentation order was **randomised per participant**, so +that column does not identify the stimulus. No item-level reliability analysis +can be run from the published clean data alone. + +Item identity survives only in the raw PsychoPy exports, in the `image` and +`desc_type` columns. This package rebuilds it: + +``` +item = image + "|" + desc_type # 40 images x 2 sources = 80 items +``` + +Each participant rated all 80 items across two blocks of 40, with the description +source flipping between blocks. The design is fully crossed: 139 participants +x 80 items = 11,120 ratings, one rating per cell, no missing data. + +**Fixing this in the release** — adding `image` and `desc_type` to +`correctness_per_trial.csv` — would make this reconstruction step unnecessary and +close a genuine reproducibility gap. + +--- + +## Quick start + +```bash +pip install -r requirements.txt + +# 1. Point config.yaml at your copy of the extracted ArchiveGPT release, +# or override it on the command line: +python run_all.py --data-root /path/to/archive_gpt +``` + +That writes three files to `output/`: + +| file | contents | +|---|---| +| `reliability_report.txt` | every table in the note, as plain text | +| `results.json` | the same numbers, machine-readable | +| `ratings_long.csv` | the rebuilt trial-level data **with item identity** | + +`ratings_long.csv` is the artefact worth keeping: it is what the published clean +data would look like if it carried item identity. + +--- + +## Settings + +Everything lives in **`config.yaml`**. The only value you normally change is: + +```yaml +paths: + data_root: "D:/McDuf/Downloads/archive_gpt" +``` + +This must point at the extracted Zenodo release ([record 15640284](https://zenodo.org/records/15640284)) — the folder +containing `analysis/data/raw_data/` and `analysis/data/clean_data/`. Forward +slashes work on Windows. Use `"."` if you drop this package inside the release. + +Other sections, all documented inline: + +- **`columns`** — column names in the raw exports. The two slider columns per + measure (`slider_q1.response`, `slider_q1_2.response`) correspond to the two + presentation blocks. +- **`design`** — what an item is, how many to expect, which `desc_type` value + means archivist-written vs model-generated, which `expert_status` value marks + each panel. +- **`analysis`** — `primary_measure` (`accuracy` or `usefulness`), `alpha`, + Monte-Carlo `homogeneity_sims`, `seed`, and how many weak raters to drop in + the final adjustment rung. +- **`validation`** — whether to abort when the reconstruction disagrees with the + published data. **Leave `enforce: true`.** + +CLI flags override the file: `--data-root`, `--output`, `--measure`, `--quick` +(200 simulations instead of 2000), `--no-validate`. + +--- + +## What gets checked before anything is computed + +The whole analysis rests on having reassembled item identity correctly, so that +is verified against the release before any coefficient is reported: + +| check | what it proves | +|---|---| +| fully crossed design | 139 x 80 ratings, exactly one per participant x item | +| two balanced blocks | 40 items per block per participant | +| counterbalancing exact | description source flips between blocks for every image | +| ratings match `correctness_per_trial.csv` | max abs difference **0.0** | +| means match `accuracy_usefulness.csv` | max abs difference **1.4e-14** | +| independent ICC cross-check (pingouin) | max discrepancy **< 1e-15** | + +The first comparison is on each participant's *sorted* ratings, because the +published `trial` column carries no item identity to align on. Passing both means +the reconstruction contains precisely the published ratings, and the only thing +added is the item label. + +The pingouin cross-check is skipped with a notice if pingouin is not installed. +Everything else works without it. + +--- + +## Layout + +``` +archivegpt_reliability/ +├── config.yaml <- the only file you normally edit +├── requirements.txt +├── run_all.py <- entry point +├── README.md +└── src/ + ├── config.py loads and validates settings, resolves paths + ├── build_dataset.py rebuilds item identity from the raw exports + ├── validate.py checks the rebuild against the published data + ├── icc.py variance components + the whole ICC family + CIs + ├── diagnostics.py every Part II diagnostic + └── report.py renders the tables, writes the outputs +``` + +### `icc.py` + +Implements the two-way random-effects decomposition + +``` +X_ij = mu + t_i + r_j + e_ij +``` + +and every coefficient derived from it. Confidence intervals follow McGraw & +Wong (1996) Table 7, and are validated against pingouin. + +| function | purpose | +|---|---| +| `variance_components` | sigma^2_item, sigma^2_rater, sigma^2_error + mean squares | +| `icc_estimates` | ICC(1,1), (A,1), (C,1) and their k-rater forms | +| `icc_confidence_intervals` | exact-F and Satterthwaite intervals | +| `single_rater_validity` | sqrt(ICC(C,1)) — one rater's correlation with truth | +| `descriptive_scale` | the components translated into rating points | +| `spearman_brown`, `raters_needed` | panel-size planning | +| `crosscheck_pingouin` | independent re-estimation | + +### `diagnostics.py` + +| function | question it answers | +|---|---| +| `pairwise_stats` | distribution of rater-pair correlations | +| `homogeneity_test` | are raters interchangeable, or does quality genuinely differ? | +| `rater_validity` | each rater vs the rest, disattenuated | +| `severity_stats` | how differently raters use the scale | +| `cumulative_adjustments` | the coefficient as design effects are modelled | +| `block_effects` | first-block to second-block drift | +| `source_split` | agreement **and** true spread, by description source | +| `range_restriction_counterfactual` | what a lower ICC in a subset actually means | +| `cross_panel_agreement` | do the two panels rank the descriptions alike? | +| `substantive_effect` | the human-vs-model comparison, items fixed then random | +| `replication_available` | did anyone rate an item twice? | +| `image_taste_test` | is a rater's deviation stable across the same image? | +| `residual_split_consequences` | what turns on the unidentified residual split | +| `within_rater_measure_coherence` | is each rater internally self-consistent? | + +--- + +## Three things the code is careful about + +**1. Which axis a "within-item SD" is taken over.** `data.std(axis=1)` is the SD +across raters within an item; `axis=0` is the SD across items within a rater. +They are different quantities and easy to confuse — an earlier draft of the note +reported the second under the first's name. `source_split` now returns both, +named explicitly. + +**2. Comparing like with like in the counterfactual.** The range-restriction +counterfactual must be built in the same ICC form as the observed value it is +compared against. Substituting a consistency-form counterfactual for an +absolute-agreement observation silently overstates the gap (0.171 instead of +0.105). + +**3. The block-order artifact in `image_taste_test`.** The two descriptions of an +image necessarily fall in different blocks, so the ~4-point block drift induces a +spurious *negative* correlation if it is not removed first. Run with +`control_block=False` to see the artifact (median r = -0.053, p = 0.014) and +`True` for the real answer (median r = -0.017, p = 0.293). `run_all.py` prints +both. + +--- + +## Known limitations + +- **No replication.** Exactly one rating per participant x item, so the residual + confounds stable rater x item interaction with momentary noise. They are not + identified, and the interpretation of the whole result turns on the split. + `residual_split_consequences` shows the range; only a retest resolves it. +- **The two-way model has no block term.** The ~4-point block drift is absorbed + into the residual. `cumulative_adjustments` shows the effect of removing it + (0.156 -> 0.165). +- **`image_taste_test` is a proxy, not a retest.** It rules out stable taste for + the depicted *image*; it says nothing about taste for the specific + *description*, which is the component that matters. +- **Sample.** 29 participants carry `expert_status == "expert"`; an external + reanalysis reported 28. If you are comparing figures directly, resolve the + one-participant discrepancy first. +- **"Expert" is heterogeneous**: 14 archaeology or history students, 9 archive + students, 4 archivists, 1 archaeologist, 1 hobby historian. + +--- + +## Requirements + +Python 3.9+. `numpy`, `pandas`, `scipy`, `PyYAML`. Optionally `pingouin` for the +independent ICC cross-check. + +## References + +- Shrout, P. E., & Fleiss, J. L. (1979). *Psychological Bulletin*, 86(2), 420–428. +- McGraw, K. O., & Wong, S. P. (1996). *Psychological Methods*, 1(1), 30–46. +- Koo, T. K., & Li, M. Y. (2016). *Journal of Chiropractic Medicine*, 15(2), 155–163. +- Bornmann, L., Mutz, R., & Daniel, H.-D. (2010). *PLoS ONE*, 5(12), e14331. +- Erosheva, E. A., Martinková, P., & Lee, C. J. (2021). *JRSS-A*, 184(3), 904–919. + +## Data + +Abele, L., Anders, G., Aydın, T., Buder, J., Fischer, H., Kimmel, D., & Huff, M. +*ArchiveGPT.* Zenodo, record 15640284. https://zenodo.org/records/15640284 diff --git a/archivegpt_reliability/archivegpt_reliability/config.yaml b/archivegpt_reliability/archivegpt_reliability/config.yaml new file mode 100644 index 0000000..9a669df --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/config.yaml @@ -0,0 +1,98 @@ +# ============================================================================== +# ArchiveGPT inter-rater reliability -- analysis settings +# ============================================================================== +# The only value you normally need to change is paths.data_root. +# +# Everything else describes the structure of the ArchiveGPT release +# (Zenodo record 15640284) and the parameters of the analysis. Override any +# value at the command line with, e.g.: +# python run_all.py --data-root /path/to/archive_gpt +# ============================================================================== + +paths: + # --------------------------------------------------------------------------- + # Root of the extracted ArchiveGPT release. Must contain + # analysis/data/raw_data/ (the per-session PsychoPy exports) + # analysis/data/clean_data/ (the published aggregates) + # Forward slashes work on Windows too. Use "." if this package sits inside + # the release folder. + # --------------------------------------------------------------------------- + data_root: "D:/McDuf/Downloads/archive_gpt" + + # Locations within data_root. Change only if you reorganised the release. + raw_data_dir: "analysis/data/raw_data" + clean_data_dir: "analysis/data/clean_data" + + # Published clean-data files. These are used ONLY to (a) define which + # participants were retained in the published sample and (b) validate the + # reconstruction. The reliability analysis itself reads the raw exports, + # because the published trial-level file does not carry item identity. + participants_file: "dprimes_willingness_trust_expertise.csv" + trial_file: "correctness_per_trial.csv" + aggregate_file: "accuracy_usefulness.csv" + + # Results are written here. Relative paths resolve against this config file. + output_dir: "output" + +columns: + # Column names in the raw PsychoPy exports. + participant: "Participant" + image: "image" + desc_type: "desc_type" + + # The experiment presents 80 descriptions in two blocks of 40. Each block + # writes to its own slider column; a trial has exactly one of them non-null. + # List them in block order. + accuracy: ["slider_q1.response", "slider_q1_2.response"] + usefulness: ["slider_q2.response", "slider_q2_2.response"] + + # Column names in the published participants file. + expert_status: "expert_status" + expertise_detail: "archive_expertise_status" + +design: + # An ITEM is one description: one image shown with one description source. + # 40 images x 2 sources = 80 items. Presentation order was randomised, so the + # `trial` column in the published clean data is serial position, NOT item + # identity. This is why the analysis must start from the raw exports. + n_images_expected: 40 + n_items_expected: 80 + + # Values of desc_type identifying who wrote the description. + source_human: "expert" # written by an archivist + source_model: "ai" # generated by InternVL2 + + # Values of expert_status identifying the two rater panels. + panel_expert: "expert" + panel_nonexpert: "non_expert" + + # Slider bounds, used for reporting only. + rating_scale: [0, 100] + +analysis: + # Which rating to treat as primary. "usefulness" reruns everything on the + # second slider. + primary_measure: "accuracy" + + # Confidence level for all interval estimates. + alpha: 0.05 + + # Monte Carlo replications for the rater-homogeneity test (Section 5 of the + # note). 2000 is ample; 200 is enough for a quick look. + homogeneity_sims: 2000 + + # Seed for every stochastic step, so runs are reproducible. + seed: 20240607 + + # How many lowest-validity raters to drop in the final rung of the + # cumulative-adjustment ladder (Section 2 of the note). + drop_weakest_n: 3 + +validation: + # Abort if the reconstruction disagrees with the published clean data. + # Leave this on: it is the check that proves the rebuilt item identity is + # correct. + enforce: true + + # Maximum tolerated absolute difference against the published files. + tolerance: 1.0e-6 diff --git a/archivegpt_reliability/archivegpt_reliability/requirements.txt b/archivegpt_reliability/archivegpt_reliability/requirements.txt new file mode 100644 index 0000000..2a85b83 --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/requirements.txt @@ -0,0 +1,9 @@ +# Core (required) +numpy>=1.23 +pandas>=1.5 +scipy>=1.9 +PyYAML>=6.0 + +# Optional: enables the independent ICC cross-check in run_all.py. +# Everything works without it; the check is simply skipped. +pingouin>=0.5.3 diff --git a/archivegpt_reliability/archivegpt_reliability/run_all.py b/archivegpt_reliability/archivegpt_reliability/run_all.py new file mode 100644 index 0000000..b65acc1 --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/run_all.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +"""Reproduce the full ArchiveGPT inter-rater reliability analysis. + + python run_all.py # uses config.yaml as written + python run_all.py --data-root /path/to/archive_gpt + python run_all.py --measure usefulness + python run_all.py --quick # fewer simulations, for a fast look + +Writes reliability_report.txt, results.json and ratings_long.csv to the output +directory named in config.yaml. +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +sys.path.insert(0, str(Path(__file__).resolve().parent)) + +from src.build_dataset import ReconstructionError, build_long_frame, rating_matrix +from src.config import ConfigError, load_config +from src.icc import crosscheck_pingouin +from src.report import build_report, write_outputs +from src.validate import ValidationError, validate_reconstruction + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser( + description="ArchiveGPT inter-rater reliability analysis.", + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + parser.add_argument("--config", default=None, help="settings file (default: config.yaml)") + parser.add_argument("--data-root", default=None, help="override paths.data_root") + parser.add_argument("--output", default=None, help="override paths.output_dir") + parser.add_argument( + "--measure", + choices=["accuracy", "usefulness"], + default=None, + help="override analysis.primary_measure", + ) + parser.add_argument( + "--quick", action="store_true", help="use 200 homogeneity simulations instead of the configured number" + ) + parser.add_argument("--no-validate", action="store_true", help="skip the validation checks") + return parser.parse_args() + + +def main() -> int: + args = parse_args() + + try: + cfg = load_config(args.config, args.data_root, args.output) + except ConfigError as exc: + print(f"\nConfiguration problem:\n {exc}\n", file=sys.stderr) + return 2 + + if args.measure: + cfg.raw["analysis"]["primary_measure"] = args.measure + if args.quick: + cfg.raw["analysis"]["homogeneity_sims"] = 200 + if args.no_validate: + cfg.raw["validation"]["enforce"] = False + + print(cfg.describe()) + + print("\nRebuilding item identity from the raw exports ...") + try: + long = build_long_frame(cfg) + except ReconstructionError as exc: + print(f"\nReconstruction failed:\n {exc}\n", file=sys.stderr) + return 3 + + print("\nValidating against the published clean data ...") + try: + report = validate_reconstruction(long, cfg) + except ValidationError as exc: + print(f"\n{exc}\n", file=sys.stderr) + return 4 + print(report.render()) + + matrix = rating_matrix( + long, measure=cfg.analysis["primary_measure"], panel=cfg.design["panel_expert"] + ) + discrepancies = crosscheck_pingouin(matrix.to_numpy()) + if discrepancies is None: + print("\n [SKIP] pingouin not installed; ICC cross-check not run " + "(pip install pingouin)") + else: + worst = max(discrepancies.values()) + status = "PASS" if worst < 1e-9 else "FAIL" + print(f"\n [{status}] independent ICC cross-check (pingouin): " + f"max discrepancy {worst:.2e}") + if status == "FAIL": + return 5 + + print("\nRunning the analysis ...") + text, results = build_report(long, cfg) + print(text) + + written = write_outputs(text, results, long, cfg) + print(f"\n{'=' * 78}\nWrote:") + for path in written: + print(f" {path}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/archivegpt_reliability/archivegpt_reliability/src/__init__.py b/archivegpt_reliability/archivegpt_reliability/src/__init__.py new file mode 100644 index 0000000..c9aa282 --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/__init__.py @@ -0,0 +1,2 @@ +"""ArchiveGPT inter-rater reliability analysis.""" +__version__ = "1.0.0" diff --git a/archivegpt_reliability/archivegpt_reliability/src/build_dataset.py b/archivegpt_reliability/archivegpt_reliability/src/build_dataset.py new file mode 100644 index 0000000..01f1f29 --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/build_dataset.py @@ -0,0 +1,211 @@ +"""Rebuild trial-level ratings with item identity from the raw PsychoPy exports. + +Why this module exists +---------------------- +The published file ``clean_data/correctness_per_trial.csv`` assigns its ``trial`` +column as ``rep(1:80, ...)`` -- serial position. Because presentation order was +randomised per participant, that column does NOT identify the stimulus, and no +item-level reliability analysis can be run from the published clean data alone. + +Item identity survives only in the raw exports, in the ``image`` and +``desc_type`` columns. An item is one description: + + item = image + "|" + desc_type (40 images x 2 sources = 80 items) + +Each participant sees all 80 across two blocks of 40, with the description +source flipping between blocks. The design is fully crossed. +""" + +from __future__ import annotations + +import re +from pathlib import Path + +import numpy as np +import pandas as pd + +from .config import Config + +VERSION_PATTERN = re.compile(r"Exp1b?_(v[0-9a-z]+)", re.IGNORECASE) + + +class ReconstructionError(RuntimeError): + """Raised when the rebuilt dataset does not have the expected structure.""" + + +def _session_version(filename: str) -> str: + match = VERSION_PATTERN.search(filename) + return match.group(1) if match else "unknown" + + +def _read_session(path: Path, cfg: Config) -> pd.DataFrame | None: + """Extract one session's trials as long rows, or None if it holds no data. + + Partial and aborted sessions are common in the release (many raw files lack + the slider columns entirely); those are skipped here and, in any case, + excluded later by the restriction to the published sample. + """ + cols = cfg.columns + try: + frame = pd.read_csv(path, low_memory=False) + except Exception: + return None + + needed = {cols["image"], cols["desc_type"], cols["participant"]} + if not needed.issubset(frame.columns): + return None + + participants = frame[cols["participant"]].dropna() + if participants.empty: + return None + participant = participants.iloc[0] + + acc_cols = cfg.measure_columns("accuracy") + use_cols = cfg.measure_columns("usefulness") + version = _session_version(path.name) + + blocks = [] + for block, (acc_col, use_col) in enumerate(zip(acc_cols, use_cols), start=1): + if acc_col not in frame.columns: + continue + rows = frame[frame[acc_col].notna()] + if rows.empty: + continue + block_frame = pd.DataFrame( + { + "participant": participant, + "version": version, + "block": block, + "image": rows[cols["image"]].astype(str).values, + "desc_type": rows[cols["desc_type"]].astype(str).values, + "accuracy": rows[acc_col].astype(float).values, + "usefulness": ( + rows[use_col].astype(float).values + if use_col in frame.columns + else np.nan + ), + "source_file": path.name, + } + ) + blocks.append(block_frame) + + if not blocks: + return None + return pd.concat(blocks, ignore_index=True) + + +def build_long_frame(cfg: Config, verbose: bool = True) -> pd.DataFrame: + """Read every raw export and return one row per participant x item. + + The returned frame is restricted to the participants retained in the + published sample and carries their panel membership. + """ + files = sorted(cfg.raw_data_dir.glob("*.csv")) + if not files: + raise ReconstructionError(f"No CSV files found in {cfg.raw_data_dir}") + + sessions = [frame for frame in (_read_session(path, cfg) for path in files) if frame is not None] + if not sessions: + raise ReconstructionError( + f"None of the {len(files)} files in {cfg.raw_data_dir} contained " + f"usable trials. Check columns.* in config.yaml." + ) + + long = pd.concat(sessions, ignore_index=True) + long = long.dropna(subset=["participant", "image", "desc_type", "accuracy"]) + long["participant"] = long["participant"].astype(int) + long["item"] = long["image"] + "|" + long["desc_type"] + + if verbose: + print( + f" read {len(files)} raw files -> {len(sessions)} with trials, " + f"{len(long):,} ratings, {long['participant'].nunique()} participants" + ) + + long = _restrict_to_published_sample(long, cfg, verbose=verbose) + _check_structure(long, cfg) + return long + + +def _restrict_to_published_sample( + long: pd.DataFrame, cfg: Config, verbose: bool = True +) -> pd.DataFrame: + """Keep only participants in the published sample and attach panel labels.""" + cols = cfg.columns + published = pd.read_csv(cfg.participants_file) + + keep_cols = [cols["participant"], cols["expert_status"]] + if cols["expertise_detail"] in published.columns: + keep_cols.append(cols["expertise_detail"]) + + published = published[keep_cols].rename( + columns={ + cols["participant"]: "participant", + cols["expert_status"]: "panel", + cols["expertise_detail"]: "expertise_detail", + } + ) + + merged = long.merge(published, on="participant", how="inner") + if verbose: + counts = merged.groupby("panel")["participant"].nunique().to_dict() + print( + f" restricted to published sample: {len(merged):,} ratings, " + f"{merged['participant'].nunique()} participants {counts}" + ) + return merged + + +def _check_structure(long: pd.DataFrame, cfg: Config) -> None: + """Assert the rebuilt design is what the experiment specifies.""" + design = cfg.design + + n_images = long["image"].nunique() + n_items = long["item"].nunique() + if n_images != design["n_images_expected"]: + raise ReconstructionError( + f"Expected {design['n_images_expected']} images, found {n_images}." + ) + if n_items != design["n_items_expected"]: + raise ReconstructionError( + f"Expected {design['n_items_expected']} items " + f"(image x desc_type), found {n_items}." + ) + + per_cell = long.groupby(["participant", "item"]).size() + if per_cell.max() != 1: + raise ReconstructionError( + f"Design is not one rating per participant x item " + f"(max = {per_cell.max()}). Repeated measurements would change " + f"the model in Part I of the note." + ) + + per_participant = long.groupby("participant")["item"].nunique() + incomplete = per_participant[per_participant != design["n_items_expected"]] + if not incomplete.empty: + raise ReconstructionError( + f"{len(incomplete)} participants do not have all " + f"{design['n_items_expected']} items: " + f"{incomplete.head().to_dict()}" + ) + + +def rating_matrix( + long: pd.DataFrame, + measure: str = "accuracy", + panel: str | None = None, + desc_type: str | None = None, +) -> pd.DataFrame: + """Pivot to an items x raters matrix, dropping any incomplete rows. + + Returns a DataFrame so item and participant labels survive; call + ``.to_numpy()`` for the array the ICC routines expect. + """ + subset = long + if panel is not None: + subset = subset[subset["panel"] == panel] + if desc_type is not None: + subset = subset[subset["desc_type"] == desc_type] + + matrix = subset.pivot_table(index="item", columns="participant", values=measure) + return matrix.dropna() diff --git a/archivegpt_reliability/archivegpt_reliability/src/config.py b/archivegpt_reliability/archivegpt_reliability/src/config.py new file mode 100644 index 0000000..964fe0e --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/config.py @@ -0,0 +1,162 @@ +"""Load, resolve and validate the analysis settings in config.yaml.""" + +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +try: + import yaml +except ImportError as exc: # pragma: no cover + raise SystemExit( + "PyYAML is required to read config.yaml.\n" + " pip install pyyaml (or: pip install -r requirements.txt)" + ) from exc + +PACKAGE_ROOT = Path(__file__).resolve().parent.parent +DEFAULT_CONFIG = PACKAGE_ROOT / "config.yaml" + + +class ConfigError(RuntimeError): + """Raised when the settings file is unusable or points somewhere wrong.""" + + +@dataclass +class Config: + """Resolved settings. Paths are absolute; everything else is as written.""" + + raw: dict[str, Any] + config_path: Path + + data_root: Path + raw_data_dir: Path + clean_data_dir: Path + participants_file: Path + trial_file: Path + aggregate_file: Path + output_dir: Path + + # --- convenience accessors ------------------------------------------------- + @property + def columns(self) -> dict[str, Any]: + return self.raw["columns"] + + @property + def design(self) -> dict[str, Any]: + return self.raw["design"] + + @property + def analysis(self) -> dict[str, Any]: + return self.raw["analysis"] + + @property + def validation(self) -> dict[str, Any]: + return self.raw["validation"] + + def measure_columns(self, measure: str) -> list[str]: + """Slider columns for 'accuracy' or 'usefulness', in block order.""" + try: + return list(self.columns[measure]) + except KeyError as exc: + raise ConfigError( + f"Unknown measure {measure!r}. " + f"Expected one of: accuracy, usefulness." + ) from exc + + def describe(self) -> str: + return ( + f"config : {self.config_path}\n" + f"data_root : {self.data_root}\n" + f"raw exports : {self.raw_data_dir}\n" + f"clean data : {self.clean_data_dir}\n" + f"output : {self.output_dir}\n" + f"measure : {self.analysis['primary_measure']}\n" + f"seed : {self.analysis['seed']}" + ) + + +def _resolve(base: Path, value: str) -> Path: + """Resolve `value` against `base` unless it is already absolute.""" + path = Path(str(value)).expanduser() + return path if path.is_absolute() else (base / path) + + +def load_config( + config_path: str | Path | None = None, + data_root_override: str | Path | None = None, + output_dir_override: str | Path | None = None, +) -> Config: + """Read config.yaml, resolve every path, and check the data is reachable. + + Parameters + ---------- + config_path + Path to the settings file. Defaults to config.yaml beside this package. + data_root_override + If given, replaces paths.data_root (the --data-root CLI flag). + output_dir_override + If given, replaces paths.output_dir (the --output CLI flag). + """ + config_path = Path(config_path or DEFAULT_CONFIG).expanduser().resolve() + if not config_path.is_file(): + raise ConfigError(f"Settings file not found: {config_path}") + + with config_path.open("r", encoding="utf-8") as handle: + raw = yaml.safe_load(handle) + + for section in ("paths", "columns", "design", "analysis", "validation"): + if section not in raw: + raise ConfigError(f"config.yaml is missing the '{section}' section.") + + paths = raw["paths"] + here = config_path.parent + + data_root = _resolve(here, data_root_override or paths["data_root"]).resolve() + if not data_root.is_dir(): + raise ConfigError( + f"data_root does not exist: {data_root}\n" + f"Edit paths.data_root in {config_path}, or pass --data-root." + ) + + raw_data_dir = _resolve(data_root, paths["raw_data_dir"]) + clean_data_dir = _resolve(data_root, paths["clean_data_dir"]) + + for label, directory in (("raw_data_dir", raw_data_dir), ("clean_data_dir", clean_data_dir)): + if not directory.is_dir(): + raise ConfigError( + f"{label} does not exist: {directory}\n" + f"Does data_root point at the extracted ArchiveGPT release?" + ) + + participants_file = _resolve(clean_data_dir, paths["participants_file"]) + trial_file = _resolve(clean_data_dir, paths["trial_file"]) + aggregate_file = _resolve(clean_data_dir, paths["aggregate_file"]) + + if not participants_file.is_file(): + raise ConfigError( + f"Participants file not found: {participants_file}\n" + f"This file defines the retained sample and is required." + ) + + output_dir = _resolve(here, output_dir_override or paths["output_dir"]) + output_dir.mkdir(parents=True, exist_ok=True) + + measure = raw["analysis"].get("primary_measure", "accuracy") + if measure not in raw["columns"]: + raise ConfigError( + f"analysis.primary_measure is {measure!r}, which has no entry under " + f"'columns'. Expected 'accuracy' or 'usefulness'." + ) + + return Config( + raw=raw, + config_path=config_path, + data_root=data_root, + raw_data_dir=raw_data_dir, + clean_data_dir=clean_data_dir, + participants_file=participants_file, + trial_file=trial_file, + aggregate_file=aggregate_file, + output_dir=output_dir, + ) diff --git a/archivegpt_reliability/archivegpt_reliability/src/diagnostics.py b/archivegpt_reliability/archivegpt_reliability/src/diagnostics.py new file mode 100644 index 0000000..d2a849b --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/diagnostics.py @@ -0,0 +1,416 @@ +"""Every diagnostic reported in Part II of the technical note. + +Each function is independent and takes either the long frame or an items x +raters matrix, so they can be run piecemeal from a notebook. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from scipy import stats + +from .config import Config +from .icc import icc_estimates, variance_components + + +# --------------------------------------------------------------------------- +# Pairwise agreement and whether the panel is homogeneous +# --------------------------------------------------------------------------- +def pairwise_stats(matrix) -> dict[str, float]: + """Distribution of the correlations between every pair of raters.""" + data = np.asarray(matrix, dtype=float) + corr = np.corrcoef(data.T) + pairs = corr[np.triu_indices(data.shape[1], k=1)] + return { + "n_pairs": int(pairs.size), + "median": float(np.median(pairs)), + "mean": float(pairs.mean()), + "sd": float(pairs.std(ddof=1)), + "min": float(pairs.min()), + "max": float(pairs.max()), + "share_below_0.10": float((pairs < 0.10).mean()), + "share_negative": float((pairs < 0).mean()), + } + + +def homogeneity_test(matrix, n_sims: int = 2000, seed: int = 0) -> dict[str, float]: + """Are the raters interchangeable, or does their quality genuinely differ? + + Simulates panels of interchangeable raters with the observed variance + components and compares the dispersion of pairwise correlations. If the + observed dispersion exceeds the simulated range, rater quality is + heterogeneous and a single pooled coefficient hides real structure. + """ + comp = variance_components(matrix) + observed = pairwise_stats(matrix)["sd"] + rng = np.random.default_rng(seed) + n, k = comp.n_items, comp.n_raters + sd_item = np.sqrt(max(comp.var_item, 0.0)) + sd_error = np.sqrt(comp.var_error) + + simulated = np.empty(n_sims) + for i in range(n_sims): + truth = rng.normal(0, sd_item, (n, 1)) + noise = rng.normal(0, sd_error, (n, k)) + corr = np.corrcoef((truth + noise).T) + simulated[i] = corr[np.triu_indices(k, k=1)].std(ddof=1) + + lo, hi = np.percentile(simulated, [2.5, 97.5]) + excess = max(observed**2 - simulated.mean() ** 2, 0.0) + return { + "observed_sd": float(observed), + "expected_sd": float(simulated.mean()), + "expected_lo": float(lo), + "expected_hi": float(hi), + "heterogeneous": bool(observed > hi), + "implied_true_sd": float(np.sqrt(excess)), + "n_sims": int(n_sims), + } + + +# --------------------------------------------------------------------------- +# Individual raters +# --------------------------------------------------------------------------- +def rater_validity(matrix) -> pd.DataFrame: + """Each rater's correlation with the mean of all others, disattenuated. + + The comparison mean is itself unreliable, so the raw correlation understates + a rater's validity. Dividing by the square root of that mean's reliability + gives the implied correlation with the item's true score. + """ + data = np.asarray(matrix, dtype=float) + k = data.shape[1] + c1 = icc_estimates(data)["ICC(C,1)"] + reliability_rest = (k - 1) * c1 / (1 + (k - 2) * c1) if c1 > 0 else np.nan + + raw = np.array( + [np.corrcoef(data[:, j], np.delete(data, j, axis=1).mean(axis=1))[0, 1] for j in range(k)] + ) + validity = raw / np.sqrt(reliability_rest) if reliability_rest > 0 else np.full(k, np.nan) + + labels = matrix.columns if isinstance(matrix, pd.DataFrame) else np.arange(k) + return pd.DataFrame( + {"rater": labels, "corr_with_rest": raw, "implied_validity": validity} + ).sort_values("corr_with_rest", ascending=False, ignore_index=True) + + +def severity_stats(matrix) -> dict[str, float]: + """How differently raters use the scale, independent of their ordering.""" + data = np.asarray(matrix, dtype=float) + means = data.mean(axis=0) + sds = data.std(axis=0, ddof=1) + return { + "rater_mean_min": float(means.min()), + "rater_mean_median": float(np.median(means)), + "rater_mean_max": float(means.max()), + "rater_mean_spread": float(means.max() - means.min()), + "rater_mean_sd": float(means.std(ddof=1)), + "within_rater_sd_min": float(sds.min()), + "within_rater_sd_median": float(np.median(sds)), + "within_rater_sd_max": float(sds.max()), + } + + +# --------------------------------------------------------------------------- +# Where the coefficient goes as design and rater effects are modelled +# --------------------------------------------------------------------------- +def cumulative_adjustments( + long: pd.DataFrame, cfg: Config, measure: str, panel: str +) -> list[tuple[str, float]]: + """The ladder from the reported figure to a fully specified one.""" + from .build_dataset import rating_matrix + + subset = long[long["panel"] == panel].copy() + matrix = rating_matrix(subset, measure=measure, panel=panel) + estimates = icc_estimates(matrix.to_numpy()) + + rungs = [ + ("as reported: ICC(A,1), severity treated as error", estimates["ICC(A,1)"]), + ("+ severity modelled: ICC(C,1)", estimates["ICC(C,1)"]), + ] + + # Remove each rater's block-specific mean: the block-order drift is a design + # effect the two-way model has no term for, so it lands in the residual. + subset["adjusted"] = ( + subset[measure] + - subset.groupby(["participant", "block"])[measure].transform("mean") + + subset.groupby("participant")[measure].transform("mean") + ) + drift_removed = rating_matrix(subset, measure="adjusted", panel=panel) + rungs.append( + ("+ block-order drift removed within rater", icc_estimates(drift_removed.to_numpy())["ICC(C,1)"]) + ) + + # Drop the weakest raters. + n_drop = int(cfg.analysis.get("drop_weakest_n", 0)) + if n_drop > 0: + ranked = rater_validity(drift_removed) + keep = ranked["rater"].to_numpy()[: len(ranked) - n_drop] + trimmed = drift_removed[keep] + rungs.append( + ( + f"+ {n_drop} lowest-validity raters excluded", + icc_estimates(trimmed.to_numpy())["ICC(C,1)"], + ) + ) + return rungs + + +def block_effects(long: pd.DataFrame, measure: str, panel: str | None = None) -> dict[str, float]: + """Mean drift between the first and second block of 40 descriptions.""" + subset = long if panel is None else long[long["panel"] == panel] + means = subset.groupby("block")[measure].mean() + first, second = float(means.get(1, np.nan)), float(means.get(2, np.nan)) + return {"block1_mean": first, "block2_mean": second, "drift": first - second} + + +# --------------------------------------------------------------------------- +# Restriction of range: splitting by who wrote the description +# --------------------------------------------------------------------------- +def source_split(long: pd.DataFrame, cfg: Config, measure: str, panel: str) -> pd.DataFrame: + """Compare the two description sources on agreement AND on true spread. + + A lower ICC for one source does not mean raters agreed less about it. It can + equally mean that source produced more uniform descriptions, leaving less + true variance to detect. Both quantities are reported side by side. + """ + from .build_dataset import rating_matrix + + rows = [] + for source in (cfg.design["source_human"], cfg.design["source_model"]): + matrix = rating_matrix(long, measure=measure, panel=panel, desc_type=source) + data = matrix.to_numpy() + comp = variance_components(data) + item_means = data.mean(axis=1) + + within = [ + np.abs(row[:, None] - row[None, :])[np.triu_indices(len(row), k=1)].mean() + for row in data + ] + rows.append( + { + "source": source, + "n_items": comp.n_items, + "mean_disagreement_pts": float(np.mean(within)), + # SD across RATERS within an item, averaged over items. Note the + # axis: axis=1 is within-item; axis=0 would be within-rater + # across items, which is a different quantity entirely. + "within_item_sd_pts": float(data.std(axis=1, ddof=1).mean()), + "within_rater_sd_pts": float(data.std(axis=0, ddof=1).mean()), + "var_rater": comp.var_rater, + "var_error": comp.var_error, + "item_mean_sd_pts": float(item_means.std(ddof=1)), + "item_mean_min": float(item_means.min()), + "item_mean_max": float(item_means.max()), + "var_item": comp.var_item, + "ICC(A,1)": icc_estimates(data)["ICC(A,1)"], + "mean_rating": float(data.mean()), + } + ) + return pd.DataFrame(rows) + + +def range_restriction_counterfactual(split: pd.DataFrame, cfg: Config) -> dict[str, float]: + """What the human-written subset's ICC would be with the model subset's spread. + + The counterfactual must be built in the SAME form as the observed value it is + compared against -- here ICC(A,1) -- so the human subset's rater variance stays + in the denominator. Mixing an absolute-agreement observation with a + consistency-form counterfactual silently overstates the gap. + """ + human = split[split["source"] == cfg.design["source_human"]].iloc[0] + model = split[split["source"] == cfg.design["source_model"]].iloc[0] + + # Substitute the model subset's true variance; keep the human subset's + # rater and residual variance exactly as observed. + counterfactual = model["var_item"] / ( + model["var_item"] + human["var_rater"] + human["var_error"] + ) + return { + "observed_human": float(human["ICC(A,1)"]), + "counterfactual_human": float(counterfactual), + # From the ANOVA variance components. + "item_variance_ratio": float(model["var_item"] / human["var_item"]), + # From the observed item means: includes sampling noise, so it differs + # slightly from sqrt(item_variance_ratio). Both are reported because the + # first is the estimand and the second is what a reader can see. + "item_sd_ratio_observed": float( + model["item_mean_sd_pts"] / human["item_mean_sd_pts"] + ), + "item_sd_ratio_components": float(np.sqrt(model["var_item"] / human["var_item"])), + } + + +# --------------------------------------------------------------------------- +# Does the aggregate hold up, and do the panels agree with each other? +# --------------------------------------------------------------------------- +def cross_panel_agreement(long: pd.DataFrame, cfg: Config, measure: str) -> dict[str, float]: + """Correlation between the two panels' item means. + + Individuals within a panel may correlate weakly while the panels themselves + converge -- the signature of real signal recovered from noisy judgements. + """ + expert = long[long["panel"] == cfg.design["panel_expert"]].groupby("item")[measure].mean() + novice = long[long["panel"] == cfg.design["panel_nonexpert"]].groupby("item")[measure].mean() + joined = pd.concat([expert, novice], axis=1, keys=["expert", "nonexpert"]).dropna() + return { + "r": float(joined["expert"].corr(joined["nonexpert"])), + "n_items": int(len(joined)), + } + + +def substantive_effect(long: pd.DataFrame, cfg: Config, measure: str, panel: str) -> dict[str, float]: + """The human-vs-model comparison, with items fixed and then items random. + + Treating items as fixed generalises only to these 80 descriptions. Treating + them as random generalises to new descriptions, which is usually the claim + a paper actually wants to make. + """ + subset = long[long["panel"] == panel] + human, model = cfg.design["source_human"], cfg.design["source_model"] + + wide = subset.pivot_table(index="participant", columns="desc_type", values=measure) + diff = (wide[human] - wide[model]).dropna() + t_participants = stats.ttest_1samp(diff, 0.0) + + item_means = subset.groupby(["item", "desc_type"])[measure].mean().reset_index() + a = item_means.loc[item_means["desc_type"] == human, measure] + b = item_means.loc[item_means["desc_type"] == model, measure] + t_items = stats.ttest_ind(a, b) + pooled_sd = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2) + + return { + "mean_human": float(wide[human].mean()), + "mean_model": float(wide[model].mean()), + "difference": float(diff.mean()), + "n_raters": int(len(diff)), + "t_participants": float(t_participants.statistic), + "df_participants": int(len(diff) - 1), + "p_participants": float(t_participants.pvalue), + "dz": float(diff.mean() / diff.std(ddof=1)), + "t_items": float(t_items.statistic), + "df_items": int(len(a) + len(b) - 2), + "p_items": float(t_items.pvalue), + "d_items": float((a.mean() - b.mean()) / pooled_sd), + } + + +# --------------------------------------------------------------------------- +# The replication question +# --------------------------------------------------------------------------- +def replication_available(long: pd.DataFrame) -> dict[str, int]: + """Did any participant rate any item more than once?""" + per_cell = long.groupby(["participant", "item"]).size() + return { + "min_per_cell": int(per_cell.min()), + "max_per_cell": int(per_cell.max()), + "n_cells": int(len(per_cell)), + "has_replication": bool(per_cell.max() > 1), + } + + +def image_taste_test( + long: pd.DataFrame, cfg: Config, measure: str, panel: str, control_block: bool = True +) -> dict[str, float]: + """Is a rater's deviation from consensus stable across the same image? + + Each rater saw both descriptions of every image, roughly forty trials apart. + If their deviation reflected stable taste for the depicted image, the two + deviations would correlate. + + ``control_block`` matters. The two descriptions of any image necessarily + fall in different blocks, so the block-order drift induces a spurious + NEGATIVE correlation if it is not removed first. Run with + ``control_block=False`` to see the artifact. + """ + subset = long[long["panel"] == panel].copy() + grouping = ["participant", "block"] if control_block else ["participant"] + + subset["residual"] = ( + subset[measure] + - subset.groupby("item")[measure].transform("mean") + - subset.groupby(grouping)[measure].transform("mean") + + subset[measure].mean() + ) + + wide = ( + subset.pivot_table(index=["participant", "image"], columns="desc_type", values="residual") + .dropna() + .reset_index() + ) + human, model = cfg.design["source_human"], cfg.design["source_model"] + per_rater = wide.groupby("participant").apply( + lambda g: g[model].corr(g[human]), include_groups=False + ).dropna() + + test = stats.ttest_1samp(per_rater, 0.0) + ci = stats.t.interval(0.95, len(per_rater) - 1, per_rater.mean(), per_rater.sem()) + return { + "median_r": float(per_rater.median()), + "mean_r": float(per_rater.mean()), + "pooled_r": float(wide[model].corr(wide[human])), + "n_raters": int(len(per_rater)), + "t": float(test.statistic), + "df": int(len(per_rater) - 1), + "p": float(test.pvalue), + "ci_lo": float(ci[0]), + "ci_hi": float(ci[1]), + "block_controlled": control_block, + } + + +def residual_split_consequences(matrix, shares=(0.0, 0.25, 0.50, 0.75)) -> pd.DataFrame: + """What the coefficient becomes under each possible split of the residual. + + The residual confounds stable rater x item interaction with momentary noise. + One rating per cell cannot separate them, and the interpretation of the whole + result turns on the split -- which is what a retest would settle. + """ + comp = variance_components(matrix) + rows = [] + for share in shares: + stable = comp.var_error * share + noise = comp.var_error * (1 - share) + rows.append( + { + "stable_share": share, + "consensus_referenced_ICC": comp.var_item / (comp.var_item + noise), + "interaction_free_ICC": ( + comp.var_item / (comp.var_item + stable) if stable > 0 else float("nan") + ), + } + ) + return pd.DataFrame(rows) + + +def within_rater_measure_coherence(long: pd.DataFrame, panel: str) -> dict[str, float]: + """Correlation between a rater's two sliders, raw and consensus-residualised. + + A high value means each rater is internally self-consistent. Because both + judgements were made on the same screen, it demonstrates stability WITHIN an + encounter, not across encounters -- which is exactly the gap a retest closes. + """ + subset = long[long["panel"] == panel].copy() + if subset["usefulness"].isna().all(): + return {} + + raw = subset.groupby("participant").apply( + lambda g: g["accuracy"].corr(g["usefulness"]), include_groups=False + ) + for measure in ("accuracy", "usefulness"): + subset[f"res_{measure}"] = ( + subset[measure] + - subset.groupby("item")[measure].transform("mean") + - subset.groupby("participant")[measure].transform("mean") + + subset[measure].mean() + ) + residual = subset.groupby("participant").apply( + lambda g: g["res_accuracy"].corr(g["res_usefulness"]), include_groups=False + ) + return { + "raw_median": float(raw.median()), + "raw_q25": float(raw.quantile(0.25)), + "raw_q75": float(raw.quantile(0.75)), + "residual_median": float(residual.median()), + } diff --git a/archivegpt_reliability/archivegpt_reliability/src/icc.py b/archivegpt_reliability/archivegpt_reliability/src/icc.py new file mode 100644 index 0000000..5a1f8fe --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/icc.py @@ -0,0 +1,277 @@ +"""Variance components and the full intraclass-correlation family. + +Implements the two-way random-effects decomposition + + X_ij = mu + t_i + r_j + e_ij + +on a complete items x raters matrix, and every ICC form derived from it. +Confidence intervals follow McGraw & Wong (1996), Table 7. + +Nothing here depends on anything beyond numpy and scipy, so the estimates can +be checked against an independent implementation (see ``crosscheck_pingouin``). + +References +---------- +Shrout, P. E., & Fleiss, J. L. (1979). Psychological Bulletin, 86(2), 420-428. +McGraw, K. O., & Wong, S. P. (1996). Psychological Methods, 1(1), 30-46. +""" + +from __future__ import annotations + +from dataclasses import dataclass, asdict + +import numpy as np +from scipy import stats + + +@dataclass +class Components: + """Variance components and mean squares from the two-way ANOVA.""" + + n_items: int + n_raters: int + grand_mean: float + ms_item: float + ms_rater: float + ms_error: float + ms_within: float # one-way pooled error, for ICC(1,*) + var_item: float # sigma^2_t -- true differences between descriptions + var_rater: float # sigma^2_r -- rater severity (a main effect) + var_error: float # sigma^2_e -- residual: interaction + noise + + @property + def var_total(self) -> float: + return self.var_item + self.var_rater + self.var_error + + @property + def shares(self) -> dict[str, float]: + total = self.var_total + return { + "item": self.var_item / total, + "rater": self.var_rater / total, + "error": self.var_error / total, + } + + def as_dict(self) -> dict: + out = asdict(self) + out["var_total"] = self.var_total + out["shares"] = self.shares + return out + + +def variance_components(matrix) -> Components: + """Two-way random-effects ANOVA on a complete items x raters matrix.""" + data = np.asarray(matrix, dtype=float) + if data.ndim != 2: + raise ValueError("Expected a 2-D items x raters matrix.") + if np.isnan(data).any(): + raise ValueError("Matrix contains missing values; pass a complete matrix.") + n, k = data.shape + if n < 2 or k < 2: + raise ValueError(f"Need at least 2 items and 2 raters, got {n} x {k}.") + + grand = data.mean() + row_means = data.mean(axis=1) + col_means = data.mean(axis=0) + + ss_item = k * ((row_means - grand) ** 2).sum() + ss_rater = n * ((col_means - grand) ** 2).sum() + ss_total = ((data - grand) ** 2).sum() + ss_error = ss_total - ss_item - ss_rater + + ms_item = ss_item / (n - 1) + ms_rater = ss_rater / (k - 1) + ms_error = ss_error / ((n - 1) * (k - 1)) + ms_within = (ss_rater + ss_error) / (n * (k - 1)) + + return Components( + n_items=n, + n_raters=k, + grand_mean=float(grand), + ms_item=float(ms_item), + ms_rater=float(ms_rater), + ms_error=float(ms_error), + ms_within=float(ms_within), + var_item=float((ms_item - ms_error) / k), + var_rater=float((ms_rater - ms_error) / n), + var_error=float(ms_error), + ) + + +def icc_estimates(matrix) -> dict[str, float]: + """All ICC point estimates. + + Keys use the McGraw & Wong names, with Shrout & Fleiss equivalents: + ICC(1,1) = one-way single -- rater and residual inseparable + ICC(A,1) = ICC(2,1) absolute agreement, single rater + ICC(C,1) = ICC(3,1) consistency, single rater + ICC(1,k), ICC(A,k), ICC(C,k) -- the same for the mean of k raters + """ + c = variance_components(matrix) + n, k = c.n_items, c.n_raters + msr, msc, mse, msw = c.ms_item, c.ms_rater, c.ms_error, c.ms_within + + return { + "ICC(1,1)": (msr - msw) / (msr + (k - 1) * msw), + "ICC(A,1)": (msr - mse) / (msr + (k - 1) * mse + k * (msc - mse) / n), + "ICC(C,1)": (msr - mse) / (msr + (k - 1) * mse), + "ICC(1,k)": (msr - msw) / msr, + "ICC(A,k)": (msr - mse) / (msr + (msc - mse) / n), + "ICC(C,k)": (msr - mse) / msr, + } + + +def icc_confidence_intervals(matrix, alpha: float = 0.05) -> dict[str, tuple[float, float]]: + """Confidence intervals for ICC(A,1), ICC(C,1), ICC(A,k), ICC(C,k). + + Follows McGraw & Wong (1996), Table 7. ICC(1,*) is omitted because the + design here is fully crossed and the one-way form should not be reported. + """ + c = variance_components(matrix) + n, k = c.n_items, c.n_raters + msr, msc, mse = c.ms_item, c.ms_rater, c.ms_error + point = icc_estimates(matrix) + q = 1.0 - alpha / 2.0 + + out: dict[str, tuple[float, float]] = {} + + # --- consistency forms: a plain F ratio ------------------------------------ + f_obs = msr / mse + df1, df2 = n - 1, (n - 1) * (k - 1) + f_lower = f_obs / stats.f.ppf(q, df1, df2) + f_upper = f_obs * stats.f.ppf(q, df2, df1) + + out["ICC(C,1)"] = ( + (f_lower - 1) / (f_lower + k - 1), + (f_upper - 1) / (f_upper + k - 1), + ) + out["ICC(C,k)"] = (1 - 1 / f_lower, 1 - 1 / f_upper) + + # --- absolute-agreement forms: Satterthwaite approximation ------------------ + def _satterthwaite(coef_a: float, coef_b: float) -> float: + num = (coef_a * msc + coef_b * mse) ** 2 + den = (coef_a * msc) ** 2 / (k - 1) + (coef_b * mse) ** 2 / ((n - 1) * (k - 1)) + return num / den + + r1 = point["ICC(A,1)"] + if 0 < r1 < 1: + a1 = k * r1 / (n * (1 - r1)) + b1 = 1 + k * r1 * (n - 1) / (n * (1 - r1)) + v1 = _satterthwaite(a1, b1) + fl = stats.f.ppf(q, n - 1, v1) + fu = stats.f.ppf(q, v1, n - 1) + out["ICC(A,1)"] = ( + n * (msr - fl * mse) / (fl * (k * msc + (k * n - k - n) * mse) + n * msr), + n * (fu * msr - mse) / (k * msc + (k * n - k - n) * mse + n * fu * msr), + ) + else: + out["ICC(A,1)"] = (float("nan"), float("nan")) + + rk = point["ICC(A,k)"] + if 0 < rk < 1: + ak = rk / (n * (1 - rk)) + bk = 1 + rk * (n - 1) / (n * (1 - rk)) + vk = _satterthwaite(ak, bk) + fl = stats.f.ppf(q, n - 1, vk) + fu = stats.f.ppf(q, vk, n - 1) + out["ICC(A,k)"] = ( + n * (msr - fl * mse) / (fl * (msc - mse) + n * msr), + n * (fu * msr - mse) / (msc - mse + n * fu * msr), + ) + else: + out["ICC(A,k)"] = (float("nan"), float("nan")) + + return out + + +def single_rater_validity(matrix) -> float: + """Correlation between one rating and the item's true score. + + ICC(C,1) is that correlation SQUARED, so this is its square root. It is the + fair way to describe an individual rater: a median pairwise correlation of + 0.14 between two noisy raters corresponds to each rater correlating 0.39 + with the truth. + """ + c1 = icc_estimates(matrix)["ICC(C,1)"] + return float(np.sqrt(c1)) if c1 > 0 else float("nan") + + +def descriptive_scale(matrix) -> dict[str, float]: + """Translate the components into expected differences in rating points. + + For a normal variable with SD s, the expected absolute difference between + two independent draws is 2s/sqrt(pi). + """ + c = variance_components(matrix) + single_error = np.sqrt(max(c.var_rater, 0) + c.var_error) + return { + "true_difference_between_items": 2 * np.sqrt(max(c.var_item, 0)) / np.sqrt(np.pi), + "difference_between_rater_means": 2 * np.sqrt(max(c.var_rater, 0)) / np.sqrt(np.pi), + "difference_on_same_item": 2 * np.sqrt(c.var_error) / np.sqrt(np.pi), + "band_single_rating": 1.96 * single_error, + "band_panel_mean": 1.96 * single_error / np.sqrt(c.n_raters), + "signal_to_noise": np.sqrt(max(c.var_item, 0)) / single_error, + } + + +def spearman_brown(icc_single: float, k: int) -> float: + """Reliability of the mean of k raters, given the single-rater value.""" + return k * icc_single / (1 + (k - 1) * icc_single) + + +def raters_needed(icc_single: float, target: float) -> float: + """How many raters the mean needs to reach `target` reliability.""" + return target * (1 - icc_single) / (icc_single * (1 - target)) + + +def crosscheck_pingouin(matrix, long_frame=None) -> dict[str, float] | None: + """Re-estimate with pingouin and return the absolute discrepancies. + + Returns None when pingouin is not installed. Any discrepancy larger than + floating-point noise indicates a bug in this module and should be treated + as a failure. + """ + try: + import pandas as pd + import pingouin as pg + except ImportError: + return None + + data = np.asarray(matrix, dtype=float) + n, k = data.shape + tidy = pd.DataFrame( + { + "item": np.repeat(np.arange(n), k), + "rater": np.tile(np.arange(k), n), + "rating": data.ravel(), + } + ) + theirs = pg.intraclass_corr( + data=tidy, targets="item", raters="rater", ratings="rating" + ).set_index("Type")["ICC"] + + # pingouin >= 0.5.3 uses the McGraw & Wong names directly; older releases + # used ICC1/ICC2/ICC3. Accept either. + legacy = { + "ICC1": "ICC(1,1)", + "ICC2": "ICC(A,1)", + "ICC3": "ICC(C,1)", + "ICC1k": "ICC(1,k)", + "ICC2k": "ICC(A,k)", + "ICC3k": "ICC(C,k)", + } + ours = icc_estimates(matrix) + discrepancies: dict[str, float] = {} + for label in ours: + if label in theirs.index: + discrepancies[label] = abs(float(theirs[label]) - ours[label]) + else: + for old, new in legacy.items(): + if new == label and old in theirs.index: + discrepancies[label] = abs(float(theirs[old]) - ours[label]) + if not discrepancies: + raise RuntimeError( + "pingouin returned no recognisable ICC labels " + f"(saw {list(theirs.index)}); cross-check could not be run." + ) + return discrepancies diff --git a/archivegpt_reliability/archivegpt_reliability/src/report.py b/archivegpt_reliability/archivegpt_reliability/src/report.py new file mode 100644 index 0000000..05787ad --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/report.py @@ -0,0 +1,302 @@ +"""Render every table in the technical note as plain text, and as JSON.""" + +from __future__ import annotations + +import json +from datetime import date +from pathlib import Path + +import numpy as np +import pandas as pd + +from . import diagnostics as dx +from .build_dataset import rating_matrix +from .config import Config +from .icc import ( + descriptive_scale, + icc_confidence_intervals, + icc_estimates, + raters_needed, + single_rater_validity, + variance_components, +) + +RULE = "=" * 78 +THIN = "-" * 78 + + +def _h(title: str) -> str: + return f"\n{RULE}\n{title}\n{RULE}" + + +def build_report(long: pd.DataFrame, cfg: Config) -> tuple[str, dict]: + """Run the full analysis and return (text report, machine-readable results).""" + measure = cfg.analysis["primary_measure"] + alpha = float(cfg.analysis["alpha"]) + expert = cfg.design["panel_expert"] + nonexpert = cfg.design["panel_nonexpert"] + + out: list[str] = [] + results: dict = { + "generated": date.today().isoformat(), + "measure": measure, + "data_root": str(cfg.data_root), + } + + out.append(f"ArchiveGPT inter-rater reliability -- {measure}") + out.append(f"generated {date.today().isoformat()} | data: {cfg.data_root}") + + # -- sample --------------------------------------------------------------- + out.append(_h("SAMPLE")) + panels = long.groupby("panel")["participant"].nunique().to_dict() + out.append(f" {len(long):,} ratings, {long['participant'].nunique()} participants, " + f"{long['item'].nunique()} items ({long['image'].nunique()} images x 2 sources)") + out.append(f" panels: {panels}") + if "expertise_detail" in long.columns: + detail = ( + long.drop_duplicates("participant") + .groupby("expertise_detail")["participant"] + .nunique() + .sort_values(ascending=False) + .to_dict() + ) + out.append(f" expertise detail: {detail}") + results["sample"] = {"n_ratings": len(long), "panels": panels} + + # -- Part I: components --------------------------------------------------- + out.append(_h("PART I -- VARIANCE COMPONENTS")) + results["components"] = {} + for panel in (expert, nonexpert): + matrix = rating_matrix(long, measure=measure, panel=panel) + comp = variance_components(matrix.to_numpy()) + shares = comp.shares + out.append(f"\n {panel} ({comp.n_items} items x {comp.n_raters} raters, " + f"grand mean {comp.grand_mean:.1f})") + out.append(f" {'component':<26s}{'variance':>10s}{'SD':>8s}{'share':>9s}") + out.append(f" {THIN[:53]}") + for label, var in ( + ("item (signal)", comp.var_item), + ("rater (severity)", comp.var_rater), + ("residual", comp.var_error), + ): + key = label.split()[0] + out.append( + f" {label:<26s}{var:10.1f}{np.sqrt(max(var,0)):8.1f}" + f"{shares[key if key != 'residual' else 'error']:9.1%}" + ) + out.append(f" {'total (one rating)':<26s}{comp.var_total:10.1f}" + f"{np.sqrt(comp.var_total):8.1f}{1.0:9.1%}") + results["components"][panel] = comp.as_dict() + + # -- Part I: the ICC family ---------------------------------------------- + out.append(_h("PART I -- THE ICC FAMILY")) + results["icc"] = {} + for panel in (expert, nonexpert): + matrix = rating_matrix(long, measure=measure, panel=panel) + data = matrix.to_numpy() + est = icc_estimates(data) + cis = icc_confidence_intervals(data, alpha=alpha) + k = data.shape[1] + out.append(f"\n {panel} (k = {k})") + out.append(f" {'form':<12s}{'estimate':>10s} {int((1-alpha)*100)}% CI") + out.append(f" {THIN[:44]}") + for form in ("ICC(1,1)", "ICC(A,1)", "ICC(C,1)", "ICC(A,k)", "ICC(C,k)"): + ci = cis.get(form) + ci_text = f"[{ci[0]:.3f}, {ci[1]:.3f}]" if ci else "--" + out.append(f" {form:<12s}{est[form]:10.3f} {ci_text}") + validity = single_rater_validity(data) + out.append(f" single-rater validity = sqrt(ICC(C,1)) = {validity:.2f}" + f" (a median pairwise r of {validity**2:.2f} is this value squared)") + results["icc"][panel] = {"estimates": est, "ci": {k_: list(v) for k_, v in cis.items()}} + + # -- Part I: descriptive scale ------------------------------------------- + out.append(_h("PART I -- THE COMPONENTS IN RATING POINTS (expert panel)")) + matrix = rating_matrix(long, measure=measure, panel=expert) + scale = descriptive_scale(matrix.to_numpy()) + out.append(f" true difference between two random descriptions {scale['true_difference_between_items']:6.1f}") + out.append(f" difference between two raters' overall averages {scale['difference_between_rater_means']:6.1f}") + out.append(f" difference between two raters on the SAME item {scale['difference_on_same_item']:6.1f}") + out.append(f" 95% band, one rating around an item's true score +/-{scale['band_single_rating']:5.1f}") + out.append(f" 95% band, panel mean around an item's true score +/-{scale['band_panel_mean']:5.1f}") + out.append(f" signal-to-noise for a single rating {scale['signal_to_noise']:6.2f}") + results["descriptive_scale"] = scale + + single = icc_estimates(matrix.to_numpy())["ICC(A,1)"] + out.append(f"\n raters needed for an aggregate reliability of " + f"0.70 / 0.80 / 0.90: " + f"{raters_needed(single, 0.70):.0f} / {raters_needed(single, 0.80):.0f} / " + f"{raters_needed(single, 0.90):.0f}") + + # -- Part II: pairwise + homogeneity ------------------------------------- + out.append(_h("PART II -- PAIRWISE AGREEMENT AND PANEL HOMOGENEITY")) + results["pairwise"] = {} + for panel in (expert, nonexpert): + data = rating_matrix(long, measure=measure, panel=panel).to_numpy() + pw = dx.pairwise_stats(data) + out.append(f"\n {panel}: {pw['n_pairs']} pairs, median r = {pw['median']:.3f}, " + f"SD = {pw['sd']:.3f}") + out.append(f" below r = 0.10: {pw['share_below_0.10']:.1%} negative: {pw['share_negative']:.1%}") + results["pairwise"][panel] = pw + + data = rating_matrix(long, measure=measure, panel=expert).to_numpy() + hom = dx.homogeneity_test( + data, n_sims=int(cfg.analysis["homogeneity_sims"]), seed=int(cfg.analysis["seed"]) + ) + out.append(f"\n Homogeneity test ({expert} panel, {hom['n_sims']} simulations):") + out.append(f" observed SD of pairwise r {hom['observed_sd']:.3f}") + out.append(f" expected if raters interchangeable {hom['expected_sd']:.3f} " + f"(95% range {hom['expected_lo']:.3f}-{hom['expected_hi']:.3f})") + out.append(f" verdict: {'HETEROGENEOUS -- rater quality genuinely differs' if hom['heterogeneous'] else 'consistent with interchangeable raters'}") + if hom["heterogeneous"]: + out.append(f" implied SD of true pairwise agreement {hom['implied_true_sd']:.3f}") + results["homogeneity"] = hom + + # -- Part II: individual raters ------------------------------------------ + out.append(_h("PART II -- INDIVIDUAL RATER VALIDITY (expert panel)")) + matrix = rating_matrix(long, measure=measure, panel=expert) + val = dx.rater_validity(matrix) + qs = [0, 25, 50, 75, 100] + out.append(f" {'':<10s}" + "".join(f"{f'p{q}':>10s}" for q in qs)) + out.append(f" {'corr rest':<10s}" + "".join( + f"{np.percentile(val['corr_with_rest'], q):+10.2f}" for q in qs)) + out.append(f" {'validity':<10s}" + "".join( + f"{np.percentile(val['implied_validity'], q):+10.2f}" for q in qs)) + n_dead = int((val["corr_with_rest"] <= 0.10).sum()) + n_good = int((val["implied_validity"] >= 0.50).sum()) + out.append(f"\n contributing essentially nothing (r <= 0.10): {n_dead} of {len(val)}") + out.append(f" reaching validity >= 0.50: {n_good} of {len(val)}") + results["rater_validity"] = { + "table": val.to_dict(orient="records"), + "n_noncontributing": n_dead, + "n_solid": n_good, + } + + sev = dx.severity_stats(matrix.to_numpy()) + out.append(f"\n Severity (scale use, independent of ordering):") + out.append(f" rater means span {sev['rater_mean_min']:.1f} to {sev['rater_mean_max']:.1f} " + f"({sev['rater_mean_spread']:.1f} points), SD {sev['rater_mean_sd']:.1f}") + out.append(f" within-rater SD spans {sev['within_rater_sd_min']:.1f} to " + f"{sev['within_rater_sd_max']:.1f}") + results["severity"] = sev + + # -- Part II: cumulative ladder ------------------------------------------ + out.append(_h("PART II -- COEFFICIENT AS DESIGN AND RATER EFFECTS ARE MODELLED")) + ladder = dx.cumulative_adjustments(long, cfg, measure, expert) + for label, value in ladder: + out.append(f" {label:<52s} {value:.3f}") + blocks = dx.block_effects(long, measure, expert) + out.append(f"\n block-order drift: block 1 mean {blocks['block1_mean']:.1f}, " + f"block 2 mean {blocks['block2_mean']:.1f} ({blocks['drift']:+.1f})") + results["cumulative"] = ladder + results["block_effects"] = blocks + + # -- Part II: source split ----------------------------------------------- + out.append(_h("PART II -- SPLIT BY DESCRIPTION SOURCE (expert panel)")) + split = dx.source_split(long, cfg, measure, expert) + human, model = cfg.design["source_human"], cfg.design["source_model"] + rows = [ + ("mean disagreement between two raters (pts)", "mean_disagreement_pts", "{:.1f}"), + ("within-item SD, across raters (pts)", "within_item_sd_pts", "{:.1f}"), + ("residual variance", "var_error", "{:.1f}"), + ("SD of true item means (pts)", "item_mean_sd_pts", "{:.1f}"), + ("item variance", "var_item", "{:.1f}"), + ("ICC(A,1)", "ICC(A,1)", "{:.3f}"), + ("mean rating", "mean_rating", "{:.1f}"), + ] + out.append(f" {'':<44s}{human:>14s}{model:>14s}") + out.append(f" {THIN[:72]}") + for label, key, fmt in rows: + a = split.loc[split["source"] == human, key].iloc[0] + b = split.loc[split["source"] == model, key].iloc[0] + out.append(f" {label:<44s}{fmt.format(a):>14s}{fmt.format(b):>14s}") + + cf = dx.range_restriction_counterfactual(split, cfg) + out.append(f"\n RANGE RESTRICTION CHECK") + out.append(f" A lower ICC for one source does not mean raters agreed less about it.") + out.append(f" observed ICC(A,1) for {human}-written: {cf['observed_human']:.3f}") + out.append(f" same raters, with the {model} subset's true spread: {cf['counterfactual_human']:.3f}") + out.append(f" true variance ratio {model}/{human}: {cf['item_variance_ratio']:.2f}x " + f"({cf['item_sd_ratio_observed']:.2f}x on the observed item-mean SD scale)") + results["source_split"] = split.to_dict(orient="records") + results["range_restriction"] = cf + + # -- Part II: does the aggregate hold ------------------------------------ + out.append(_h("PART II -- THE AGGREGATE, AND THE SUBSTANTIVE EFFECT")) + for panel in (expert, nonexpert): + data = rating_matrix(long, measure=measure, panel=panel).to_numpy() + est = icc_estimates(data) + out.append(f" {panel:<12s} single rater {est['ICC(A,1)']:.3f} " + f"mean of {data.shape[1]:>3d} raters {est['ICC(A,k)']:.3f}") + xp = dx.cross_panel_agreement(long, cfg, measure) + out.append(f"\n the two panels rank the {xp['n_items']} descriptions at r = {xp['r']:.3f}") + results["cross_panel"] = xp + + results["substantive"] = {} + for panel in (expert, nonexpert): + eff = dx.substantive_effect(long, cfg, measure, panel) + out.append(f"\n {panel}: {human}-written {eff['mean_human']:.1f} vs " + f"{model}-written {eff['mean_model']:.1f} ({eff['difference']:+.2f} pts)") + out.append(f" items FIXED : t({eff['df_participants']}) = {eff['t_participants']:.2f}, " + f"p = {eff['p_participants']:.2g}, dz = {eff['dz']:.2f}") + out.append(f" items RANDOM : t({eff['df_items']}) = {eff['t_items']:.2f}, " + f"p = {eff['p_items']:.2g}, d = {eff['d_items']:.2f}") + results["substantive"][panel] = eff + + # -- Part II: replication ------------------------------------------------- + out.append(_h("PART II -- REPLICATION AND THE RETEST QUESTION")) + rep = dx.replication_available(long) + out.append(f" ratings per participant x item: min {rep['min_per_cell']}, " + f"max {rep['max_per_cell']} across {rep['n_cells']:,} cells") + out.append(f" replication available: {rep['has_replication']}") + results["replication"] = rep + + out.append(f"\n Proxy test -- is a rater's deviation stable across the same IMAGE?") + for controlled in (False, True): + taste = dx.image_taste_test(long, cfg, measure, expert, control_block=controlled) + tag = "block drift controlled" if controlled else "block drift NOT controlled (artifact)" + out.append(f" {tag:<38s} median r = {taste['median_r']:+.3f}, " + f"t({taste['df']}) = {taste['t']:+.2f}, p = {taste['p']:.3f}") + if controlled: + out.append(f" {'':<38s} 95% CI [{taste['ci_lo']:+.3f}, {taste['ci_hi']:+.3f}]") + results["image_taste"] = taste + + coherence = dx.within_rater_measure_coherence(long, expert) + if coherence: + out.append(f"\n Within-rater coherence across the two sliders:") + out.append(f" raw median r = {coherence['raw_median']:.2f}, " + f"consensus-residualised median r = {coherence['residual_median']:.2f}") + out.append(f" (same screen seconds apart: stability WITHIN an encounter, " + f"not across encounters)") + results["measure_coherence"] = coherence + + out.append(f"\n Consequences of the unidentified split of the residual:") + consequences = dx.residual_split_consequences( + rating_matrix(long, measure=measure, panel=expert).to_numpy() + ) + out.append(f" {'stable share':>14s}{'consensus ICC':>16s}{'interaction-free ICC':>22s}") + for _, row in consequences.iterrows(): + free = row["interaction_free_ICC"] + out.append(f" {row['stable_share']:>13.0%}{row['consensus_referenced_ICC']:>16.3f}" + f"{(f'{free:.3f}' if np.isfinite(free) else '--'):>22s}") + results["residual_split"] = consequences.to_dict(orient="records") + + return "\n".join(out), results + + +def write_outputs(text: str, results: dict, long: pd.DataFrame, cfg: Config) -> list[Path]: + """Write the text report, the JSON results and the rebuilt dataset.""" + written = [] + + report_path = cfg.output_dir / "reliability_report.txt" + report_path.write_text(text, encoding="utf-8") + written.append(report_path) + + json_path = cfg.output_dir / "results.json" + json_path.write_text(json.dumps(results, indent=2, default=float), encoding="utf-8") + written.append(json_path) + + data_path = cfg.output_dir / "ratings_long.csv" + long.to_csv(data_path, index=False) + written.append(data_path) + + return written diff --git a/archivegpt_reliability/archivegpt_reliability/src/validate.py b/archivegpt_reliability/archivegpt_reliability/src/validate.py new file mode 100644 index 0000000..85b43fd --- /dev/null +++ b/archivegpt_reliability/archivegpt_reliability/src/validate.py @@ -0,0 +1,186 @@ +"""Prove the rebuilt item identity is correct by checking it against the release. + +The reliability analysis stands or falls on whether ``image`` + ``desc_type`` +was reassembled correctly from the raw exports. Two independent checks are +available, because the published clean data -- while it cannot identify items -- +does contain the same ratings in a different arrangement: + +1. Every rating must appear in ``correctness_per_trial.csv``. Comparison is on + each participant's sorted values, since the published ``trial`` column is + serial position and carries no item identity to align on. +2. The implied participant x desc_type means must reproduce + ``accuracy_usefulness.csv`` exactly. + +If both pass, the reconstruction contains precisely the published ratings, and +the only thing added is the item label. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field + +import numpy as np +import pandas as pd + +from .config import Config + + +class ValidationError(RuntimeError): + """Raised when the reconstruction disagrees with the published data.""" + + +@dataclass +class ValidationReport: + checks: list[tuple[str, bool, str]] = field(default_factory=list) + + def add(self, name: str, passed: bool, detail: str) -> None: + self.checks.append((name, passed, detail)) + + @property + def ok(self) -> bool: + return all(passed for _, passed, _ in self.checks) + + def render(self) -> str: + lines = [] + for name, passed, detail in self.checks: + mark = "PASS" if passed else "FAIL" + lines.append(f" [{mark}] {name}: {detail}") + return "\n".join(lines) + + +def validate_reconstruction(long: pd.DataFrame, cfg: Config) -> ValidationReport: + """Run every available check and return the report.""" + report = ValidationReport() + tolerance = float(cfg.validation["tolerance"]) + measure = cfg.analysis["primary_measure"] + + _check_design(long, cfg, report) + _check_trial_file(long, cfg, report, tolerance, measure) + _check_aggregate_file(long, cfg, report, tolerance) + + if cfg.validation.get("enforce", True) and not report.ok: + raise ValidationError( + "Reconstruction does not match the published clean data:\n" + + report.render() + + "\n\nSet validation.enforce: false in config.yaml to continue anyway." + ) + return report + + +def _check_design(long: pd.DataFrame, cfg: Config, report: ValidationReport) -> None: + design = cfg.design + n_part = long["participant"].nunique() + n_items = long["item"].nunique() + expected_rows = n_part * design["n_items_expected"] + + report.add( + "fully crossed design", + len(long) == expected_rows and long.groupby(["participant", "item"]).size().max() == 1, + f"{len(long):,} ratings = {n_part} participants x {n_items} items, " + f"one rating per cell", + ) + + per_block = long.groupby(["participant", "block"])["item"].nunique() + report.add( + "two balanced blocks", + set(per_block.unique()) <= {design["n_items_expected"] // 2}, + f"each participant rated {design['n_items_expected'] // 2} items per block", + ) + + flip = ( + long.pivot_table( + index=["participant", "image"], columns="block", values="desc_type", aggfunc="first" + ) + .dropna() + ) + if flip.shape[1] == 2: + same = (flip.iloc[:, 0] == flip.iloc[:, 1]).sum() + report.add( + "counterbalancing exact", + same == 0, + f"description source flips between blocks for every image " + f"({same} exceptions)", + ) + + +def _check_trial_file( + long: pd.DataFrame, + cfg: Config, + report: ValidationReport, + tolerance: float, + measure: str, +) -> None: + """Compare per-participant sorted ratings with the published trial file.""" + if not cfg.trial_file.is_file(): + report.add("published trial file", True, f"not present, skipped ({cfg.trial_file.name})") + return + + published = pd.read_csv(cfg.trial_file) + published_col = {"accuracy": "perc_acc", "usefulness": "perc_use"}.get(measure) + if published_col not in published.columns: + report.add( + "published trial file", + True, + f"no '{published_col}' column, skipped", + ) + return + + worst = 0.0 + mismatched = 0 + for participant, ours in long.groupby("participant")[measure]: + theirs = published.loc[published["Participant"] == participant, published_col].dropna() + if len(theirs) != len(ours): + mismatched += 1 + continue + diff = np.abs(np.sort(ours.to_numpy()) - np.sort(theirs.to_numpy())).max() + worst = max(worst, float(diff)) + + report.add( + f"ratings match {cfg.trial_file.name}", + worst <= tolerance and mismatched == 0, + f"max abs difference {worst:.2e} across " + f"{long['participant'].nunique()} participants, {mismatched} length mismatches", + ) + + +def _check_aggregate_file( + long: pd.DataFrame, cfg: Config, report: ValidationReport, tolerance: float +) -> None: + """Compare implied participant x desc_type means with the published aggregate.""" + if not cfg.aggregate_file.is_file(): + report.add( + "published aggregate file", True, f"not present, skipped ({cfg.aggregate_file.name})" + ) + return + + published = pd.read_csv(cfg.aggregate_file) + pairs = [("accuracy", "mean_perc_acc"), ("usefulness", "mean_perc_use")] + worst = 0.0 + compared = 0 + + for measure, published_col in pairs: + if published_col not in published.columns or long[measure].isna().all(): + continue + ours = ( + long.groupby(["participant", "desc_type"])[measure] + .mean() + .rename("ours") + .reset_index() + ) + merged = ours.merge( + published[["Participant", "desc_type", published_col]], + left_on=["participant", "desc_type"], + right_on=["Participant", "desc_type"], + how="inner", + ) + if merged.empty: + continue + diff = (merged["ours"] - merged[published_col]).abs().max() + worst = max(worst, float(diff)) + compared += len(merged) + + report.add( + f"means match {cfg.aggregate_file.name}", + worst <= tolerance, + f"max abs difference {worst:.2e} across {compared} participant x source cells", + )