Experiments DB Dev
152 lines
13 KiB
Markdown
152 lines
13 KiB
Markdown
# MATLAB tDCS Analysis Port — Implementation Plan
|
||
|
||
> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax.
|
||
|
||
**Goal:** Port the Python tDCS GLM analysis (`analysis/tdcs_glm.py`) to MATLAB under `analysis/matlab/`, driven by raw app exports, with 6 self-contained `switch/case` scenarios that print raw GLM output + interpretation, plus a `matlab.unittest` suite.
|
||
|
||
**Architecture:** `fitglme` (GLMM, per-subject random intercept) replaces the Python GEE. Two raw matrices (Success, Total attempts) are parsed and joined to a long table; a curated 12-animal map defines groups; each scenario applies a merge map + day window, fits the models, and reports.
|
||
|
||
**Tech Stack:** MATLAB R2025b + Statistics and Machine Learning Toolbox (installed at `~/matlab`).
|
||
|
||
## Global Constraints
|
||
|
||
- **MATLAB is NOT runnable in this environment until the user's license is activated** (`~/matlab/bin/activate_matlab.sh`). Write all code and tests; they are executed later via `~/matlab/bin/matlab -batch "run_tests"`. In each task, "run the test" means: **the test is authored and self-consistent; note it is unrun pending license.** Do not claim tests pass.
|
||
- Reference implementation: `analysis/tdcs_glm.py` and `analysis/METHODS.md`. Port equivalent behavior; do not change the Python.
|
||
- Reference group = **`Electrode-Box-B2`**; `day_c` = mean-centered day; `day_c2 = day_c.^2`.
|
||
- Analysis restricted to **day ≥ 0**. The **rate model only** drops `total==0` sessions; the count model keeps them.
|
||
- **Curated 12-animal group map** (by animal name; all other animals in the raw files are dropped):
|
||
- Naive: `Vu-vuong, Khoai-lang-2, Khoai-tay-2, OM-2`
|
||
- Electrode-Box-A: `Khoai-tay-1`
|
||
- Electrode-Box-A2: `Banh-mi-2, Egg-tart-2, Root-beer-2`
|
||
- Electrode-Box-B2: `Banh-mi-1, Egg-tart-1, Root-beer-1`
|
||
- Right-Electrode: `Khoai-lang-1`
|
||
- 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}`; `unmerged = {}` (5 groups, classify unknowns).
|
||
- Raw files already generated at `analysis/matlab/raw/raw_success.csv` and `raw_total.csv` (app export format: `Data,<metric>` / blank / `Group,…` / `# Days Reach,<subjects>` / `<day>,<values>`). Contain all 15 animals incl. negative days.
|
||
- Load-check cells (post-join, day≥0): `Vu-vuong` day 0 → success=18, total=61; `Banh-mi-1` day 0 → success=13, total=84.
|
||
|
||
---
|
||
|
||
### Task 1: Data loading + curated map
|
||
|
||
**Files:**
|
||
- Create: `analysis/matlab/tdcs_config.m` (constants: analysis group map, merge maps, ref group)
|
||
- Create: `analysis/matlab/tdcs_load_data.m`
|
||
- Create: `analysis/matlab/tests/tLoadData.m`
|
||
|
||
**Interfaces:**
|
||
- Produces `cfg = tdcs_config()` returning a struct: `cfg.groups` (containers.Map name→group for the 12), `cfg.ref='Electrode-Box-B2'`, `cfg.mergeA2`, `cfg.mergeB2` (containers.Map), `cfg.anchorLow='Electrode-Box-A2'`, `cfg.anchorHigh='Electrode-Box-B2'`, `cfg.unknowns={'Electrode-Box-A','Right-Electrode'}`.
|
||
- Produces `T = tdcs_load_data()` returning a table with variables `subject` (string), `group` (categorical), `day` (double), `success` (double), `total` (double); only the 12 curated animals; only `day >= 0`.
|
||
|
||
**Steps:**
|
||
- [ ] **Step 1: Write `tdcs_config.m`** — return the struct above using the Global Constraints values. `containers.Map` for the name→group map (all 12 entries) and the two merge maps.
|
||
- [ ] **Step 2: Write the failing test `tests/tLoadData.m`** (matlab.unittest). Assertions:
|
||
- `T` has 12 unique subjects; `height(T)` matches the non-blank day≥0 rows.
|
||
- Group subject counts: Naive=4, Electrode-Box-A2=3, Electrode-Box-B2=3, Electrode-Box-A=1, Right-Electrode=1.
|
||
- `T(T.subject=="Vu-vuong" & T.day==0, {'success','total'})` equals `[18 61]`.
|
||
- `T(T.subject=="Banh-mi-1" & T.day==0, {'success','total'})` equals `[13 84]`.
|
||
- No rows with `day < 0`.
|
||
- [ ] **Step 3: Write `tdcs_load_data.m`.** Port the CSV-matrix parsing:
|
||
- Read `raw/raw_success.csv` and `raw/raw_total.csv` as raw lines (`readlines` or `fileread`+`splitlines`).
|
||
- Line 1 = `Data,<metric>` (skip); line 2 blank (skip); line 3 = `Group,<…>` (skip — grouping comes from the curated map, not this row); line 4 = `# Days Reach,<subject headers>` → subject names (columns 2:end); lines 5+ = `<day>,<value…>`.
|
||
- Build long form: for each (day, subject) cell that is non-empty, emit a row. Parse success matrix and total matrix separately, then `innerjoin`/`outerjoin` on (subject, day). A cell present in one metric but blank in the other → treat blank as missing; keep the row only if success present (matches Python: rows come from Success; total is fully populated where success is).
|
||
- Apply the curated map: keep only subjects in `cfg.groups`; assign `group` from the map (overriding the file's Group row). Drop `day < 0`. Return the table.
|
||
- [ ] **Step 4: (Deferred) run** `~/matlab/bin/matlab -batch "run_tests"` once licensed; expect `tLoadData` green. For now, self-review the parsing against the two raw files and the load-check cells by reading.
|
||
- [ ] **Step 5: Commit** `git add analysis/matlab && git commit -m "matlab(tdcs): data loading + curated group map"`
|
||
|
||
---
|
||
|
||
### Task 2: Scenario derivation
|
||
|
||
**Files:**
|
||
- Create: `analysis/matlab/tdcs_scenario_data.m`
|
||
- Create: `analysis/matlab/tests/tScenarioData.m`
|
||
|
||
**Interfaces:**
|
||
- Consumes `tdcs_load_data`, `tdcs_config`.
|
||
- Produces `[S, meta] = tdcs_scenario_data(scenario)` where `scenario` ∈ {`unmerged_full`,`unmerged_d0_10`,`mergeA2_full`,`mergeA2_d0_10`,`mergeB2_full`,`mergeB2_d0_10`}. `S` = scenario table (group remapped, day window applied, `day_c`/`day_c2` added, `group` categorical with `cfg.ref` as the FIRST/reference category). `meta` = struct with `.mergeKey`, `.maxDay`, `.isUnmerged`. Also **saves** `S` to `analysis/matlab/derived/<scenario>.csv` (create the folder if missing).
|
||
|
||
**Steps:**
|
||
- [ ] **Step 1: Write the failing test `tests/tScenarioData.m`.** For each scenario, assert the per-group subject counts match this table, and that `max(S.day) <= maxDay`:
|
||
|
||
| scenario | subj/group | maxDay |
|
||
|---|---|---|
|
||
| unmerged_full | A=1,A2=3,B2=3,Naive=4,Right=1 | 26 |
|
||
| unmerged_d0_10 | A=1,A2=3,B2=3,Naive=4,Right=1 | 10 |
|
||
| mergeA2_full | A2=4,B2=4,Naive=4 | 26 |
|
||
| mergeA2_d0_10 | A2=4,B2=4,Naive=4 | 10 |
|
||
| mergeB2_full | A2=3,B2=5,Naive=4 | 26 |
|
||
| mergeB2_d0_10 | A2=3,B2=5,Naive=4 | 10 |
|
||
|
||
Also assert `derived/mergeB2_full.csv` is written and reloads to the same subject/group counts, and that the reference level of `S.group` is `Electrode-Box-B2`.
|
||
- [ ] **Step 2: Write `tdcs_scenario_data.m`.** Parse `scenario` into `(mergeKey, maxDay)` by splitting on `_full`/`_d0_10`. Load `T`, apply the merge `containers.Map` (relabel groups; `unmerged` = identity), filter `day <= maxDay`, center day, add `day_c2`, set `group` categorical with `reordercats` so `cfg.ref` is first. Write `derived/<scenario>.csv` via `writetable`.
|
||
- [ ] **Step 3: (Deferred) run tests once licensed.** Self-review the counts against the table by reading.
|
||
- [ ] **Step 4: Commit** `git commit -am "matlab(tdcs): scenario derivation + derived data export"`
|
||
|
||
---
|
||
|
||
### Task 3: Models (GLMM + per-animal + classification)
|
||
|
||
**Files:**
|
||
- Create: `analysis/matlab/tdcs_models.m`
|
||
- Create: `analysis/matlab/tests/tModels.m`
|
||
|
||
**Interfaces:**
|
||
- Consumes `tdcs_scenario_data`, `tdcs_config`.
|
||
- Produces `R = tdcs_models(S, meta, cfg)` returning a struct of results:
|
||
- `R.countGLME`, `R.rateGLME`, `R.learnGLME` — the `GeneralizedLinearMixedModel` objects.
|
||
- `R.perAnimal` — struct with `.countWelchP, .countMWUp, .rateWelchP, .rateMWUp, .nB, .nA` (B2 vs A2).
|
||
- `R.h2` — struct with `.countIRR, .countOneSidedP, .rateOR, .rateOneSidedP` (Box-A2 vs B2 from the GLMEs).
|
||
- `R.classify` — (unmerged only) struct per unknown with vs-A2 and vs-B2 ratios/p and a label; empty otherwise.
|
||
|
||
**Steps:**
|
||
- [ ] **Step 1: Write the failing test `tests/tModels.m`** using these Python golden values (per-animal tests are pure stats and must match to 3 decimals; GLME checks are loose direction/magnitude):
|
||
|
||
Per-animal (B2 vs A2), `assertEqual(..., 'AbsTol', 0.005)`:
|
||
| scenario | nB,nA | rateWelchP | rateMWUp | countWelchP | countMWUp |
|
||
|---|---|---|---|---|---|
|
||
| unmerged_full | 3,3 | 0.072 | 0.050 | 0.055 | 0.050 |
|
||
| mergeA2_full | 4,4 | 0.056 | 0.029 | 0.035 | 0.029 |
|
||
| mergeB2_full | 5,3 | 0.044 | 0.018 | 0.020 | 0.018 |
|
||
| unmerged_d0_10 | 3,3 | 0.070 | 0.050 | 0.041 | 0.050 |
|
||
| mergeA2_d0_10 | 4,4 | 0.069 | 0.057 | 0.023 | 0.014 |
|
||
| mergeB2_d0_10 | 5,3 | 0.094 | 0.071 | 0.023 | 0.018 |
|
||
|
||
GLME sanity (loose): for `mergeB2_full`, `R.h2.rateOR` in `[1.3, 1.7]` and `R.h2.rateOneSidedP < 0.05` and `R.h2.countIRR > 1`. For `unmerged_full`, `R.classify` has entries for both unknowns.
|
||
- [ ] **Step 2: Write `tdcs_models.m`.** Port the models:
|
||
- Count: `fitglme(S, 'success ~ group + day_c + day_c2 + (1|subject)', 'Distribution','Poisson','Link','log')`.
|
||
- Rate: drop `total==0` rows, then `fitglme(Sr, 'success ~ group + day_c + day_c2 + (1|subject)', 'Distribution','Binomial','BinomialSize',Sr.total,'Link','logit')`.
|
||
- Learning: `fitglme(S, 'success ~ group*day_c + day_c2 + (1|subject)', 'Distribution','Poisson')`; joint test on the interaction terms via `coefTest`.
|
||
- H2: from each GLME read the `[T.Electrode-Box-A2]` coefficient (vs B2 ref); IRR/OR = `exp(coef)`; one-sided p = `normcdf(coef/se)` (tests coef<0). Report `1/exp(coef)` as the B2 advantage.
|
||
- Per-animal: build per-subject summaries (mean success; pooled Σsuccess/Σtotal), then Welch via `ttest2(b,a,'Vartype','unequal')` and one-sided MWU via `ranksum(b,a,'tail','right')`. (Direction: B2 vs A2.)
|
||
- Classification (unmerged): for each unknown, contrast vs A2 and vs B2 (use `coefTest` with a contrast vector, or read coefficients + covariance from `R.<model>.CoefficientCovariance`), label A2-like/B2-like/ambiguous per METHODS.md.
|
||
- [ ] **Step 3: (Deferred) run tests once licensed.** Self-review the estimator calls and the per-animal math against the Python.
|
||
- [ ] **Step 4: Commit** `git commit -am "matlab(tdcs): GLMM + per-animal + classification models"`
|
||
|
||
---
|
||
|
||
### Task 4: Report, entry point, runners
|
||
|
||
**Files:**
|
||
- Create: `analysis/matlab/tdcs_report.m`, `analysis/matlab/tdcs_glm.m`, `analysis/matlab/run_all.m`, `analysis/matlab/run_tests.m`
|
||
- Create: `analysis/matlab/README.md`
|
||
|
||
**Interfaces:**
|
||
- Consumes all prior. Produces `tdcs_glm(scenario)` — the `switch/case` entry that validates the scenario name, calls `tdcs_scenario_data` → `tdcs_models` → `tdcs_report`, prints to console and writes `results/<scenario>.txt`.
|
||
|
||
**Steps:**
|
||
- [ ] **Step 1: Write `tdcs_report.m`** — given `S`, `R`, `meta`, `cfg`, print: (1) descriptives (per group: n_subj, n_sessions, mean success, pooled rate); (2) raw GLM — `disp(R.countGLME)` and `disp(R.rateGLME)` fixed-effects tables + the learning joint test; (3) interpretation mirroring `verdicts()`/`classify()` in the Python (H2 verdict with IRR/OR + one-sided p; per-animal p's; classification for unmerged; the standing caveats). Use `fprintf`; also capture to `results/<scenario>.txt` via `diary` or an output string.
|
||
- [ ] **Step 2: Write `tdcs_glm.m`** with an explicit `switch scenario` over the 6 names (each case is a self-contained call; `otherwise` errors with the valid list). Even though cases share the pipeline, enumerate all 6 so each is independently invocable and greppable, per the spec.
|
||
- [ ] **Step 3: Write `run_all.m`** (loop the 6 scenarios calling `tdcs_glm`) and `run_tests.m` (`results = runtests('analysis/matlab/tests'); disp(results);` with a nonzero exit on failure for `-batch`).
|
||
- [ ] **Step 4: Write `README.md`** — how to activate the license, run one scenario (`matlab -batch "tdcs_glm('mergeB2_full')"`), run all, run tests; note results were validated against the Python golden values and that GLMM ≠ GEE (interpretation matches, exact numbers differ).
|
||
- [ ] **Step 5: (Deferred) once licensed**, run `matlab -batch "run_tests"` and `matlab -batch "run_all"`, capture output, fix failures. Until then, self-review by reading.
|
||
- [ ] **Step 6: Commit** `git commit -am "matlab(tdcs): report, switch/case entry, runners, README"`
|
||
|
||
---
|
||
|
||
## Self-Review Notes
|
||
|
||
- **Spec coverage:** raw files in website format (Task 1 input, pre-generated) · parse+join+curated map (Task 1) · day windows + merge scenarios saved (Task 2) · `fitglme` count/rate/learning + per-animal + classification (Task 3) · switch/case per scenario + raw GLM + interpretation (Task 4) · `matlab.unittest` suite with golden values (Tasks 1–3). All covered.
|
||
- **Runtime caveat is pervasive:** every task defers actual execution to post-licensing; reviews are by reading + golden-value cross-checks, and no task may claim tests pass.
|
||
- **Golden values** are authoritative (computed from the Python this session). Per-animal stats must match to 3 decimals; GLME magnitudes are checked loosely because GLMM ≠ GEE.
|
||
- **Out of scope (per spec):** GEE/sandwich SEs, MATLAB plots, Bayesian MAP model, editing the Python or the DB.
|