analysis(tdcs): GLM hypothesis tests + methods for the reaching study

Poisson GEE / binomial rate / learning-rate models, subject random-intercept
sensitivity, unknown-group classification and merge scenarios, plus METHODS.md
and learning-curve plots. Data reconstructed from the experiment DB and
verified against the exported matrix.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

analysis/METHODS.md

@@ -0,0 +1,158 @@
+# tDCS reaching study — GLM methods
+
+This document explains the statistics in `tdcs_glm.py`: the models, the exact
+formulas, how each subject's progress is accounted for, and the caveats.
+
+## The question
+
+Two **known** conditions anchor the performance scale and differ from each other:
+
+- **H2 (anchor check):** `Electrode-Box-B2` performs **better** than `Electrode-Box-A2`.
+
+Two **unknown** conditions are then **classified** against those anchors — for
+each of `Electrode-Box-A` and `Right-Electrode`, is it **A2-like** or **B2-like**?
+We do *not* assume either belongs to B2; each unknown is compared to *both*
+anchors, and the anchor it cannot be distinguished from is its likely class.
+
+There is also a **merged-assumption** mode (`--merge`): assume the two unknowns
+resolve as `Right-Electrode == Box-B2` and `Electrode-Box-A == Box-A2`, fold them
+into the anchors, and re-estimate everything with 4 subjects per group.
+
+## Data and outcome
+
+- **Outcome:** `success` = successful reaches in a session
+  (`analysis_summary.counts.Success`), a non-negative **count**.
+- **Attempts:** `total` = reach attempts in the session (`analysis_summary.total`);
+  used as the denominator for the rate model. `success_rate == success / total`.
+- **Time:** `day` = the "# Days Reach" field (training day, 0..26). Analyses use
+  days ≥ 0 (pre-training negative days and zero-attempt sessions are excluded;
+  zero-attempt sessions are undefined for the rate model).
+- **Groups (subjects):** Naive (4), Box-A2 (3), Box-B2 (3), Box-A (1),
+  Right-Electrode (1). Under `--merge`: Box-A2 (4), Box-B2 (4), Naive (4).
+
+**Provenance:** `tdcs_reach_data.csv` was reconstructed from the experiment
+database and verified cell-by-cell against the exported matrix (71/71
+unambiguous cells on days 0–5 matched exactly).
+
+## The three models
+
+All models use `Electrode-Box-B2` as the **reference** group, so each group term
+is that group's contrast *versus Box-B2*. `day_c` is the centered training day and
+`day_c2 = day_c²` captures the rise-then-plateau of the learning curve.
+
+### (A) Level — count (primary)
+
+Poisson GEE on the success counts, clustered by subject:
+
+```
+success ~ C(group, Treatment('Electrode-Box-B2')) + day_c + day_c2
+family        = Poisson (log link)
+groups        = subject           # repeated-measures cluster
+cov_struct    = Exchangeable      # working within-subject correlation
+SE            = robust (sandwich)
+```
+
+`exp(coef)` for a group term is an **incidence-rate ratio (IRR)**: expected
+successes relative to Box-B2.
+
+### (B) Level — rate
+
+Binomial GLM on successes-out-of-attempts, with cluster-robust SEs by subject:
+
+```
+cbind(success, total - success) ~ C(group, Treatment('Electrode-Box-B2')) + day_c + day_c2
+family   = Binomial (logit link)
+cov_type = cluster (groups = subject)
+```
+
+`exp(coef)` is an **odds ratio** for a successful reach relative to Box-B2. This
+controls for differing numbers of attempts, so it answers "who is more accurate
+per attempt?" rather than "who attempts more?"
+
+### (C) Learning rate
+
+Poisson GEE with a **group × day** interaction, to ask whether groups improve at
+different *rates* (not just different levels):
+
+```
+success ~ C(group, Treatment('Electrode-Box-B2')) * day_c + day_c2
+```
+
+Each `group[T.X]:day_c` term is the difference in log-slope versus Box-B2; a joint
+Wald test asks whether *any* group's slope differs. Large p ⇒ parallel learning.
+
+### Anchors-only model
+
+The A2-vs-B2 comparison, refit on **just the two anchor groups** so nothing else
+influences the shared day terms or the dispersion/correlation nuisance:
+
+```
+# subset to {Box-A2, Box-B2}
+success ~ C(group, Treatment('Electrode-Box-B2')) + day_c + day_c2     # count
+cbind(success, total-success) ~ C(group, ...) + day_c + day_c2         # rate
+```
+
+With B2 as reference the single group term is the A2-vs-B2 effect; H2 is a
+one-sided test that this coefficient is below zero.
+
+## How each subject's progress is accounted for
+
+Two distinct pieces:
+
+1. **Progress over training** — the `day_c + day_c2` fixed terms model the average
+   learning curve, so groups are compared at comparable points in training rather
+   than being confounded by *when* each was measured.
+
+2. **Repeated measures / individual baselines** — each subject contributes many
+   correlated sessions and has its own baseline. Handled two ways:
+
+   - **Primary (GEE):** subject is the cluster; an exchangeable working
+     correlation plus robust (sandwich) SEs give *population-average* group
+     effects whose inference is valid under within-subject correlation and
+     Poisson overdispersion.
+
+   - **Sensitivity (mixed model):** a Poisson model with a **per-subject random
+     intercept** (`(1 | subject)`) so each animal gets its own baseline level;
+     the group effects are estimated after allowing for that individual variation.
+     statsmodels has no frequentist Poisson GLMM, so this is the **MAP/Laplace**
+     fit (`fit_vb` diverges on these large counts). Its p-values are approximate —
+     read it as a direction/magnitude check that should agree with GEE.
+
+## Classification logic (unknowns)
+
+For each unknown group, `diff_contrast` builds a linear contrast of that group
+against **each anchor** (both coded vs the B2 reference) and tests it:
+
+- indistinguishable from B2 **and** different from A2 → **B2-like**
+- indistinguishable from A2 **and** different from B2 → **A2-like**
+- indistinguishable from both → **ambiguous** (report the numerically nearer one)
+- different from both → **unlike both**
+
+## Caveats
+
+- **Tiny groups.** Box-A2/B2 have 3 subjects, Naive 4, and each unknown has
+  **n = 1 subject**. Classifying a one-subject condition is weak: "matches
+  anchor X" means "not statistically distinguishable from X," **not** proof of
+  equivalence.
+- **Single-cluster fragility.** Cluster-robust/GEE inference with one cluster in a
+  group can produce artificially small SEs — most visibly the learning-rate
+  interaction for the n=1 groups; do not read those p-values literally.
+- **Mixed-model confounding.** With a per-subject random intercept, a group made
+  of one subject is partly confounded with that subject's random intercept, so its
+  fixed effect is shrunk.
+- **Rate vs count.** "Success" alone is a count; the rate model (success/attempts)
+  is the fairer accuracy comparison when attempt counts differ.
+
+## Files and usage
+
+- `tdcs_glm.py` — the analysis (run it directly).
+- `tdcs_reach_data.csv` — the verified long-format data (`subject, group, day, success, total`).
+- `tdcs_learning_curves*.png` — per-group learning curves (count and rate panels).
+
+```
+python3 analysis/tdcs_glm.py               # all training days, 5 groups
+python3 analysis/tdcs_glm.py --max-day 10  # restrict to days 0–10
+python3 analysis/tdcs_glm.py --merge       # assume Right==B2 and Box-A==A2 (3 groups)
+```
+
+Requires: pandas, numpy, scipy, statsmodels, matplotlib.

analysis/tdcs_glm.py

@@ -0,0 +1,440 @@
+#!/usr/bin/env python3
+"""
+tDCS reaching study — GLM hypothesis tests.
+
+QUESTION
+    Two KNOWN conditions anchor the scale and differ from each other:
+        H2 (anchor check): Electrode-Box-B2 performs BETTER than Electrode-Box-A2.
+    Two UNKNOWN conditions then get CLASSIFIED against those anchors:
+        For "Electrode-Box-A" and for "Right-Electrode", is each one
+        A2-like or B2-like? (We do not assume either belongs to B2.)
+    Each unknown is compared to BOTH anchors; a group it cannot be
+    distinguished from is its likely class.
+
+Each comparison is examined through THREE complementary GLMs:
+
+    (A) LEVEL / count      Poisson GEE on `success` (successful reaches per
+                           session), adjusting for training day. Answers
+                           "who makes more successful reaches overall?"
+
+    (B) LEVEL / rate       Binomial GLM on success / attempts
+                           (`success` out of `total`), cluster-robust.
+                           Controls for differing numbers of attempts —
+                           "who is more accurate per attempt?"
+
+    (C) LEARNING RATE      Poisson GEE with a group x day interaction.
+                           Answers "do the groups improve at different
+                           RATES?" (slope of the learning curve).
+
+All models adjust for time with day (and day^2 for the plateau) and account
+for repeated measures on each subject via GEE (subject = cluster,
+exchangeable correlation, robust SE) or cluster-robust standard errors.
+Group is coded with Electrode-Box-B2 as the REFERENCE, so H1 and H2 read
+directly off the group terms:
+    Right-Electrode term  -> H1 (expect ~0 / non-significant)
+    Electrode-Box-A2 term -> H2 (expect < 0: B2 above A2)
+exp(coef) is an incidence-rate ratio (count model) or odds ratio (rate
+model) relative to Box-B2.
+
+CAVEATS
+    * Tiny groups: Box-A2 n=3, Box-B2 n=3, Naive n=4, and each UNKNOWN
+      (Electrode-Box-A, Right-Electrode) has only n=1 SUBJECT. Classifying a
+      1-subject condition is weak: "matches anchor X" means "not statistically
+      distinguishable from X", NOT proof of equivalence, and cluster-robust
+      inference with a single cluster in a group is fragile (can show
+      artificially small SEs, especially in the learning-rate interaction).
+
+DATA PROVENANCE
+    analysis/tdcs_reach_data.csv was reconstructed from the experiment
+    database (experiment "tDCS", outcome = analysis_summary.counts.Success,
+    attempts = analysis_summary.total, x-axis = the "# Days Reach" field) and
+    verified cell-by-cell against the exported matrix (71/71 unambiguous
+    cells on days 0-5 matched exactly).
+
+USAGE
+    python3 analysis/tdcs_glm.py
+    Requires: pandas, numpy, scipy, statsmodels, matplotlib.
+    Writes: analysis/tdcs_learning_curves.png
+"""
+import os
+import argparse
+import numpy as np
+import pandas as pd
+import patsy
+import statsmodels.api as sm
+import statsmodels.formula.api as smf
+from scipy import stats
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+
+REF = "Electrode-Box-B2"
+# Two KNOWN anchor conditions that differ (Box-B2 > Box-A2, see H2).
+ANCHOR_LOW = "Electrode-Box-A2"
+ANCHOR_HIGH = "Electrode-Box-B2"
+# Two UNKNOWN conditions to classify: is each one A2-like or B2-like?
+CANON_UNKNOWNS = ["Electrode-Box-A", "Right-Electrode"]
+ALL_GROUPS = ["Naive", ANCHOR_LOW, ANCHOR_HIGH] + CANON_UNKNOWNS
+# The active group config is set at runtime (after any --merge/--assign remap).
+UNKNOWNS = list(CANON_UNKNOWNS)
+MAIN = list(ALL_GROUPS)
+OTHERS = [g for g in MAIN if g != REF]
+# --merge preset: assume each unknown resolves into an anchor.
+MERGE_MAP = {"Electrode-Box-A": ANCHOR_LOW, "Right-Electrode": ANCHOR_HIGH}
+HERE = os.path.dirname(os.path.abspath(__file__))
+DATA = os.path.join(HERE, "tdcs_reach_data.csv")
+PLOT = os.path.join(HERE, "tdcs_learning_curves.png")
+RHS = f"C(group, Treatment('{REF}')) + day_c + day_c2"
+
+
+def gterm(name):
+    """statsmodels/patsy column name for a group level (vs the reference)."""
+    return f"C(group, Treatment('{REF}'))[T.{name}]"
+
+
+def xterm(name):
+    """Interaction column name (group level x day slope)."""
+    return f"C(group, Treatment('{REF}'))[T.{name}]:day_c"
+
+
+def load(max_day=None, mapping=None):
+    df = pd.read_csv(DATA)
+    if mapping:                                      # remap unknown -> target group
+        df["group"] = df["group"].map(lambda g: mapping.get(g, g))
+    df = df[df["group"].isin(ALL_GROUPS)].copy()    # keep canonical groups
+    df = df[df["day"] >= 0].copy()                   # training days only (day 0..)
+    if max_day is not None:
+        df = df[df["day"] <= max_day].copy()         # restrict the analysis window
+    df["day_c"] = df["day"] - df["day"].mean()       # center day
+    df["day_c2"] = df["day_c"] ** 2
+    present = [g for g in ALL_GROUPS if g != REF and g in set(df["group"])]
+    df["group"] = pd.Categorical(df["group"], categories=[REF] + present)
+    return df
+
+
+def contrast(res, tname):
+    """coef, robust SE, p, effect=exp(coef) and its 95% CI for one term."""
+    ci = res.conf_int()
+    lo, hi = ci.loc[tname]
+    return dict(coef=res.params[tname], se=res.bse[tname], p=res.pvalues[tname],
+                eff=np.exp(res.params[tname]), lo=np.exp(lo), hi=np.exp(hi))
+
+
+def one_sided_p_below(c):
+    """One-sided p for H: coef < 0 (i.e. Box-B2 above the compared group)."""
+    return stats.norm.cdf(c["coef"] / c["se"])
+
+
+def diff_contrast(res, a, b):
+    """Compare group `a` vs group `b` (both coded relative to REF) via a linear
+    contrast. Returns coef, robust SE, p, effect=exp(coef)=ratio a/b, and CI.
+    Works for any pair, including b == REF (then it is just the `a` term)."""
+    names = list(res.params.index)
+    v = np.zeros(len(names))
+    if a != REF:
+        v[names.index(gterm(a))] += 1.0
+    if b != REF:
+        v[names.index(gterm(b))] -= 1.0
+    tt = res.t_test(v)
+    coef = float(np.ravel(tt.effect)[0])
+    se = float(np.ravel(tt.sd)[0])
+    p = float(np.ravel(tt.pvalue)[0])
+    lo, hi = np.ravel(tt.conf_int())[:2]
+    return dict(coef=coef, se=se, p=p, eff=np.exp(coef),
+                lo=np.exp(lo), hi=np.exp(hi))
+
+
+# ----------------------------------------------------------------------------- models
+def fit_count_level(df):
+    """(A) Poisson GEE on success counts; group main effects vs Box-B2."""
+    return smf.gee(f"success ~ {RHS}", groups="subject", data=df,
+                   family=sm.families.Poisson(),
+                   cov_struct=sm.cov_struct.Exchangeable()).fit()
+
+
+def fit_rate_level(df):
+    """(B) Binomial GLM on success/attempts with cluster-robust SE.
+    Zero-attempt sessions (total==0) carry no rate information and are dropped."""
+    d = df[df["total"] > 0].copy()
+    dropped = len(df) - len(d)
+    if dropped:
+        print(f"[rate model] dropped {dropped} zero-attempt session(s) (undefined rate)")
+    X = patsy.dmatrix(RHS, d, return_type="dataframe")
+    endog = np.column_stack([d["success"].values,
+                             (d["total"] - d["success"]).values])
+    return sm.GLM(endog, X, family=sm.families.Binomial()).fit(
+        cov_type="cluster", cov_kwds={"groups": d["subject"].values})
+
+
+def fit_learning_rate(df):
+    """(C) Poisson GEE with group x day interaction (slope differences)."""
+    return smf.gee(
+        f"success ~ C(group, Treatment('{REF}')) * day_c + day_c2",
+        groups="subject", data=df, family=sm.families.Poisson(),
+        cov_struct=sm.cov_struct.Exchangeable()).fit()
+
+
+def fit_mixed_count(df):
+    """Subject random-intercept Poisson mixed model (MAP / Laplace).
+
+    Adds a per-subject random intercept so each subject has its own baseline
+    level; the group fixed effects are then estimated after allowing for that
+    individual variation. A subject-specific complement to the population-
+    average GEE. statsmodels has no frequentist Poisson GLMM, so this is the
+    posterior-mode (MAP) fit with a Laplace covariance for the fixed effects
+    (fit_vb failed to converge on these large counts; MAP is stable)."""
+    import warnings
+    from statsmodels.genmod.bayes_mixed_glm import PoissonBayesMixedGLM
+    vc = {"subject": "0 + C(subject)"}          # random intercept per subject
+    model = PoissonBayesMixedGLM.from_formula(f"success ~ {RHS}", vc, df)
+    with warnings.catch_warnings():
+        # MAP stops at |gradient|~5e-5 (effectively converged); silence the
+        # over-strict "did not converge" notice.
+        warnings.filterwarnings("ignore", message="Laplace fitting did not converge")
+        return model.fit_map()
+
+
+def mixed_contrast(res, tname):
+    """coef, SD, approx p and IRR + 95% interval from a BayesMixedGLM fit."""
+    names = list(res.model.exog_names)
+    i = names.index(tname)
+    mean, sd = float(res.fe_mean[i]), float(res.fe_sd[i])
+    p = 2.0 * stats.norm.sf(abs(mean / sd))
+    return dict(coef=mean, se=sd, p=p, eff=np.exp(mean),
+                lo=np.exp(mean - 1.96 * sd), hi=np.exp(mean + 1.96 * sd))
+
+
+# ----------------------------------------------------------------------------- reports
+def describe(df):
+    print("=" * 78)
+    print("DESCRIPTIVES")
+    print("=" * 78)
+    grp = df.groupby("group", observed=True)
+    tbl = pd.DataFrame({
+        "n_subj": grp["subject"].nunique(),
+        "n_sessions": grp.size(),
+        "mean_success": grp["success"].mean().round(1),
+        "mean_rate": (grp["success"].sum() / grp["total"].sum()).round(3),
+        "max_day": grp["day"].max(),
+    })
+    print(tbl.to_string())
+
+
+def report_level(title, res, effect_label):
+    print("\n" + "=" * 78)
+    print(title)
+    print("=" * 78)
+    print(f"Effect vs {REF} ({effect_label} relative to Box-B2):")
+    h = f"{'group':20s} {'effect':>7s} {'95% CI':>16s} {'coef':>8s} {'SE':>7s} {'p':>9s}"
+    print(h)
+    print("-" * len(h))
+    out = {}
+    for g in OTHERS:
+        c = contrast(res, gterm(g))
+        out[g] = c
+        ci = f"[{c['lo']:.2f}, {c['hi']:.2f}]"
+        print(f"{g:20s} {c['eff']:7.3f} {ci:>16s} {c['coef']:8.3f} {c['se']:7.3f} {c['p']:9.4f}")
+    return out
+
+
+def report_learning(res):
+    print("\n" + "=" * 78)
+    print("(C) LEARNING RATE — Poisson GEE, group x day interaction")
+    print("=" * 78)
+    base = res.params["day_c"]
+    print(f"Box-B2 learning slope: {np.exp(base):.3f}x successes per training day "
+          f"(baseline).")
+    print("Slope DIFFERENCE vs Box-B2 (exp(coef) = per-day multiplier on the rate ratio):")
+    h = f"{'group':20s} {'slopeΔ/day':>10s} {'coef':>8s} {'SE':>7s} {'p':>9s}"
+    print(h)
+    print("-" * len(h))
+    for g in OTHERS:
+        c = contrast(res, xterm(g))
+        print(f"{g:20s} {c['eff']:10.3f} {c['coef']:8.3f} {c['se']:7.3f} {c['p']:9.4f}")
+    # Joint Wald test: do ANY groups differ from Box-B2 in learning rate?
+    names = list(res.params.index)
+    idx = [names.index(xterm(g)) for g in OTHERS]
+    R = np.zeros((len(idx), len(names)))
+    for i, j in enumerate(idx):
+        R[i, j] = 1.0
+    jt = res.wald_test(R, scalar=True)
+    print(f"\nJoint test (all group x day interactions = 0): "
+          f"chi2={float(jt.statistic):.2f}, df={len(idx)}, p={float(jt.pvalue):.4f}")
+    print("  -> small p = groups improve at DIFFERENT rates; large p = parallel learning.")
+
+
+def report_anchor_only(df):
+    """A2-vs-B2 ONLY: refit on just the two anchors so nothing else influences
+    the shared day terms / dispersion. B2 is the reference, so the single group
+    term is the A2-vs-B2 comparison."""
+    print("\n" + "=" * 78)
+    print("ANCHORS ONLY — Box-A2 vs Box-B2 (two-group models)")
+    print("=" * 78)
+    d = df[df["group"].isin([ANCHOR_LOW, ANCHOR_HIGH])].copy()
+    d["day_c"] = d["day"] - d["day"].mean()
+    d["day_c2"] = d["day_c"] ** 2
+    d["group"] = pd.Categorical(d["group"], categories=[ANCHOR_HIGH, ANCHOR_LOW])
+    n_sub = d.groupby("group", observed=True)["subject"].nunique().to_dict()
+    print(f"subjects: Box-B2={n_sub.get(ANCHOR_HIGH)}, Box-A2={n_sub.get(ANCHOR_LOW)}"
+          f"   sessions: {len(d)}")
+    print(f"formula: success ~ C(group, Treatment('{REF}')) + day_c + day_c2")
+    for label, res in [("count/level (Poisson GEE)", fit_count_level(d)),
+                       ("rate/level  (Binomial cluster-robust)", fit_rate_level(d))]:
+        c = contrast(res, gterm(ANCHOR_LOW))
+        p1 = one_sided_p_below(c)
+        tag = "SUPPORTED" if (c["coef"] < 0 and p1 < 0.05) else "not supported"
+        print(f"  [{label}] A2/B2 = {c['eff']:.2f} [{c['lo']:.2f}, {c['hi']:.2f}]  "
+              f"=> B2 = {1.0 / c['eff']:.2f}x A2, one-sided p={p1:.4f}  -> H2 {tag}")
+
+
+def report_mixed(df):
+    """Subject random-intercept Poisson mixed model (sensitivity vs GEE)."""
+    print("\n" + "=" * 78)
+    print("SENSITIVITY — Poisson MIXED model, per-subject random intercept (MAP/Laplace)")
+    print("=" * 78)
+    try:
+        res = fit_mixed_count(df)
+    except Exception as exc:                       # VB can be finicky; degrade gracefully
+        print(f"  (mixed model skipped: {type(exc).__name__}: {exc})")
+        return
+    print("Group effect vs Box-B2 (IRR; each subject given its own baseline):")
+    h = f"{'group':20s} {'IRR':>7s} {'~95% CI':>16s} {'mean':>8s} {'SD':>7s} {'~p':>9s}"
+    print(h)
+    print("-" * len(h))
+    for g in OTHERS:
+        c = mixed_contrast(res, gterm(g))
+        ci = f"[{c['lo']:.2f}, {c['hi']:.2f}]"
+        print(f"{g:20s} {c['eff']:7.3f} {ci:>16s} {c['coef']:8.3f} {c['se']:7.3f} {c['p']:9.4f}")
+    print(f"  subject random-intercept SD = {float(np.exp(res.vcp_mean[0])):.3f} "
+          "(log scale); Laplace SEs, treat p-values as approximate.")
+    print("  Compare directions/magnitudes with the GEE table (they should agree).")
+    print("  Single-subject groups (Box-A, Right-Electrode) are partly confounded")
+    print("  with their own random intercept here, so their effects are shrunk.")
+
+
+def plot_curves(df, plot_path):
+    fig, ax = plt.subplots(1, 2, figsize=(12, 4.8))
+    for g in MAIN:
+        s = df[df["group"] == g]
+        d = s.groupby("day").apply(
+            lambda x: pd.Series({"succ": x["success"].mean(),
+                                 "rate": x["success"].sum() / x["total"].sum()}),
+            include_groups=False)
+        ax[0].plot(d.index, d["succ"], marker="o", ms=3, label=g)
+        ax[1].plot(d.index, d["rate"], marker="o", ms=3, label=g)
+    ax[0].set(xlabel="# Days Reach (training day)", ylabel="mean successful reaches",
+              title="(A) Success count learning curves")
+    ax[1].set(xlabel="# Days Reach (training day)", ylabel="success rate (success/attempts)",
+              title="(B) Success-rate learning curves")
+    for a in ax:
+        a.legend(fontsize=8)
+        a.grid(alpha=0.3)
+    fig.tight_layout()
+    fig.savefig(plot_path, dpi=120)
+    print(f"\nSaved learning-curve plot -> {plot_path}")
+
+
+def classify_one(res, unknown):
+    """Compare an unknown group against BOTH anchors; return a label + detail."""
+    vlo = diff_contrast(res, unknown, ANCHOR_LOW)   # unknown vs Box-A2
+    vhi = diff_contrast(res, unknown, ANCHOR_HIGH)  # unknown vs Box-B2
+    like_lo = vlo["p"] >= 0.05   # indistinguishable from Box-A2
+    like_hi = vhi["p"] >= 0.05   # indistinguishable from Box-B2
+    if like_hi and not like_lo:
+        verdict = "B2-like (differs from A2, matches B2)"
+    elif like_lo and not like_hi:
+        verdict = "A2-like (differs from B2, matches A2)"
+    elif like_lo and like_hi:
+        # tie: point toward whichever ratio is closer to 1 on the log scale
+        nearer = ANCHOR_HIGH if abs(vhi["coef"]) < abs(vlo["coef"]) else ANCHOR_LOW
+        verdict = f"AMBIGUOUS (matches both; numerically nearer {nearer})"
+    else:
+        verdict = "UNLIKE BOTH (differs from A2 and B2)"
+    return vlo, vhi, verdict
+
+
+def verdicts(count, rate):
+    print("\n" + "=" * 78)
+    print("VERDICTS")
+    print("=" * 78)
+
+    print("\nAnchor check — H2: Box-B2 BETTER than Box-A2 (the two anchors must differ)")
+    for label, res in [("count/level", count), ("rate/level ", rate)]:
+        c = contrast(res, gterm(ANCHOR_LOW))
+        p1 = one_sided_p_below(c)
+        tag = "SUPPORTED" if (c["coef"] < 0 and p1 < 0.05) else "not supported"
+        print(f"  [{label}] Box-B2 = {1.0 / c['eff']:.2f}x Box-A2  "
+              f"(one-sided p={p1:.4f})  -> {tag}")
+
+    if not UNKNOWNS:
+        print("\n(No unknown groups — they were merged into the anchors; "
+              "only the H2 anchor contrast applies.)")
+        return
+    print("\nClassification — is each UNKNOWN condition A2-like or B2-like?")
+    print("(ratio >1 = above that anchor; p = differs from that anchor)")
+    for u in UNKNOWNS:
+        print(f"\n  {u}:")
+        for label, res in [("count/level", count), ("rate/level ", rate)]:
+            vlo, vhi, verdict = classify_one(res, u)
+            print(f"    [{label}] vs Box-A2: {vlo['eff']:.2f}x "
+                  f"[{vlo['lo']:.2f},{vlo['hi']:.2f}] p={vlo['p']:.3f}   "
+                  f"vs Box-B2: {vhi['eff']:.2f}x [{vhi['lo']:.2f},{vhi['hi']:.2f}] "
+                  f"p={vhi['p']:.3f}")
+            print(f"                 -> {verdict}")
+    print("\n  NOTE: each unknown has ONLY 1 subject. 'Matches' means 'not statistically")
+    print("  distinguishable' — weak evidence at n=1, not proof of equivalence.")
+
+
+def main():
+    ap = argparse.ArgumentParser(description="tDCS GLM hypothesis tests")
+    ap.add_argument("--max-day", type=int, default=None,
+                    help="restrict analysis to training days 0..MAX_DAY (default: all)")
+    ap.add_argument("--merge", action="store_true",
+                    help="preset: Electrode-Box-A->Box-A2 and Right-Electrode->Box-B2")
+    ap.add_argument("--assign", type=str, default=None,
+                    help="custom merges 'UNKNOWN=TARGET,...', e.g. "
+                         "'Electrode-Box-A=Electrode-Box-B2,Right-Electrode=Electrode-Box-B2'")
+    args = ap.parse_args()
+
+    mapping = {}
+    if args.merge:
+        mapping.update(MERGE_MAP)
+    if args.assign:
+        for pair in args.assign.split(","):
+            k, v = pair.split("=")
+            mapping[k.strip()] = v.strip()
+
+    df = load(args.max_day, mapping)
+
+    global MAIN, OTHERS, UNKNOWNS
+    present = set(df["group"])
+    MAIN = [g for g in ALL_GROUPS if g in present]
+    OTHERS = [g for g in MAIN if g != REF]
+    UNKNOWNS = [g for g in CANON_UNKNOWNS if g in present]
+
+    window = f"days 0-{args.max_day}" if args.max_day is not None else "all training days (0+)"
+    tag = f"_d0-{args.max_day}" if args.max_day is not None else ""
+    tag += "_merged" if mapping else ""
+    plot_path = os.path.join(HERE, f"tdcs_learning_curves{tag}.png")
+    mode = ("assign " + ", ".join(f"{k.replace('Electrode-','')}->{v.replace('Electrode-','')}"
+                                  for k, v in mapping.items())) if mapping \
+        else "5 groups (anchors + unknowns)"
+    print(f"ANALYSIS WINDOW: {window}   |   MODE: {mode}   |   observations: {len(df)}")
+    describe(df)
+    count = fit_count_level(df)
+    rate = fit_rate_level(df)
+    learn = fit_learning_rate(df)
+    report_level("(A) LEVEL / COUNT — Poisson GEE on successful reaches",
+                 count, "incidence-rate ratio")
+    report_level("(B) LEVEL / RATE — Binomial GLM on success/attempts (cluster-robust)",
+                 rate, "odds ratio")
+    report_learning(learn)
+    report_anchor_only(df)
+    report_mixed(df)
+    verdicts(count, rate)
+    plot_curves(df, plot_path)
+    print("\nDone. See the module docstring for modeling choices and caveats.")
+
+
+if __name__ == "__main__":
+    main()

analysis/tdcs_reach_data.csv

@@ -0,0 +1,194 @@
+subject,group,day,success,total
+Vu-vuong,Naive,0,18,61
+Vu-vuong,Naive,1,35,91
+Vu-vuong,Naive,2,61,127
+Vu-vuong,Naive,3,34,146
+Vu-vuong,Naive,4,61,141
+Vu-vuong,Naive,5,37,151
+Vu-vuong,Naive,6,49,131
+Vu-vuong,Naive,7,65,149
+Vu-vuong,Naive,8,67,141
+Vu-vuong,Naive,9,73,140
+Vu-vuong,Naive,10,72,131
+Vu-vuong,Naive,11,89,158
+Vu-vuong,Naive,12,97,154
+Vu-vuong,Naive,13,93,145
+Vu-vuong,Naive,14,75,147
+Vu-vuong,Naive,15,68,142
+Vu-vuong,Naive,16,96,162
+Vu-vuong,Naive,17,68,144
+Vu-vuong,Naive,18,81,134
+Vu-vuong,Naive,19,66,144
+Vu-vuong,Naive,20,84,127
+Vu-vuong,Naive,21,66,133
+Khoai-tay-1,Electrode-Box-A,-2,1,18
+Khoai-tay-1,Electrode-Box-A,-1,3,24
+Khoai-tay-1,Electrode-Box-A,0,14,62
+Khoai-tay-1,Electrode-Box-A,1,22,83
+Khoai-tay-1,Electrode-Box-A,2,6,79
+Khoai-tay-1,Electrode-Box-A,3,28,97
+Khoai-tay-1,Electrode-Box-A,4,60,134
+Khoai-tay-1,Electrode-Box-A,5,75,138
+Khoai-tay-1,Electrode-Box-A,6,81,137
+Khoai-tay-1,Electrode-Box-A,7,78,147
+Khoai-tay-1,Electrode-Box-A,8,91,132
+Khoai-tay-1,Electrode-Box-A,9,99,146
+Khoai-tay-1,Electrode-Box-A,10,94,143
+Khoai-tay-1,Electrode-Box-A,11,110,156
+Khoai-tay-1,Electrode-Box-A,12,105,143
+Khoai-tay-1,Electrode-Box-A,13,106,153
+Khoai-tay-1,Electrode-Box-A,14,91,152
+Banh-mi-2,Electrode-Box-A2,-1,2,42
+Banh-mi-2,Electrode-Box-A2,0,25,97
+Banh-mi-2,Electrode-Box-A2,1,22,101
+Banh-mi-2,Electrode-Box-A2,2,22,119
+Banh-mi-2,Electrode-Box-A2,3,29,118
+Banh-mi-2,Electrode-Box-A2,4,27,136
+Banh-mi-2,Electrode-Box-A2,5,43,146
+Banh-mi-2,Electrode-Box-A2,6,74,146
+Banh-mi-2,Electrode-Box-A2,7,70,148
+Banh-mi-2,Electrode-Box-A2,8,65,130
+Banh-mi-2,Electrode-Box-A2,9,79,151
+Banh-mi-2,Electrode-Box-A2,10,93,152
+Egg-tart-2,Electrode-Box-A2,0,9,32
+Egg-tart-2,Electrode-Box-A2,1,2,38
+Egg-tart-2,Electrode-Box-A2,2,31,93
+Egg-tart-2,Electrode-Box-A2,3,44,101
+Egg-tart-2,Electrode-Box-A2,4,54,131
+Egg-tart-2,Electrode-Box-A2,5,84,139
+Egg-tart-2,Electrode-Box-A2,6,85,145
+Egg-tart-2,Electrode-Box-A2,7,79,143
+Egg-tart-2,Electrode-Box-A2,8,76,131
+Egg-tart-2,Electrode-Box-A2,9,88,149
+Egg-tart-2,Electrode-Box-A2,10,81,151
+Egg-tart-2,Electrode-Box-A2,11,78,152
+Egg-tart-2,Electrode-Box-A2,12,96,155
+Egg-tart-2,Electrode-Box-A2,13,84,155
+Root-beer-2,Electrode-Box-A2,0,22,74
+Root-beer-2,Electrode-Box-A2,1,31,87
+Root-beer-2,Electrode-Box-A2,2,49,134
+Root-beer-2,Electrode-Box-A2,3,31,89
+Root-beer-2,Electrode-Box-A2,4,60,140
+Root-beer-2,Electrode-Box-A2,5,84,147
+Banh-mi-1,Electrode-Box-B2,0,13,84
+Banh-mi-1,Electrode-Box-B2,1,35,86
+Banh-mi-1,Electrode-Box-B2,2,47,110
+Banh-mi-1,Electrode-Box-B2,3,65,140
+Banh-mi-1,Electrode-Box-B2,4,70,127
+Banh-mi-1,Electrode-Box-B2,5,102,142
+Banh-mi-1,Electrode-Box-B2,6,90,131
+Banh-mi-1,Electrode-Box-B2,7,109,148
+Banh-mi-1,Electrode-Box-B2,8,104,137
+Banh-mi-1,Electrode-Box-B2,9,119,150
+Banh-mi-1,Electrode-Box-B2,10,121,158
+Banh-mi-1,Electrode-Box-B2,11,121,148
+Banh-mi-1,Electrode-Box-B2,12,120,149
+Banh-mi-1,Electrode-Box-B2,13,135,154
+Egg-tart-1,Electrode-Box-B2,0,7,56
+Egg-tart-1,Electrode-Box-B2,1,16,78
+Egg-tart-1,Electrode-Box-B2,2,23,103
+Egg-tart-1,Electrode-Box-B2,3,63,120
+Egg-tart-1,Electrode-Box-B2,4,69,132
+Egg-tart-1,Electrode-Box-B2,5,83,136
+Egg-tart-1,Electrode-Box-B2,6,71,142
+Egg-tart-1,Electrode-Box-B2,7,79,138
+Egg-tart-1,Electrode-Box-B2,8,98,142
+Egg-tart-1,Electrode-Box-B2,9,89,139
+Egg-tart-1,Electrode-Box-B2,10,96,143
+Egg-tart-1,Electrode-Box-B2,11,96,148
+Egg-tart-1,Electrode-Box-B2,12,101,156
+Egg-tart-1,Electrode-Box-B2,13,103,152
+Egg-tart-1,Electrode-Box-B2,14,97,152
+Root-beer-1,Electrode-Box-B2,0,11,85
+Root-beer-1,Electrode-Box-B2,1,18,76
+Root-beer-1,Electrode-Box-B2,2,40,105
+Root-beer-1,Electrode-Box-B2,3,55,134
+Root-beer-1,Electrode-Box-B2,4,75,136
+Root-beer-1,Electrode-Box-B2,5,64,133
+Root-beer-1,Electrode-Box-B2,6,104,139
+Root-beer-1,Electrode-Box-B2,7,98,148
+Root-beer-1,Electrode-Box-B2,8,81,145
+Root-beer-1,Electrode-Box-B2,9,89,156
+Root-beer-1,Electrode-Box-B2,10,105,158
+Khoai-lang-2,Naive,0,0,0
+Khoai-lang-2,Naive,1,0,0
+Khoai-lang-2,Naive,2,10,47
+Khoai-lang-2,Naive,3,11,52
+Khoai-lang-2,Naive,4,9,56
+Khoai-lang-2,Naive,5,34,95
+Khoai-lang-2,Naive,6,21,72
+Khoai-lang-2,Naive,7,23,99
+Khoai-lang-2,Naive,8,64,136
+Khoai-lang-2,Naive,9,75,131
+Khoai-lang-2,Naive,10,63,134
+Khoai-lang-2,Naive,11,59,139
+Khoai-lang-2,Naive,12,51,129
+Khoai-lang-2,Naive,13,73,143
+Khoai-lang-2,Naive,14,82,136
+Khoai-lang-2,Naive,15,70,145
+Khoai-lang-2,Naive,16,76,135
+Khoai-lang-2,Naive,17,76,150
+Khoai-lang-2,Naive,18,63,122
+Khoai-lang-2,Naive,19,48,116
+Khoai-lang-2,Naive,20,65,134
+Khoai-lang-2,Naive,21,75,131
+Khoai-lang-2,Naive,22,98,146
+Khoai-lang-2,Naive,23,88,139
+Khoai-lang-2,Naive,24,94,148
+Khoai-lang-2,Naive,25,56,102
+Khoai-lang-2,Naive,26,75,143
+Khoai-tay-2,Naive,1,21,68
+Khoai-tay-2,Naive,2,6,79
+Khoai-tay-2,Naive,3,0,83
+Khoai-tay-2,Naive,4,28,87
+Khoai-tay-2,Naive,5,31,125
+Khoai-tay-2,Naive,6,62,144
+Khoai-tay-2,Naive,7,87,141
+Khoai-tay-2,Naive,8,105,152
+Khoai-tay-2,Naive,9,75,148
+Khoai-tay-2,Naive,10,101,149
+Khoai-tay-2,Naive,11,102,154
+Khoai-tay-2,Naive,12,95,144
+Khoai-tay-2,Naive,13,77,149
+Khoai-tay-2,Naive,14,91,153
+OM-2,Naive,-2,0,0
+OM-2,Naive,-1,5,18
+OM-2,Naive,0,1,7
+OM-2,Naive,1,11,33
+OM-2,Naive,2,19,62
+OM-2,Naive,3,29,117
+OM-2,Naive,4,73,126
+OM-2,Naive,5,53,135
+OM-2,Naive,6,73,138
+OM-2,Naive,7,80,131
+OM-2,Naive,8,91,141
+OM-2,Naive,9,90,135
+OM-2,Naive,10,95,142
+OM-2,Naive,11,60,133
+OM-2,Naive,12,58,142
+Khoai-lang-1,Right-Electrode,-4,0,0
+Khoai-lang-1,Right-Electrode,-3,2,7
+Khoai-lang-1,Right-Electrode,-1,0,0
+Khoai-lang-1,Right-Electrode,0,3,69
+Khoai-lang-1,Right-Electrode,1,18,97
+Khoai-lang-1,Right-Electrode,2,27,84
+Khoai-lang-1,Right-Electrode,3,46,121
+Khoai-lang-1,Right-Electrode,4,65,143
+Khoai-lang-1,Right-Electrode,5,76,141
+Khoai-lang-1,Right-Electrode,6,79,146
+Khoai-lang-1,Right-Electrode,7,81,153
+Khoai-lang-1,Right-Electrode,8,98,144
+Khoai-lang-1,Right-Electrode,9,101,149
+Khoai-lang-1,Right-Electrode,10,103,150
+Khoai-lang-1,Right-Electrode,11,113,154
+Khoai-lang-1,Right-Electrode,12,117,148
+Khoai-lang-1,Right-Electrode,13,114,146
+Khoai-lang-1,Right-Electrode,14,119,132
+Khoai-lang-1,Right-Electrode,15,111,156
+Khoai-lang-1,Right-Electrode,16,74,156
+Khoai-lang-1,Right-Electrode,17,109,139
+Khoai-lang-1,Right-Electrode,18,88,138
+Khoai-lang-1,Right-Electrode,19,98,140
+Khoai-lang-1,Right-Electrode,20,118,153
+Khoai-lang-1,Right-Electrode,21,109,146
+Khoai-lang-1,Right-Electrode,22,106,141