initial commit

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.

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

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

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

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())

archivegpt_reliability/archivegpt_reliability/src/__init__.py

@@ -0,0 +1,2 @@
+"""ArchiveGPT inter-rater reliability analysis."""
+__version__ = "1.0.0"

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()

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,
+    )

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()),
+    }

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

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

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",
+    )