experiments-database/analysis/METHODS.md

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.