docs(spec): MATLAB port of the tDCS GLM analysis

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
2026-07-20 07:14:59 -04:00
parent 5314071c6b
commit 2e6464cd14
@@ -0,0 +1,143 @@
# MATLAB port of the tDCS GLM analysis — Design
**Date:** 2026-07-20
**Status:** Approved (design)
**Location:** `analysis/matlab/`
**Ports:** the Python analysis in `analysis/` (`tdcs_glm.py`, `METHODS.md`) to MATLAB.
## Goal
Reproduce the tDCS reaching analysis in MATLAB: consume raw exports from the app,
derive the day-window × merge-scenario variations, run each scenario
independently via `switch/case`, and print the raw GLM output plus a
plain-language interpretation. Include a `matlab.unittest` test suite.
## Environment
- MATLAB **R2025b** + **Statistics and Machine Learning Toolbox** installed at
`~/matlab` (via `mpm`). Requires the user's license to be **activated** before
anything runs (`~/matlab/bin/activate_matlab.sh`, or a license file in
`~/matlab/licenses/`).
- Headless execution: `~/matlab/bin/matlab -batch "<cmd>"`. Tests run via
`matlab -batch "run_tests"`. **The MATLAB code and tests cannot be executed in
the dev environment until the license is active** — they are written first and
run once licensing is in place.
## Modeling engine
MATLAB has no GEE; the port uses **`fitglme`** (generalized linear mixed model)
with a **per-subject random intercept** as the single primary estimator — this
folds the Python GEE + random-intercept models into one native tool. The Python
mixed model already agreed with its GEE, so results should track. Reference group
is **Electrode-Box-B2**; `day_c` is the mean-centered training day and
`day_c2 = day_c^2` captures the plateau.
Models per scenario (on the scenario's data table):
- **Count:** `success ~ group + day_c + day_c2 + (1|subject)`, `Distribution='Poisson'` → IRRs vs B2.
- **Rate:** `success ~ group + day_c + day_c2 + (1|subject)`, `Distribution='Binomial'`, `BinomialSize=total` → odds ratios. Zero-attempt sessions (`total==0`, undefined rate) are dropped for this model only; the count model keeps them (they are genuine 0-success observations).
- **Learning rate:** `success ~ group*day_c + day_c2 + (1|subject)`, Poisson → group×day slope differences + a joint `coefTest`.
- **Per-animal:** collapse each animal to one number (mean count; pooled rate = Σsuccess/Σtotal), then `ttest2` (Welch) and `ranksum` (Mann-Whitney, one-sided) — the honest small-n check.
- **`unmerged_*` scenarios only:** classify Box-A and Right-Electrode against **both** anchors via `coefTest` contrasts → A2-like / B2-like / ambiguous.
## Data flow
```
raw/raw_success.csv ┐
raw/raw_total.csv ┴─ tdcs_load_data → long table {subject,group,day,success,total}
→ tdcs_scenario_data(scenario)
→ filtered+merged table (saved: derived/<scenario>.csv)
→ fitglme models → console + results/<scenario>.txt
```
### Raw files (produced now, in the website's export format)
Two app-style matrices, identical to the app's **Export Data** output with
X = "# Days Reach", Group by = "Group":
- `raw/raw_success.csv` — Data = Success
- `raw/raw_total.csv` — Data = Total attempts
Format (both):
```
Data,<metric>
<blank>
Group,<group per subject column>
# Days Reach,<subject headers>
<day>,<value per subject>
...
```
These are generated from the running app (via the same export logic) and are
byte-for-byte what the UI would export; the user can regenerate them from the UI.
### `tdcs_load_data.m`
Parses the matrix format (skips the `Data,` line and blank line, reads the
`Group` row, the `# Days Reach` + subject-header row, then the day rows),
reshapes both matrices to long, and joins on (subject, day) to produce one table
with `subject, group, day, success, total`. Drops blank cells; applies day ≥ 0.
Returns the full long table (all 5 groups, Box-A included).
## Scenarios (`switch/case`, 6 self-contained cases)
`tdcs_glm(scenario)` dispatches on the scenario name; each case runs end-to-end
and can be run alone. Cases:
| case | merge | window |
|---|---|---|
| `unmerged_full` | 5 groups (classify Box-A & Right) | days 026 |
| `unmerged_d0_10` | 5 groups | days 010 |
| `mergeA2_full` | Box-A→A2, Right→B2 | 026 |
| `mergeA2_d0_10` | Box-A→A2, Right→B2 | 010 |
| `mergeB2_full` | Box-A→B2, Right→B2 | 026 |
| `mergeB2_d0_10` | Box-A→B2, Right→B2 | 010 |
Merge maps: `mergeA2 = {Electrode-Box-A→Electrode-Box-A2, Right-Electrode→Electrode-Box-B2}`;
`mergeB2 = {Electrode-Box-A→Electrode-Box-B2, Right-Electrode→Electrode-Box-B2}`.
Each case saves its derived scenario table to `derived/<scenario>.csv`.
`run_all.m` loops all six; `run_tests.m` runs the test suite.
## Printed output per scenario
To console and `results/<scenario>.txt`:
1. **Descriptives** — per group: n_subjects, n_sessions, mean success count, pooled rate.
2. **Raw GLM** — the literal `fitglme` fixed-effects table (Estimate, SE, tStat, pValue, CI) for count and rate models, plus the learning-rate joint test statistic/p.
3. **Interpretation** — plain language mirroring the Python verdicts: H2 (B2 vs A2) as IRR/OR + one-sided p + SUPPORTED/not; per-animal p-values; classification (unmerged only); and the standing caveats (n=34/group, single-subject fragility, count vs rate).
## Files
```
analysis/matlab/
raw/raw_success.csv, raw/raw_total.csv # generated inputs (website format)
tdcs_load_data.m # parse + reshape + join → long table
tdcs_scenario_data.m # apply merge map + day window; save derived
tdcs_models.m # fit GLMMs, per-animal tests, classification, contrasts
tdcs_report.m # format descriptives + raw GLM + interpretation
tdcs_glm.m # switch/case entry over the 6 scenarios
run_all.m, run_tests.m
derived/ results/ # outputs (created at runtime)
tests/
tLoadData.m # join correctness (verified cells), day filter, zero-attempt drop
tScenarioData.m # per-scenario groups/subject counts, window caps
tModels.m # golden-value/direction checks vs the Python numbers (loose tolerance)
```
## Testing
`matlab.unittest`, run headless via `matlab -batch "run_tests"`:
- **Data:** `tdcs_load_data` join correct — verified cells (e.g. Vu-vuong day0 success=18, total=61; Banh-mi-1 day0 success=13) and per-group subject counts; day≥0 and zero-attempt handling.
- **Scenario derivation:** each scenario yields the right groups/subjects (`mergeB2` → B2=5 subjects, A2=3; `mergeA2` → 4/4/4; `d0_10` caps at day 10) and saves `derived/<scenario>.csv`.
- **Model sanity / golden values:** GLM results land near the Python numbers within a loose tolerance and match direction — e.g. `mergeB2_full` rate H2 odds ≈ 1.41.6 with B2>A2; anchor H2 supported on rate; `unmerged` classification labels. GLMM ≠ GEE exactly, so magnitude tolerances are wide but sign/ordering checks are strict.
Executed once the license is active; written but unrun until then (flagged in the report).
## Out of scope (YAGNI)
- Re-implementing GEE / cluster-robust sandwich SEs (using native `fitglme` instead).
- The learning-curve plots (the Python `.png`s remain; not re-created in MATLAB unless asked).
- Changing the Python analysis (it stays as the reference).
- Bayesian mixed model (the Python MAP sensitivity has no direct `fitglme` analog and is dropped; the GLMM *is* the random-intercept model).