Compare commits
8 Commits
5314071c6b
...
13414d31fc
| Author | SHA1 | Date | |
|---|---|---|---|
 Experiments DB Dev
|
13414d31fc | ||
|
Experiments DB Dev
|
01f59a5265 | ||
|
Experiments DB Dev
|
27d388530a | ||
|
Experiments DB Dev
|
75ea8a031e | ||
|
Experiments DB Dev
|
e398237c76 | ||
|
Experiments DB Dev
|
33121e3a8f | ||
|
Experiments DB Dev
|
808040a522 | ||
|
Experiments DB Dev
|
2e6464cd14 |
@@ -0,0 +1,109 @@
|
|||||||
|
# tDCS reaching study — MATLAB GLM port
|
||||||
|
|
||||||
|
This directory is a MATLAB port of the Python tDCS GLM analysis
|
||||||
|
(`analysis/tdcs_glm.py`, documented in `analysis/METHODS.md`), driven by raw
|
||||||
|
app-export matrices rather than the pre-joined `tdcs_reach_data.csv`. It
|
||||||
|
reproduces the same three models (count level, rate level, learning rate) and
|
||||||
|
the same H2 anchor check / unknown-condition classification, but replaces the
|
||||||
|
Python's population-average GEE with a `fitglme` GLMM (Poisson/Binomial,
|
||||||
|
subject random intercept) — see "GLMM vs GEE" below.
|
||||||
|
|
||||||
|
## Requirements
|
||||||
|
|
||||||
|
- MATLAB R2025b + Statistics and Machine Learning Toolbox.
|
||||||
|
- A license activated at `~/matlab/licenses/license.dat` (already in place in
|
||||||
|
this environment). If MATLAB ever reports it is unlicensed, activate with:
|
||||||
|
|
||||||
|
```
|
||||||
|
~/matlab/bin/activate_matlab.sh
|
||||||
|
```
|
||||||
|
|
||||||
|
## Layout
|
||||||
|
|
||||||
|
- `tdcs_config.m` — curated 12-animal group map, reference group, merge maps.
|
||||||
|
- `tdcs_load_data.m` — parses `raw/raw_success.csv` and `raw/raw_total.csv`
|
||||||
|
(app-export matrix format), joins them, and applies the curated group map.
|
||||||
|
- `tdcs_scenario_data.m` — builds one of the 6 scenarios (merge map + day
|
||||||
|
window), and saves the scenario table to `derived/<scenario>.csv`.
|
||||||
|
- `tdcs_models.m` — fits the count/rate/learning GLMEs, the per-animal
|
||||||
|
Welch/Mann-Whitney tests, the H2 anchor check, and (unmerged scenarios
|
||||||
|
only) the unknown-condition classification.
|
||||||
|
- `tdcs_report.m` — prints descriptives, the raw GLME fixed-effects tables,
|
||||||
|
and an interpretation section (H2 verdict, per-animal p-values,
|
||||||
|
classification, caveats); also writes `results/<scenario>.txt`.
|
||||||
|
- `tdcs_glm.m` — the entry point: `tdcs_glm(scenario)` runs
|
||||||
|
`tdcs_scenario_data` → `tdcs_models` → `tdcs_report` for one named
|
||||||
|
scenario, via an explicit `switch/case` over all 6 scenario names.
|
||||||
|
- `run_all.m` — runs all 6 scenarios in one script.
|
||||||
|
- `run_tests.m` — runs the `matlab.unittest` suite in `tests/` and exits
|
||||||
|
nonzero on any failure.
|
||||||
|
- `tests/` — `matlab.unittest` test classes (`tLoadData`, `tScenarioData`,
|
||||||
|
`tModels`) with golden values cross-checked against the Python.
|
||||||
|
- `raw/` — pre-generated app-export CSVs (input; not modified by this code).
|
||||||
|
- `derived/` — scenario tables written by `tdcs_scenario_data` (generated).
|
||||||
|
- `results/` — report text files written by `tdcs_report`/`tdcs_glm`
|
||||||
|
(generated).
|
||||||
|
|
||||||
|
## Running
|
||||||
|
|
||||||
|
Run one scenario (valid names: `unmerged_full`, `unmerged_d0_10`,
|
||||||
|
`mergeA2_full`, `mergeA2_d0_10`, `mergeB2_full`, `mergeB2_d0_10`):
|
||||||
|
|
||||||
|
```
|
||||||
|
cd analysis/matlab
|
||||||
|
~/matlab/bin/matlab -batch "addpath(pwd); tdcs_glm('mergeB2_full')"
|
||||||
|
```
|
||||||
|
|
||||||
|
This prints the full report to the console and writes it to
|
||||||
|
`results/mergeB2_full.txt`.
|
||||||
|
|
||||||
|
Run all 6 scenarios:
|
||||||
|
|
||||||
|
```
|
||||||
|
~/matlab/bin/matlab -batch "addpath(pwd); run_all"
|
||||||
|
```
|
||||||
|
|
||||||
|
Run the test suite:
|
||||||
|
|
||||||
|
```
|
||||||
|
~/matlab/bin/matlab -batch "addpath(pwd); run_tests"
|
||||||
|
```
|
||||||
|
|
||||||
|
(`run_tests` calls `error(...)` on any failure, so the `matlab -batch`
|
||||||
|
process exits with a nonzero status — suitable for CI.)
|
||||||
|
|
||||||
|
## Validation against the Python
|
||||||
|
|
||||||
|
Per-animal statistics (per-subject Welch t-test / one-sided Mann-Whitney U,
|
||||||
|
Box-B2 vs Box-A2) are pure stats independent of the GLM choice, and were
|
||||||
|
matched to the Python's golden values to 3 decimals across all 6 scenarios
|
||||||
|
(see `docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md`, Task 3).
|
||||||
|
|
||||||
|
## GLMM vs GEE (read this before comparing numbers to the Python)
|
||||||
|
|
||||||
|
The Python reference uses **GEE** (population-average effects, exchangeable
|
||||||
|
working correlation, robust/sandwich SEs) for the count and learning-rate
|
||||||
|
models, and a cluster-robust Binomial GLM for the rate model. This MATLAB
|
||||||
|
port instead uses **GLMM** (`fitglme`, subject-specific effects via a
|
||||||
|
`(1|subject)` random intercept, Laplace/PL approximation) for all three
|
||||||
|
models, because MATLAB's Statistics and Machine Learning Toolbox has no GEE
|
||||||
|
implementation.
|
||||||
|
|
||||||
|
GLMM and GEE estimate related but distinct quantities (subject-specific vs.
|
||||||
|
population-average effects), so **exact coefficients, ratios, and p-values
|
||||||
|
will differ** between the two implementations — sometimes by a fair amount
|
||||||
|
for the small, unbalanced groups in this study. What should agree is the
|
||||||
|
**interpretation**: the direction of the H2 anchor effect, which unknown
|
||||||
|
condition classifies as A2-like/B2-like/ambiguous, and whether the learning
|
||||||
|
rates differ across groups. Treat the MATLAB numbers as a sensitivity check
|
||||||
|
on the Python's conclusions, not as a numerically-identical reproduction.
|
||||||
|
|
||||||
|
## Scope vs. the Python
|
||||||
|
|
||||||
|
The MATLAB port reproduces the Python's three core models (count, rate,
|
||||||
|
learning-rate), per-animal tests, unknown-group classification, and the H2
|
||||||
|
verdict (read off the main model's A2-vs-B2 term). The Python's two auxiliary
|
||||||
|
sensitivity sections — `report_anchor_only` (a two-anchor-only refit) and
|
||||||
|
`report_mixed` (a Bayesian MAP random-intercept model) — are intentionally
|
||||||
|
**not** ported: the anchor-only refit is closely tracked by the main-model H2
|
||||||
|
term, and MATLAB's native `fitglme` already *is* the random-intercept model.
|
||||||
@@ -0,0 +1,35 @@
|
|||||||
|
Data,Success
|
||||||
|
|
||||||
|
Group,—,—,Electrode-Box-A,Electrode-Box-A2,Electrode-Box-A2,Electrode-Box-A2,Electrode-Box-B2,Electrode-Box-B2,Electrode-Box-B2,Electrode-Box-B2,Left-Electrode-Cathodal,Naive,Naive,Naive,Right-Electrode
|
||||||
|
# Days Reach,Flan-1,Vu-vuong,Khoai-tay-1,Banh-mi-2,Egg-tart-2,Root-beer-2,Banh-mi-1,Egg-tart-1,Flan-2,Root-beer-1,OM-1,Khoai-lang-2,Khoai-tay-2,OM-2,Khoai-lang-1
|
||||||
|
-4,,,,,,,,,,,,,,,0
|
||||||
|
-3,,,,,,,,,,,,,,,2
|
||||||
|
-2,,,1,,,,,,,,,,,0,
|
||||||
|
-1,,,3,2,,,,,,,,,,5,0
|
||||||
|
0,,18,14,25,9,22,13,7,,11,0,0,,1,3
|
||||||
|
1,,35,22,22,2,31,35,16,,18,1,0,21,11,18
|
||||||
|
2,,61,6,22,31,49,47,23,,40,3,10,6,19,27
|
||||||
|
3,,34,28,29,44,31,65,63,,55,1,11,0,29,46
|
||||||
|
4,,61,60,27,54,60,70,69,,75,11,9,28,73,65
|
||||||
|
5,,37,75,43,84,84,102,83,,64,5,34,31,53,76
|
||||||
|
6,,49,81,74,85,,90,71,,104,9,21,62,73,79
|
||||||
|
7,,65,78,70,79,,109,79,,98,,23,87,80,81
|
||||||
|
8,,67,91,65,76,,104,98,,81,,64,105,91,98
|
||||||
|
9,,73,99,79,88,,119,89,,89,,75,75,90,101
|
||||||
|
10,,72,94,93,81,,121,96,,105,,63,101,95,103
|
||||||
|
11,,89,110,,78,,121,96,,,,59,102,60,113
|
||||||
|
12,,97,105,,96,,120,101,,,,51,95,58,117
|
||||||
|
13,,93,106,,84,,135,103,,,,73,77,,114
|
||||||
|
14,,75,91,,,,,97,,,,82,91,,119
|
||||||
|
15,,68,,,,,,,,,,70,,,111
|
||||||
|
16,,96,,,,,,,,,,76,,,74
|
||||||
|
17,,68,,,,,,,,,,76,,,109
|
||||||
|
18,,81,,,,,,,,,,63,,,88
|
||||||
|
19,,66,,,,,,,,,,48,,,98
|
||||||
|
20,,84,,,,,,,,,,65,,,118
|
||||||
|
21,,66,,,,,,,,,,75,,,109
|
||||||
|
22,,,,,,,,,,,,98,,,106
|
||||||
|
23,,,,,,,,,,,,88,,,
|
||||||
|
24,,,,,,,,,,,,94,,,
|
||||||
|
25,,,,,,,,,,,,56,,,
|
||||||
|
26,,,,,,,,,,,,75,,,
|
||||||
|
@@ -0,0 +1,35 @@
|
|||||||
|
Data,Total attempts
|
||||||
|
|
||||||
|
Group,—,—,Electrode-Box-A,Electrode-Box-A2,Electrode-Box-A2,Electrode-Box-A2,Electrode-Box-B2,Electrode-Box-B2,Electrode-Box-B2,Electrode-Box-B2,Left-Electrode-Cathodal,Naive,Naive,Naive,Right-Electrode
|
||||||
|
# Days Reach,Flan-1,Vu-vuong,Khoai-tay-1,Banh-mi-2,Egg-tart-2,Root-beer-2,Banh-mi-1,Egg-tart-1,Flan-2,Root-beer-1,OM-1,Khoai-lang-2,Khoai-tay-2,OM-2,Khoai-lang-1
|
||||||
|
-4,,,,,,,,,,,,,,,0
|
||||||
|
-3,,,,,,,,,20,,,,,,7
|
||||||
|
-2,,,18,,,,,,,,,,,0,
|
||||||
|
-1,,,24,42,,,,,,,,,,18,0
|
||||||
|
0,,61,62,97,32,74,84,56,,85,0,0,,7,69
|
||||||
|
1,,91,83,101,38,87,86,78,,76,10,0,68,33,97
|
||||||
|
2,,127,79,119,93,134,110,103,,105,5,47,79,62,84
|
||||||
|
3,,146,97,118,101,89,140,120,,134,2,52,83,117,121
|
||||||
|
4,,141,134,136,131,140,127,132,,136,44,56,87,126,143
|
||||||
|
5,,151,138,146,139,147,142,136,,133,23,95,125,135,141
|
||||||
|
6,,131,137,146,145,,131,142,,139,32,72,144,138,146
|
||||||
|
7,,149,147,148,143,,148,138,,148,,99,141,131,153
|
||||||
|
8,,141,132,130,131,,137,142,,145,,136,152,141,144
|
||||||
|
9,,140,146,151,149,,150,139,,156,,131,148,135,149
|
||||||
|
10,,131,143,152,151,,158,143,,158,,134,149,142,150
|
||||||
|
11,,158,156,,152,,148,148,,,,139,154,133,154
|
||||||
|
12,,154,143,,155,,149,156,,,,129,144,142,148
|
||||||
|
13,,145,153,,155,,154,152,,,,143,149,,146
|
||||||
|
14,,147,152,,,,,152,,,,136,153,,132
|
||||||
|
15,,142,,,,,,,,,,145,,,156
|
||||||
|
16,,162,,,,,,,,,,135,,,156
|
||||||
|
17,,144,,,,,,,,,,150,,,139
|
||||||
|
18,,134,,,,,,,,,,122,,,138
|
||||||
|
19,,144,,,,,,,,,,116,,,140
|
||||||
|
20,,127,,,,,,,,,,134,,,153
|
||||||
|
21,,133,,,,,,,,,,131,,,146
|
||||||
|
22,,,,,,,,,,,,146,,,141
|
||||||
|
23,,,,,,,,,,,,139,,,
|
||||||
|
24,,,,,,,,,,,,148,,,
|
||||||
|
25,,,,,,,,,,,,102,,,
|
||||||
|
26,,,,,,,,,,,,143,,,
|
||||||
|
@@ -0,0 +1,13 @@
|
|||||||
|
%RUN_ALL Run all 6 tDCS GLM scenarios and write results/<scenario>.txt for each.
|
||||||
|
% Usage: matlab -batch "run_all"
|
||||||
|
|
||||||
|
scenarios = {'unmerged_full', 'unmerged_d0_10', ...
|
||||||
|
'mergeA2_full', 'mergeA2_d0_10', ...
|
||||||
|
'mergeB2_full', 'mergeB2_d0_10'};
|
||||||
|
|
||||||
|
for i = 1:numel(scenarios)
|
||||||
|
fprintf('\n\n### Running scenario: %s ###\n', scenarios{i});
|
||||||
|
tdcs_glm(scenarios{i});
|
||||||
|
end
|
||||||
|
|
||||||
|
fprintf('\n\nAll %d scenarios complete. See analysis/matlab/results/*.txt.\n', numel(scenarios));
|
||||||
@@ -0,0 +1,14 @@
|
|||||||
|
%RUN_TESTS Run the full matlab.unittest suite and exit nonzero on failure.
|
||||||
|
% Usage: matlab -batch "run_tests"
|
||||||
|
% (add the analysis/matlab directory to the path first if it is not
|
||||||
|
% already on it, e.g. matlab -batch "addpath('analysis/matlab'); run_tests")
|
||||||
|
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
results = runtests(fullfile(thisDir, 'tests'));
|
||||||
|
disp(results);
|
||||||
|
|
||||||
|
fprintf('PASS=%d FAIL=%d\n', sum([results.Passed]), sum([results.Failed]));
|
||||||
|
|
||||||
|
if any([results.Failed])
|
||||||
|
error('run_tests:failed', 'One or more tests failed.');
|
||||||
|
end
|
||||||
@@ -0,0 +1,57 @@
|
|||||||
|
function cfg = tdcs_config()
|
||||||
|
%TDCS_CONFIG Curated group map and merge maps for the tDCS reaching study.
|
||||||
|
% CFG = TDCS_CONFIG() returns a struct with:
|
||||||
|
% .groups containers.Map, subject name (char) -> curated group name
|
||||||
|
% (char), for the 12 curated animals only. All other
|
||||||
|
% animals present in the raw export are not in this map
|
||||||
|
% and must be dropped by callers.
|
||||||
|
% .ref reference group used to code the GLME group factor,
|
||||||
|
% 'Electrode-Box-B2'.
|
||||||
|
% .mergeA2 containers.Map collapsing the two UNKNOWN conditions
|
||||||
|
% toward the weaker anchor: Electrode-Box-A ->
|
||||||
|
% Electrode-Box-A2, Right-Electrode -> Electrode-Box-B2.
|
||||||
|
% .mergeB2 containers.Map collapsing both UNKNOWN conditions onto
|
||||||
|
% the stronger anchor: Electrode-Box-A -> Electrode-Box-B2,
|
||||||
|
% Right-Electrode -> Electrode-Box-B2.
|
||||||
|
% .anchorLow 'Electrode-Box-A2', the weaker of the two known anchors.
|
||||||
|
% .anchorHigh 'Electrode-Box-B2', the stronger known anchor (== ref).
|
||||||
|
% .unknowns cell array of the two ungrouped conditions that get
|
||||||
|
% classified against the two anchors: {'Electrode-Box-A',
|
||||||
|
% 'Right-Electrode'}.
|
||||||
|
%
|
||||||
|
% See docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Global
|
||||||
|
% Constraints) for the source of these values.
|
||||||
|
|
||||||
|
% Curated 12-animal group map (by animal name).
|
||||||
|
subjectGroup = {
|
||||||
|
'Vu-vuong', 'Naive'
|
||||||
|
'Khoai-lang-2', 'Naive'
|
||||||
|
'Khoai-tay-2', 'Naive'
|
||||||
|
'OM-2', 'Naive'
|
||||||
|
'Khoai-tay-1', 'Electrode-Box-A'
|
||||||
|
'Banh-mi-2', 'Electrode-Box-A2'
|
||||||
|
'Egg-tart-2', 'Electrode-Box-A2'
|
||||||
|
'Root-beer-2', 'Electrode-Box-A2'
|
||||||
|
'Banh-mi-1', 'Electrode-Box-B2'
|
||||||
|
'Egg-tart-1', 'Electrode-Box-B2'
|
||||||
|
'Root-beer-1', 'Electrode-Box-B2'
|
||||||
|
'Khoai-lang-1', 'Right-Electrode'
|
||||||
|
};
|
||||||
|
|
||||||
|
cfg = struct();
|
||||||
|
cfg.groups = containers.Map(subjectGroup(:, 1), subjectGroup(:, 2));
|
||||||
|
|
||||||
|
cfg.ref = 'Electrode-Box-B2';
|
||||||
|
cfg.anchorLow = 'Electrode-Box-A2';
|
||||||
|
cfg.anchorHigh = 'Electrode-Box-B2';
|
||||||
|
cfg.unknowns = {'Electrode-Box-A', 'Right-Electrode'};
|
||||||
|
|
||||||
|
cfg.mergeA2 = containers.Map( ...
|
||||||
|
{'Electrode-Box-A', 'Right-Electrode'}, ...
|
||||||
|
{'Electrode-Box-A2', 'Electrode-Box-B2'});
|
||||||
|
|
||||||
|
cfg.mergeB2 = containers.Map( ...
|
||||||
|
{'Electrode-Box-A', 'Right-Electrode'}, ...
|
||||||
|
{'Electrode-Box-B2', 'Electrode-Box-B2'});
|
||||||
|
|
||||||
|
end
|
||||||
@@ -0,0 +1,64 @@
|
|||||||
|
function tdcs_glm(scenario)
|
||||||
|
%TDCS_GLM Run one tDCS GLM scenario end-to-end and report it.
|
||||||
|
% TDCS_GLM(SCENARIO) runs the full pipeline for one named scenario --
|
||||||
|
% TDCS_SCENARIO_DATA -> TDCS_MODELS -> TDCS_REPORT -- printing the
|
||||||
|
% report to the console and writing it to
|
||||||
|
% analysis/matlab/results/<SCENARIO>.txt.
|
||||||
|
%
|
||||||
|
% SCENARIO must be one of (case-sensitive):
|
||||||
|
% unmerged_full, unmerged_d0_10,
|
||||||
|
% mergeA2_full, mergeA2_d0_10,
|
||||||
|
% mergeB2_full, mergeB2_d0_10
|
||||||
|
%
|
||||||
|
% Each case below is a self-contained, independently invocable call of
|
||||||
|
% the same three-function pipeline; they are enumerated explicitly
|
||||||
|
% (rather than sharing one generic branch) so each scenario name is
|
||||||
|
% greppable on its own, per the project's spec.
|
||||||
|
%
|
||||||
|
% Example:
|
||||||
|
% tdcs_glm('mergeB2_full')
|
||||||
|
%
|
||||||
|
% See analysis/tdcs_glm.py (the Python this ports) and
|
||||||
|
% docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Task 4).
|
||||||
|
|
||||||
|
cfg = tdcs_config();
|
||||||
|
|
||||||
|
switch scenario
|
||||||
|
case 'unmerged_full'
|
||||||
|
[S, meta] = tdcs_scenario_data('unmerged_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'unmerged_full');
|
||||||
|
|
||||||
|
case 'unmerged_d0_10'
|
||||||
|
[S, meta] = tdcs_scenario_data('unmerged_d0_10');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'unmerged_d0_10');
|
||||||
|
|
||||||
|
case 'mergeA2_full'
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeA2_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'mergeA2_full');
|
||||||
|
|
||||||
|
case 'mergeA2_d0_10'
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeA2_d0_10');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'mergeA2_d0_10');
|
||||||
|
|
||||||
|
case 'mergeB2_full'
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeB2_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'mergeB2_full');
|
||||||
|
|
||||||
|
case 'mergeB2_d0_10'
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeB2_d0_10');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
tdcs_report(S, R, meta, cfg, 'mergeB2_d0_10');
|
||||||
|
|
||||||
|
otherwise
|
||||||
|
error('tdcs_glm:badScenario', ...
|
||||||
|
['Unrecognized scenario "%s". Valid scenarios are: ' ...
|
||||||
|
'unmerged_full, unmerged_d0_10, mergeA2_full, mergeA2_d0_10, ' ...
|
||||||
|
'mergeB2_full, mergeB2_d0_10.'], scenario);
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
function T = tdcs_load_data()
|
||||||
|
%TDCS_LOAD_DATA Load, join, and curate the raw tDCS reaching matrices.
|
||||||
|
% T = TDCS_LOAD_DATA() reads analysis/matlab/raw/raw_success.csv and
|
||||||
|
% raw_total.csv (app-exported matrix format: metric header line, blank
|
||||||
|
% line, file "Group" line, "# Days Reach" subject-header line, then one
|
||||||
|
% row per day with one value per subject column), reshapes each to long
|
||||||
|
% form, joins them on (subject, day), applies the curated 12-animal
|
||||||
|
% group map from tdcs_config (the file's own Group row is ignored),
|
||||||
|
% and restricts to day >= 0.
|
||||||
|
%
|
||||||
|
% Returns a table with variables:
|
||||||
|
% subject (string) animal name
|
||||||
|
% group (categorical) curated group (Naive, Electrode-Box-A,
|
||||||
|
% Electrode-Box-A2, Electrode-Box-B2,
|
||||||
|
% Right-Electrode)
|
||||||
|
% day (double) "# Days Reach" value, >= 0
|
||||||
|
% success (double) successful reaches in the session
|
||||||
|
% total (double) total attempts in the session
|
||||||
|
%
|
||||||
|
% Only the 12 curated animals (per tdcs_config().groups) are retained;
|
||||||
|
% any other animal present in the raw export is dropped.
|
||||||
|
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
successFile = fullfile(thisDir, 'raw', 'raw_success.csv');
|
||||||
|
totalFile = fullfile(thisDir, 'raw', 'raw_total.csv');
|
||||||
|
|
||||||
|
Tsuccess = localReadMatrix(successFile, 'success');
|
||||||
|
Ttotal = localReadMatrix(totalFile, 'total');
|
||||||
|
|
||||||
|
% Left join keyed on (subject, day): every non-blank Success cell becomes
|
||||||
|
% a row; Total is looked up for the same (subject, day). Per the plan's
|
||||||
|
% Global Constraints, Total is fully populated wherever Success is, so
|
||||||
|
% this should never introduce a missing total — but the join is a left
|
||||||
|
% join (not an inner join) specifically so that an unmatched Success row
|
||||||
|
% would still be kept (with total = NaN) rather than silently dropped.
|
||||||
|
T = outerjoin(Tsuccess, Ttotal, 'Keys', {'subject', 'day'}, ...
|
||||||
|
'MergeKeys', true, 'Type', 'left');
|
||||||
|
|
||||||
|
% Apply the curated map: keep only the 12 curated animals, and assign
|
||||||
|
% group from the map (overriding the file's own Group row).
|
||||||
|
cfg = tdcs_config();
|
||||||
|
subjectsCell = cellstr(T.subject);
|
||||||
|
keep = isKey(cfg.groups, subjectsCell);
|
||||||
|
T = T(keep, :);
|
||||||
|
subjectsCell = subjectsCell(keep);
|
||||||
|
|
||||||
|
groupVals = cell(height(T), 1);
|
||||||
|
for i = 1:height(T)
|
||||||
|
groupVals{i} = cfg.groups(subjectsCell{i});
|
||||||
|
end
|
||||||
|
T.group = categorical(groupVals);
|
||||||
|
|
||||||
|
% Restrict to day >= 0 (analysis window).
|
||||||
|
T = T(T.day >= 0, :);
|
||||||
|
|
||||||
|
% Final variable order and a deterministic row order.
|
||||||
|
T = T(:, {'subject', 'group', 'day', 'success', 'total'});
|
||||||
|
T = sortrows(T, {'subject', 'day'});
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
function Tlong = localReadMatrix(filePath, valueName)
|
||||||
|
%LOCALREADMATRIX Parse one app-exported day-by-subject matrix to long form.
|
||||||
|
% Line 1 = "Data,<metric>" (skipped), line 2 blank (skipped), line 3 =
|
||||||
|
% "Group,<...>" (skipped -- grouping comes from the curated map, not
|
||||||
|
% this row), line 4 = "# Days Reach,<subject headers>" (columns 2:end
|
||||||
|
% are the subject names), lines 5+ = "<day>,<value per subject>" with
|
||||||
|
% blank cells meaning "no session that day for that subject".
|
||||||
|
|
||||||
|
lines = readlines(filePath);
|
||||||
|
% Defensive: drop any wholly-blank trailing line (e.g. if the file ends
|
||||||
|
% with a newline). Real data rows always start with a day number and are
|
||||||
|
% never zero-length.
|
||||||
|
while numel(lines) > 0 && strlength(strip(lines(end))) == 0
|
||||||
|
lines(end) = [];
|
||||||
|
end
|
||||||
|
|
||||||
|
headerParts = strsplit(lines(4), ',', 'CollapseDelimiters', false);
|
||||||
|
subjects = string(headerParts(2:end));
|
||||||
|
nSubjects = numel(subjects);
|
||||||
|
|
||||||
|
dataLines = lines(5:end);
|
||||||
|
dataLines = dataLines(strlength(strip(dataLines)) > 0);
|
||||||
|
nRows = numel(dataLines);
|
||||||
|
|
||||||
|
day = nan(nRows * nSubjects, 1);
|
||||||
|
subjectCol = strings(nRows * nSubjects, 1);
|
||||||
|
value = nan(nRows * nSubjects, 1);
|
||||||
|
|
||||||
|
idx = 0;
|
||||||
|
for r = 1:nRows
|
||||||
|
parts = strsplit(dataLines(r), ',', 'CollapseDelimiters', false);
|
||||||
|
dayVal = str2double(parts(1));
|
||||||
|
for c = 1:nSubjects
|
||||||
|
col = c + 1;
|
||||||
|
if col <= numel(parts)
|
||||||
|
cellStr = strtrim(parts(col));
|
||||||
|
else
|
||||||
|
cellStr = "";
|
||||||
|
end
|
||||||
|
if strlength(cellStr) == 0
|
||||||
|
continue
|
||||||
|
end
|
||||||
|
idx = idx + 1;
|
||||||
|
day(idx) = dayVal;
|
||||||
|
subjectCol(idx) = subjects(c);
|
||||||
|
value(idx) = str2double(cellStr);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
day = day(1:idx);
|
||||||
|
subjectCol = subjectCol(1:idx);
|
||||||
|
value = value(1:idx);
|
||||||
|
|
||||||
|
Tlong = table(subjectCol, day, value, 'VariableNames', {'subject', 'day', valueName});
|
||||||
|
|
||||||
|
end
|
||||||
@@ -0,0 +1,216 @@
|
|||||||
|
function R = tdcs_models(S, meta, cfg)
|
||||||
|
%TDCS_MODELS Fit the tDCS GLMMs, per-animal tests, and (unmerged) classification.
|
||||||
|
% R = TDCS_MODELS(S, META, CFG) consumes one scenario table S and its
|
||||||
|
% META (from TDCS_SCENARIO_DATA) plus CFG (from TDCS_CONFIG), and
|
||||||
|
% returns a struct R:
|
||||||
|
%
|
||||||
|
% R.countGLME GeneralizedLinearMixedModel, Poisson count-level model:
|
||||||
|
% success ~ group + day_c + day_c2 + (1|subject)
|
||||||
|
% R.rateGLME GeneralizedLinearMixedModel, Binomial rate-level model
|
||||||
|
% (total==0 sessions dropped first):
|
||||||
|
% success ~ group + day_c + day_c2 + (1|subject)
|
||||||
|
% R.learnGLME GeneralizedLinearMixedModel, Poisson learning-rate
|
||||||
|
% model with a group x day_c interaction:
|
||||||
|
% success ~ group*day_c + day_c2 + (1|subject)
|
||||||
|
% R.learn struct with the joint Wald test (via COEFTEST) that
|
||||||
|
% all group x day_c interaction terms are zero:
|
||||||
|
% .jointP, .jointF, .jointDF1, .jointDF2.
|
||||||
|
% R.perAnimal per-subject Box-B2 vs Box-A2 comparison (pure stats,
|
||||||
|
% no GLME involved): .countWelchP, .countMWUp,
|
||||||
|
% .rateWelchP, .rateMWUp, .nB, .nA. "count" = per-subject
|
||||||
|
% mean success; "rate" = per-subject pooled
|
||||||
|
% sum(success)/sum(total). Welch = TTEST2(...,
|
||||||
|
% 'Vartype','unequal') (two-sided); MWU = RANKSUM(...,
|
||||||
|
% 'tail','right') (one-sided, B2 > A2).
|
||||||
|
% R.h2 anchor check (Box-A2 vs Box-B2) read off the
|
||||||
|
% count/rate GLMEs' group_<anchorLow> fixed effect
|
||||||
|
% (coded vs the cfg.ref = anchorHigh reference level):
|
||||||
|
% .countIRR, .countOneSidedP, .rateOR, .rateOneSidedP.
|
||||||
|
% IRR/OR = 1/exp(coef) is the Box-B2-over-Box-A2
|
||||||
|
% advantage (so > 1 supports H2); one-sided p =
|
||||||
|
% NORMCDF(coef/se), testing coef < 0 (B2 above A2).
|
||||||
|
% R.classify (unmerged scenarios only, i.e. meta.isUnmerged) struct
|
||||||
|
% with one field per unknown in cfg.unknowns (field name
|
||||||
|
% via matlab.lang.makeValidName), each holding .count and
|
||||||
|
% .rate sub-structs with .vsA2 (.ratio, .p), .vsB2
|
||||||
|
% (.ratio, .p), and .label (one of 'B2-like', 'A2-like',
|
||||||
|
% 'AMBIGUOUS (nearer <anchor>)', 'UNLIKE BOTH'), per
|
||||||
|
% METHODS.md's classification logic. For merged
|
||||||
|
% scenarios (meta.isUnmerged == false), R.classify is an
|
||||||
|
% empty (no-field) struct.
|
||||||
|
%
|
||||||
|
% See analysis/tdcs_glm.py (the Python this ports) and
|
||||||
|
% docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Task 3).
|
||||||
|
|
||||||
|
countFormula = 'success ~ group + day_c + day_c2 + (1|subject)';
|
||||||
|
R.countGLME = fitglme(S, countFormula, 'Distribution', 'Poisson', 'Link', 'log');
|
||||||
|
|
||||||
|
Sr = S(S.total > 0, :);
|
||||||
|
R.rateGLME = fitglme(Sr, countFormula, 'Distribution', 'Binomial', ...
|
||||||
|
'BinomialSize', Sr.total, 'Link', 'logit');
|
||||||
|
|
||||||
|
learnFormula = 'success ~ group*day_c + day_c2 + (1|subject)';
|
||||||
|
R.learnGLME = fitglme(S, learnFormula, 'Distribution', 'Poisson');
|
||||||
|
R.learn = localJointInteractionTest(R.learnGLME);
|
||||||
|
|
||||||
|
R.perAnimal = localPerAnimal(S, cfg);
|
||||||
|
R.h2 = localH2(R.countGLME, R.rateGLME, cfg);
|
||||||
|
|
||||||
|
if meta.isUnmerged
|
||||||
|
R.classify = localClassifyAll(R.countGLME, R.rateGLME, cfg);
|
||||||
|
else
|
||||||
|
R.classify = struct();
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
% --------------------------------------------------------------- per-animal
|
||||||
|
function out = localPerAnimal(S, cfg)
|
||||||
|
%LOCALPERANIMAL Per-subject Box-B2 vs Box-A2 Welch t-test / MWU (pure stats).
|
||||||
|
b = S(S.group == cfg.anchorHigh, :);
|
||||||
|
a = S(S.group == cfg.anchorLow, :);
|
||||||
|
|
||||||
|
[bCount, bRate, nB] = localSubjectSummary(b);
|
||||||
|
[aCount, aRate, nA] = localSubjectSummary(a);
|
||||||
|
|
||||||
|
[~, countWelchP] = ttest2(bCount, aCount, 'Vartype', 'unequal');
|
||||||
|
countMWUp = ranksum(bCount, aCount, 'tail', 'right');
|
||||||
|
[~, rateWelchP] = ttest2(bRate, aRate, 'Vartype', 'unequal');
|
||||||
|
rateMWUp = ranksum(bRate, aRate, 'tail', 'right');
|
||||||
|
|
||||||
|
out = struct( ...
|
||||||
|
'countWelchP', countWelchP, 'countMWUp', countMWUp, ...
|
||||||
|
'rateWelchP', rateWelchP, 'rateMWUp', rateMWUp, ...
|
||||||
|
'nB', nB, 'nA', nA);
|
||||||
|
end
|
||||||
|
|
||||||
|
function [meanSuccess, pooledRate, nSubj] = localSubjectSummary(T)
|
||||||
|
%LOCALSUBJECTSUMMARY Per-subject mean success and pooled success/total.
|
||||||
|
subj = unique(T.subject);
|
||||||
|
nSubj = numel(subj);
|
||||||
|
meanSuccess = zeros(nSubj, 1);
|
||||||
|
pooledRate = zeros(nSubj, 1);
|
||||||
|
for i = 1:nSubj
|
||||||
|
rows = T.subject == subj(i);
|
||||||
|
meanSuccess(i) = mean(T.success(rows));
|
||||||
|
pooledRate(i) = sum(T.success(rows)) / sum(T.total(rows));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% --------------------------------------------------------------------- H2
|
||||||
|
function h2 = localH2(countGLME, rateGLME, cfg)
|
||||||
|
%LOCALH2 Anchor check: Box-A2 vs Box-B2, read off each GLME's group term.
|
||||||
|
[countCoef, countSE] = localCoefAndSE(countGLME, cfg.anchorLow, cfg.ref);
|
||||||
|
[rateCoef, rateSE] = localCoefAndSE(rateGLME, cfg.anchorLow, cfg.ref);
|
||||||
|
|
||||||
|
h2 = struct( ...
|
||||||
|
'countIRR', 1 / exp(countCoef), ...
|
||||||
|
'countOneSidedP', normcdf(countCoef / countSE), ...
|
||||||
|
'rateOR', 1 / exp(rateCoef), ...
|
||||||
|
'rateOneSidedP', normcdf(rateCoef / rateSE));
|
||||||
|
end
|
||||||
|
|
||||||
|
% -------------------------------------------------------------- classify
|
||||||
|
function out = localClassifyAll(countGLME, rateGLME, cfg)
|
||||||
|
%LOCALCLASSIFYALL Classify every unknown in cfg.unknowns vs both anchors.
|
||||||
|
out = struct();
|
||||||
|
for i = 1:numel(cfg.unknowns)
|
||||||
|
unknown = cfg.unknowns{i};
|
||||||
|
fieldName = matlab.lang.makeValidName(unknown);
|
||||||
|
out.(fieldName) = struct( ...
|
||||||
|
'count', localClassifyOne(countGLME, unknown, cfg), ...
|
||||||
|
'rate', localClassifyOne(rateGLME, unknown, cfg));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function verdict = localClassifyOne(model, unknown, cfg)
|
||||||
|
%LOCALCLASSIFYONE Compare UNKNOWN vs both anchors on one GLME; label it.
|
||||||
|
% Mirrors classify_one() in analysis/tdcs_glm.py / METHODS.md's
|
||||||
|
% "Classification logic (unknowns)": indistinguishable (p >= 0.05) from
|
||||||
|
% one anchor and different from the other => that anchor's label;
|
||||||
|
% indistinguishable from both => AMBIGUOUS (report the numerically
|
||||||
|
% nearer one); different from both => UNLIKE BOTH.
|
||||||
|
[coefLo, pLo] = localContrast(model, unknown, cfg.anchorLow, cfg.ref);
|
||||||
|
[coefHi, pHi] = localContrast(model, unknown, cfg.anchorHigh, cfg.ref);
|
||||||
|
|
||||||
|
vsA2 = struct('ratio', exp(coefLo), 'p', pLo);
|
||||||
|
vsB2 = struct('ratio', exp(coefHi), 'p', pHi);
|
||||||
|
|
||||||
|
likeLo = pLo >= 0.05;
|
||||||
|
likeHi = pHi >= 0.05;
|
||||||
|
if likeHi && ~likeLo
|
||||||
|
label = 'B2-like';
|
||||||
|
elseif likeLo && ~likeHi
|
||||||
|
label = 'A2-like';
|
||||||
|
elseif likeLo && likeHi
|
||||||
|
if abs(coefHi) < abs(coefLo)
|
||||||
|
nearer = cfg.anchorHigh;
|
||||||
|
else
|
||||||
|
nearer = cfg.anchorLow;
|
||||||
|
end
|
||||||
|
label = sprintf('AMBIGUOUS (nearer %s)', nearer);
|
||||||
|
else
|
||||||
|
label = 'UNLIKE BOTH';
|
||||||
|
end
|
||||||
|
|
||||||
|
verdict = struct('vsA2', vsA2, 'vsB2', vsB2, 'label', label);
|
||||||
|
end
|
||||||
|
|
||||||
|
% ------------------------------------------------------------- contrasts
|
||||||
|
function [coef, se] = localCoefAndSE(model, aName, refName)
|
||||||
|
%LOCALCOEFANDSE Coefficient + SE for group level A_NAME vs the reference.
|
||||||
|
v = localContrastVector(model, aName, refName, refName);
|
||||||
|
beta = model.Coefficients.Estimate;
|
||||||
|
Cov = model.CoefficientCovariance;
|
||||||
|
coef = v * beta;
|
||||||
|
se = sqrt(v * Cov * v');
|
||||||
|
end
|
||||||
|
|
||||||
|
function [coef, p] = localContrast(model, aName, bName, refName)
|
||||||
|
%LOCALCONTRAST Linear contrast A_NAME vs B_NAME (both coded vs REF_NAME).
|
||||||
|
% Two-sided p via COEFTEST (an F-test on one row = the two-sided
|
||||||
|
% Wald/t-test); works for any pair, including B_NAME == REF_NAME.
|
||||||
|
v = localContrastVector(model, aName, bName, refName);
|
||||||
|
beta = model.Coefficients.Estimate;
|
||||||
|
coef = v * beta;
|
||||||
|
p = coefTest(model, v);
|
||||||
|
end
|
||||||
|
|
||||||
|
function v = localContrastVector(model, aName, bName, refName)
|
||||||
|
%LOCALCONTRASTVECTOR +1 at group_A_NAME (if not the reference), -1 at
|
||||||
|
% group_B_NAME (if not the reference), over MODEL.CoefficientNames.
|
||||||
|
cn = model.CoefficientNames;
|
||||||
|
v = zeros(1, numel(cn));
|
||||||
|
if ~strcmp(aName, refName)
|
||||||
|
v(localCoefIndex(cn, aName)) = v(localCoefIndex(cn, aName)) + 1;
|
||||||
|
end
|
||||||
|
if ~strcmp(bName, refName)
|
||||||
|
v(localCoefIndex(cn, bName)) = v(localCoefIndex(cn, bName)) - 1;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function idx = localCoefIndex(coefNames, groupName)
|
||||||
|
%LOCALCOEFINDEX Index into CoefficientNames for a "group_<groupName>" term.
|
||||||
|
idx = find(strcmp(coefNames, ['group_' groupName]), 1);
|
||||||
|
if isempty(idx)
|
||||||
|
error('tdcs_models:missingGroupTerm', ...
|
||||||
|
'No "group_%s" coefficient in this GLME (CoefficientNames: %s).', ...
|
||||||
|
groupName, strjoin(coefNames, ', '));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% -------------------------------------------------------------- learning
|
||||||
|
function out = localJointInteractionTest(learnGLME)
|
||||||
|
%LOCALJOINTINTERACTIONTEST Joint Wald test that all group x day_c
|
||||||
|
% interaction fixed effects are zero (do groups learn at different
|
||||||
|
% rates?), mirroring the GEE wald_test in analysis/tdcs_glm.py.
|
||||||
|
cn = learnGLME.CoefficientNames;
|
||||||
|
isInteraction = contains(cn, ':day_c') & ~contains(cn, 'day_c2');
|
||||||
|
idx = find(isInteraction);
|
||||||
|
H = zeros(numel(idx), numel(cn));
|
||||||
|
for i = 1:numel(idx)
|
||||||
|
H(i, idx(i)) = 1;
|
||||||
|
end
|
||||||
|
[p, F, df1, df2] = coefTest(learnGLME, H);
|
||||||
|
out = struct('jointP', p, 'jointF', F, 'jointDF1', df1, 'jointDF2', df2);
|
||||||
|
end
|
||||||
@@ -0,0 +1,172 @@
|
|||||||
|
function tdcs_report(S, R, meta, cfg, scenario)
|
||||||
|
%TDCS_REPORT Print + persist the tDCS GLM report for one scenario.
|
||||||
|
% TDCS_REPORT(S, R, META, CFG, SCENARIO) prints a three-part report to
|
||||||
|
% the console (DESCRIPTIVES, RAW GLM, INTERPRETATION) and writes the
|
||||||
|
% same text to analysis/matlab/results/<SCENARIO>.txt (the results/
|
||||||
|
% folder is created if it does not already exist).
|
||||||
|
%
|
||||||
|
% S scenario table, from TDCS_SCENARIO_DATA.
|
||||||
|
% R model results struct, from TDCS_MODELS.
|
||||||
|
% META scenario metadata struct, from TDCS_SCENARIO_DATA
|
||||||
|
% (.mergeKey, .maxDay, .isUnmerged).
|
||||||
|
% CFG config struct, from TDCS_CONFIG.
|
||||||
|
% SCENARIO scenario name (char/string), used only for the report
|
||||||
|
% header and the output filename.
|
||||||
|
%
|
||||||
|
% The INTERPRETATION section mirrors verdicts()/classify() in
|
||||||
|
% analysis/tdcs_glm.py and the "Classification logic (unknowns)" /
|
||||||
|
% "Caveats" sections of analysis/METHODS.md: the H2 anchor-check
|
||||||
|
% verdict (IRR/OR + one-sided p + SUPPORTED/not supported), the
|
||||||
|
% per-animal Welch/Mann-Whitney p-values, the unknown-condition
|
||||||
|
% classification (unmerged scenarios only), and the standing caveats
|
||||||
|
% (tiny groups, single-subject classification, count-vs-rate).
|
||||||
|
%
|
||||||
|
% See analysis/tdcs_glm.py, analysis/METHODS.md, and
|
||||||
|
% docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Task 4).
|
||||||
|
|
||||||
|
scenario = char(scenario);
|
||||||
|
txt = '';
|
||||||
|
|
||||||
|
txt = [txt localHeader(sprintf('tDCS GLM report -- scenario: %s', scenario))];
|
||||||
|
txt = [txt sprintf('merge key: %s day window: 0..%d observations: %d\n', ...
|
||||||
|
meta.mergeKey, meta.maxDay, height(S))];
|
||||||
|
|
||||||
|
txt = [txt localDescriptives(S)];
|
||||||
|
txt = [txt localRawGLM(R)];
|
||||||
|
txt = [txt localInterpretation(R, meta, cfg)];
|
||||||
|
|
||||||
|
fprintf('%s', txt);
|
||||||
|
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
resultsDir = fullfile(thisDir, 'results');
|
||||||
|
if ~exist(resultsDir, 'dir')
|
||||||
|
mkdir(resultsDir);
|
||||||
|
end
|
||||||
|
outFile = fullfile(resultsDir, [scenario '.txt']);
|
||||||
|
fid = fopen(outFile, 'w');
|
||||||
|
if fid == -1
|
||||||
|
error('tdcs_report:cannotWrite', 'Could not open "%s" for writing.', outFile);
|
||||||
|
end
|
||||||
|
cleanupObj = onCleanup(@() fclose(fid)); %#ok<NASGU>
|
||||||
|
fprintf(fid, '%s', txt);
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
% ------------------------------------------------------------------ header
|
||||||
|
function s = localHeader(title)
|
||||||
|
bar = repmat('=', 1, 78);
|
||||||
|
s = sprintf('%s\n%s\n%s\n', bar, title, bar);
|
||||||
|
end
|
||||||
|
|
||||||
|
% -------------------------------------------------------------- descriptives
|
||||||
|
function s = localDescriptives(S)
|
||||||
|
s = localHeader('DESCRIPTIVES');
|
||||||
|
h = sprintf('%-20s %8s %10s %13s %10s %8s', ...
|
||||||
|
'group', 'n_subj', 'n_sessions', 'mean_success', 'mean_rate', 'max_day');
|
||||||
|
s = [s h sprintf('\n') repmat('-', 1, length(h)) sprintf('\n')];
|
||||||
|
|
||||||
|
grps = categories(S.group);
|
||||||
|
for i = 1:numel(grps)
|
||||||
|
rows = S.group == grps{i};
|
||||||
|
nSubj = numel(unique(S.subject(rows)));
|
||||||
|
nSessions = sum(rows);
|
||||||
|
meanSuccess = mean(S.success(rows));
|
||||||
|
meanRate = sum(S.success(rows)) / sum(S.total(rows));
|
||||||
|
maxDay = max(S.day(rows));
|
||||||
|
s = [s sprintf('%-20s %8d %10d %13.1f %10.3f %8d\n', ...
|
||||||
|
grps{i}, nSubj, nSessions, meanSuccess, meanRate, maxDay)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
s = [s sprintf('\n')];
|
||||||
|
end
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- raw GLM
|
||||||
|
function s = localRawGLM(R)
|
||||||
|
s = localHeader('(A) LEVEL / COUNT -- Poisson GLMM (subject random intercept)');
|
||||||
|
s = [s localCleanDisp(R.countGLME) sprintf('\n')];
|
||||||
|
|
||||||
|
s = [s localHeader('(B) LEVEL / RATE -- Binomial GLMM (subject random intercept)')];
|
||||||
|
s = [s localCleanDisp(R.rateGLME) sprintf('\n')];
|
||||||
|
|
||||||
|
s = [s localHeader('(C) LEARNING RATE -- joint test (all group x day_c interactions = 0)')];
|
||||||
|
s = [s sprintf('F=%.3f, df1=%d, df2=%.2f, p=%.4f\n', ...
|
||||||
|
R.learn.jointF, R.learn.jointDF1, R.learn.jointDF2, R.learn.jointP)];
|
||||||
|
s = [s sprintf([' -> small p = groups improve at DIFFERENT rates; ' ...
|
||||||
|
'large p = parallel learning.\n\n'])];
|
||||||
|
end
|
||||||
|
|
||||||
|
function s = localCleanDisp(model)
|
||||||
|
%LOCALCLEANDISP disp(MODEL) captured to text, with the Command-Window-only
|
||||||
|
% <strong>...</strong> bold-markup tags (rendered by the interactive
|
||||||
|
% MATLAB terminal, but emitted as literal text by EVALC) stripped so the
|
||||||
|
% report/console/file output is plain text.
|
||||||
|
raw = evalc('disp(model)');
|
||||||
|
s = regexprep(raw, '</?strong>', '');
|
||||||
|
end
|
||||||
|
|
||||||
|
% ------------------------------------------------------------ interpretation
|
||||||
|
function s = localInterpretation(R, meta, cfg)
|
||||||
|
s = localHeader('INTERPRETATION');
|
||||||
|
s = [s sprintf(['Note: these are subject-level GLMMs (Laplace-approximated ' ...
|
||||||
|
'fitglme), not the\nPython reference''s population-average GEE -- ' ...
|
||||||
|
'directions/magnitudes should agree,\nexact ratios and p-values will ' ...
|
||||||
|
'differ.\n\n'])];
|
||||||
|
|
||||||
|
s = [s sprintf('Anchor check -- H2: Box-B2 BETTER than Box-A2 (the two anchors must differ)\n')];
|
||||||
|
s = [s localH2Line('count/level', R.h2.countIRR, R.h2.countOneSidedP)];
|
||||||
|
s = [s localH2Line('rate/level ', R.h2.rateOR, R.h2.rateOneSidedP)];
|
||||||
|
|
||||||
|
s = [s sprintf('\nPer-animal (Box-B2 n=%d vs Box-A2 n=%d), pure stats (no GLME):\n', ...
|
||||||
|
R.perAnimal.nB, R.perAnimal.nA)];
|
||||||
|
s = [s sprintf(' count (per-subject mean success): Welch two-sided p=%.4f, Mann-Whitney one-sided (B2>A2) p=%.4f\n', ...
|
||||||
|
R.perAnimal.countWelchP, R.perAnimal.countMWUp)];
|
||||||
|
s = [s sprintf(' rate (per-subject pooled success/total): Welch two-sided p=%.4f, Mann-Whitney one-sided (B2>A2) p=%.4f\n', ...
|
||||||
|
R.perAnimal.rateWelchP, R.perAnimal.rateMWUp)];
|
||||||
|
|
||||||
|
if meta.isUnmerged
|
||||||
|
s = [s sprintf('\nClassification -- is each UNKNOWN condition A2-like or B2-like?\n')];
|
||||||
|
s = [s sprintf('(ratio >1 = above that anchor; p = differs from that anchor)\n')];
|
||||||
|
for i = 1:numel(cfg.unknowns)
|
||||||
|
unknown = cfg.unknowns{i};
|
||||||
|
fieldName = matlab.lang.makeValidName(unknown);
|
||||||
|
entry = R.classify.(fieldName);
|
||||||
|
s = [s sprintf('\n %s:\n', unknown)]; %#ok<AGROW>
|
||||||
|
s = [s localClassifyLine('count/level', entry.count)]; %#ok<AGROW>
|
||||||
|
s = [s localClassifyLine('rate/level ', entry.rate)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
s = [s sprintf(['\n NOTE: each unknown has ONLY 1 subject. ''Matches'' means ' ...
|
||||||
|
'''not statistically\n distinguishable'', weak evidence at n=1, not proof ' ...
|
||||||
|
'of equivalence.\n'])];
|
||||||
|
else
|
||||||
|
s = [s sprintf(['\n(No unknown groups -- they were merged into the anchors; ' ...
|
||||||
|
'only the H2 anchor\ncontrast applies.)\n'])];
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf('\n') localHeader('CAVEATS')];
|
||||||
|
s = [s sprintf([' - Tiny groups: each arm has only n=3-4 subjects (Naive n=4; ' ...
|
||||||
|
'anchor arms n=3-5\n depending on merge), and -- in unmerged scenarios -- ' ...
|
||||||
|
'each unknown condition\n (Electrode-Box-A, Right-Electrode) has only ' ...
|
||||||
|
'n=1 subject. Treat every group\n comparison here as preliminary.\n'])];
|
||||||
|
s = [s sprintf([' - Single-subject classification: for a 1-subject unknown, ' ...
|
||||||
|
'''matches anchor X''\n means ''not statistically distinguishable from X'', ' ...
|
||||||
|
'NOT proof of equivalence;\n inference with a single subject in a group ' ...
|
||||||
|
'is fragile.\n'])];
|
||||||
|
s = [s sprintf([' - Count vs rate: ''success'' alone is a raw count; the rate ' ...
|
||||||
|
'model\n (success/attempts) is the fairer accuracy comparison when attempt ' ...
|
||||||
|
'counts differ\n between groups.\n'])];
|
||||||
|
end
|
||||||
|
|
||||||
|
function s = localH2Line(label, ratio, p)
|
||||||
|
tag = 'not supported';
|
||||||
|
if ratio > 1 && p < 0.05
|
||||||
|
tag = 'SUPPORTED';
|
||||||
|
end
|
||||||
|
s = sprintf(' [%s] Box-B2 = %.2fx Box-A2 (one-sided p=%.4f) -> %s\n', ...
|
||||||
|
label, ratio, p, tag);
|
||||||
|
end
|
||||||
|
|
||||||
|
function s = localClassifyLine(label, verdict)
|
||||||
|
s = sprintf([' [%s] vs Box-A2: %.2fx p=%.3f vs Box-B2: %.2fx p=%.3f\n' ...
|
||||||
|
' -> %s\n'], ...
|
||||||
|
label, verdict.vsA2.ratio, verdict.vsA2.p, ...
|
||||||
|
verdict.vsB2.ratio, verdict.vsB2.p, verdict.label);
|
||||||
|
end
|
||||||
@@ -0,0 +1,113 @@
|
|||||||
|
function [S, meta] = tdcs_scenario_data(scenario)
|
||||||
|
%TDCS_SCENARIO_DATA Build one merge/window scenario from the curated tDCS data.
|
||||||
|
% [S, META] = TDCS_SCENARIO_DATA(SCENARIO) loads the curated 12-animal
|
||||||
|
% table via TDCS_LOAD_DATA, applies the requested group-merge map from
|
||||||
|
% TDCS_CONFIG, restricts to day <= maxDay per the scenario's window
|
||||||
|
% suffix, adds mean-centered day covariates DAY_C / DAY_C2, and recodes
|
||||||
|
% GROUP as a categorical with CFG.REF ('Electrode-Box-B2') as the FIRST
|
||||||
|
% (reference) category via REORDERCATS. S is also written to
|
||||||
|
% analysis/matlab/derived/<SCENARIO>.csv via WRITETABLE (the derived/
|
||||||
|
% folder is created if it does not already exist).
|
||||||
|
%
|
||||||
|
% SCENARIO must be one of:
|
||||||
|
% unmerged_full, unmerged_d0_10, mergeA2_full, mergeA2_d0_10,
|
||||||
|
% mergeB2_full, mergeB2_d0_10
|
||||||
|
%
|
||||||
|
% The scenario name is split into:
|
||||||
|
% - a merge key: 'unmerged' | 'mergeA2' | 'mergeB2', selecting which
|
||||||
|
% containers.Map from TDCS_CONFIG (if any) to apply to relabel
|
||||||
|
% GROUP. 'unmerged' applies no relabeling (identity) -- it is not
|
||||||
|
% backed by a config map at all, since only the two "unknown"
|
||||||
|
% conditions (Electrode-Box-A, Right-Electrode) are ever merge
|
||||||
|
% targets, and the unmerged scenario reports them unmerged/as-is.
|
||||||
|
% - a day window: '_full' -> maxDay = 26, '_d0_10' -> maxDay = 10.
|
||||||
|
%
|
||||||
|
% META is a struct with fields:
|
||||||
|
% .mergeKey 'unmerged' | 'mergeA2' | 'mergeB2'
|
||||||
|
% .maxDay 26 or 10
|
||||||
|
% .isUnmerged true only when mergeKey == 'unmerged'
|
||||||
|
|
||||||
|
[mergeKey, maxDay] = localParseScenario(scenario);
|
||||||
|
|
||||||
|
cfg = tdcs_config();
|
||||||
|
T = tdcs_load_data();
|
||||||
|
|
||||||
|
switch mergeKey
|
||||||
|
case 'unmerged'
|
||||||
|
% Identity: no relabeling. An empty map means isKey(...) is
|
||||||
|
% always false below, so every group passes through unchanged.
|
||||||
|
mergeMap = containers.Map('KeyType', 'char', 'ValueType', 'char');
|
||||||
|
case 'mergeA2'
|
||||||
|
mergeMap = cfg.mergeA2;
|
||||||
|
case 'mergeB2'
|
||||||
|
mergeMap = cfg.mergeB2;
|
||||||
|
otherwise
|
||||||
|
error('tdcs_scenario_data:badScenario', ...
|
||||||
|
'Unrecognized scenario "%s".', scenario);
|
||||||
|
end
|
||||||
|
|
||||||
|
% Relabel GROUP per the merge map. Groups that are not a key in the map
|
||||||
|
% (e.g. Naive, Electrode-Box-A2, Electrode-Box-B2 under mergeA2/mergeB2,
|
||||||
|
% or every group under the empty 'unmerged' map) are left unchanged --
|
||||||
|
% only the two "unknown" conditions are ever merge-map keys.
|
||||||
|
groupCell = cellstr(T.group);
|
||||||
|
for i = 1:numel(groupCell)
|
||||||
|
if isKey(mergeMap, groupCell{i})
|
||||||
|
groupCell{i} = mergeMap(groupCell{i});
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
% Day window (applied after relabeling, before centering, per the plan).
|
||||||
|
dayMask = T.day <= maxDay;
|
||||||
|
S = T(dayMask, :);
|
||||||
|
groupCell = groupCell(dayMask);
|
||||||
|
|
||||||
|
% Recode GROUP as categorical with cfg.ref first (the GLME reference
|
||||||
|
% level). Only categories actually present after merge+window survive.
|
||||||
|
S.group = categorical(groupCell);
|
||||||
|
allCats = categories(S.group);
|
||||||
|
assert(ismember(cfg.ref, allCats), 'tdcs_scenario_data:refMissing', ...
|
||||||
|
'Reference group "%s" is absent from scenario "%s".', cfg.ref, char(scenario));
|
||||||
|
otherCats = allCats(~strcmp(allCats, cfg.ref));
|
||||||
|
S.group = reordercats(S.group, [{cfg.ref}; otherCats]);
|
||||||
|
|
||||||
|
% Mean-centered day covariates, computed on the post-window data.
|
||||||
|
S.day_c = S.day - mean(S.day);
|
||||||
|
S.day_c2 = S.day_c .^ 2;
|
||||||
|
|
||||||
|
% Persist the derived scenario table.
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
derivedDir = fullfile(thisDir, 'derived');
|
||||||
|
if ~exist(derivedDir, 'dir')
|
||||||
|
mkdir(derivedDir);
|
||||||
|
end
|
||||||
|
outFile = fullfile(derivedDir, [char(scenario) '.csv']);
|
||||||
|
writetable(S, outFile);
|
||||||
|
|
||||||
|
meta = struct();
|
||||||
|
meta.mergeKey = mergeKey;
|
||||||
|
meta.maxDay = maxDay;
|
||||||
|
meta.isUnmerged = strcmp(mergeKey, 'unmerged');
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
function [mergeKey, maxDay] = localParseScenario(scenario)
|
||||||
|
%LOCALPARSESCENARIO Split "<mergeKey>_full" / "<mergeKey>_d0_10" apart.
|
||||||
|
scenario = char(scenario);
|
||||||
|
if endsWith(scenario, '_full')
|
||||||
|
maxDay = 26;
|
||||||
|
mergeKey = scenario(1:end - numel('_full'));
|
||||||
|
elseif endsWith(scenario, '_d0_10')
|
||||||
|
maxDay = 10;
|
||||||
|
mergeKey = scenario(1:end - numel('_d0_10'));
|
||||||
|
else
|
||||||
|
error('tdcs_scenario_data:badScenario', ...
|
||||||
|
'Scenario "%s" does not end in "_full" or "_d0_10".', scenario);
|
||||||
|
end
|
||||||
|
|
||||||
|
if ~ismember(mergeKey, {'unmerged', 'mergeA2', 'mergeB2'})
|
||||||
|
error('tdcs_scenario_data:badScenario', ...
|
||||||
|
'Scenario "%s" has unrecognized merge key "%s".', scenario, mergeKey);
|
||||||
|
end
|
||||||
|
|
||||||
|
end
|
||||||
@@ -0,0 +1,72 @@
|
|||||||
|
classdef tLoadData < matlab.unittest.TestCase
|
||||||
|
%TLOADDATA Tests for tdcs_load_data / tdcs_config (Task 1).
|
||||||
|
%
|
||||||
|
% Authored per TDD: written before tdcs_load_data.m existed, so it is
|
||||||
|
% expected to fail until the loader is implemented. As of this
|
||||||
|
% commit MATLAB is not licensed in this environment, so this suite
|
||||||
|
% has been authored and reviewed by hand but has NOT been executed.
|
||||||
|
% Expected values (row count, group counts, load-check cells) were
|
||||||
|
% independently verified against analysis/matlab/raw/*.csv by a
|
||||||
|
% throwaway Python script (see task-1-report.md), not by running
|
||||||
|
% this test.
|
||||||
|
|
||||||
|
properties (Constant)
|
||||||
|
% Non-blank (subject, day) success cells for the 12 curated
|
||||||
|
% animals with day >= 0, counted directly from raw_success.csv.
|
||||||
|
ExpectedRowCount = 185
|
||||||
|
end
|
||||||
|
|
||||||
|
methods (Test)
|
||||||
|
function testTwelveCuratedSubjects(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
testCase.verifyEqual(numel(unique(T.subject)), 12);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testRowCountMatchesNonBlankDayGeq0Cells(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
testCase.verifyEqual(height(T), tLoadData.ExpectedRowCount);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testGroupSubjectCounts(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
subjGroup = unique(T(:, {'subject', 'group'}));
|
||||||
|
counts = countcats(subjGroup.group);
|
||||||
|
catNames = categories(subjGroup.group);
|
||||||
|
|
||||||
|
expected = containers.Map( ...
|
||||||
|
{'Naive', 'Electrode-Box-A2', 'Electrode-Box-B2', ...
|
||||||
|
'Electrode-Box-A', 'Right-Electrode'}, ...
|
||||||
|
{4, 3, 3, 1, 1});
|
||||||
|
|
||||||
|
testCase.verifyEqual(numel(catNames), double(expected.Count));
|
||||||
|
for i = 1:numel(catNames)
|
||||||
|
name = catNames{i};
|
||||||
|
testCase.verifyTrue(isKey(expected, name), ...
|
||||||
|
sprintf('Unexpected group "%s" in loaded data.', name));
|
||||||
|
testCase.verifyEqual(counts(i), expected(name), ...
|
||||||
|
sprintf('Subject count mismatch for group "%s".', name));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function testVuVuongDayZeroLoadCheck(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
row = T(T.subject == "Vu-vuong" & T.day == 0, {'success', 'total'});
|
||||||
|
testCase.verifyEqual(height(row), 1);
|
||||||
|
testCase.verifyEqual(row.success(1), 18);
|
||||||
|
testCase.verifyEqual(row.total(1), 61);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testBanhMi1DayZeroLoadCheck(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
row = T(T.subject == "Banh-mi-1" & T.day == 0, {'success', 'total'});
|
||||||
|
testCase.verifyEqual(height(row), 1);
|
||||||
|
testCase.verifyEqual(row.success(1), 13);
|
||||||
|
testCase.verifyEqual(row.total(1), 84);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testNoNegativeDayRows(testCase)
|
||||||
|
T = tdcs_load_data();
|
||||||
|
testCase.verifyTrue(all(T.day >= 0));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@@ -0,0 +1,121 @@
|
|||||||
|
classdef tModels < matlab.unittest.TestCase
|
||||||
|
%TMODELS Tests for tdcs_models (Task 3: GLMM + per-animal + classify).
|
||||||
|
%
|
||||||
|
% Authored per TDD: written before tdcs_models.m existed, so it is
|
||||||
|
% expected to fail until the model-fitting function is implemented.
|
||||||
|
% MATLAB R2025b + Statistics and Machine Learning Toolbox IS licensed
|
||||||
|
% in this environment (unlike Tasks 1-2), so this suite is actually
|
||||||
|
% run: RED before tdcs_models.m existed, GREEN after.
|
||||||
|
%
|
||||||
|
% Per-animal golden values (Box-B2 vs Box-A2, pure stats: per-subject
|
||||||
|
% mean success for "count", per-subject pooled sum(success)/sum(total)
|
||||||
|
% for "rate"; Welch two-sided TTEST2, one-sided 'right' RANKSUM) were
|
||||||
|
% independently reproduced with a throwaway Python script against
|
||||||
|
% analysis/tdcs_reach_data.csv (equal_var=False t-test,
|
||||||
|
% mannwhitneyu(..., alternative='greater')), confirming the plan's
|
||||||
|
% Task 3 golden-value table to better than 1e-3 in every cell (see
|
||||||
|
% task-3-report.md). They must match to 3 decimals (AbsTol 0.005) --
|
||||||
|
% do not loosen this tolerance.
|
||||||
|
%
|
||||||
|
% GLME-derived checks (H2 odds/incidence-rate ratios, classification)
|
||||||
|
% are checked loosely (direction/magnitude), since GLMM (subject
|
||||||
|
% random intercept) is not numerically identical to the Python's GEE
|
||||||
|
% (population-average, robust SEs) -- see the plan's Global
|
||||||
|
% Constraints and analysis/METHODS.md.
|
||||||
|
|
||||||
|
properties (TestParameter)
|
||||||
|
perAnimalCase = struct( ...
|
||||||
|
'unmerged_full', struct( ...
|
||||||
|
'scenario', 'unmerged_full', 'nB', 3, 'nA', 3, ...
|
||||||
|
'rateWelchP', 0.072, 'rateMWUp', 0.050, ...
|
||||||
|
'countWelchP', 0.055, 'countMWUp', 0.050), ...
|
||||||
|
'mergeA2_full', struct( ...
|
||||||
|
'scenario', 'mergeA2_full', 'nB', 4, 'nA', 4, ...
|
||||||
|
'rateWelchP', 0.056, 'rateMWUp', 0.029, ...
|
||||||
|
'countWelchP', 0.035, 'countMWUp', 0.029), ...
|
||||||
|
'mergeB2_full', struct( ...
|
||||||
|
'scenario', 'mergeB2_full', 'nB', 5, 'nA', 3, ...
|
||||||
|
'rateWelchP', 0.044, 'rateMWUp', 0.018, ...
|
||||||
|
'countWelchP', 0.020, 'countMWUp', 0.018), ...
|
||||||
|
'unmerged_d0_10', struct( ...
|
||||||
|
'scenario', 'unmerged_d0_10', 'nB', 3, 'nA', 3, ...
|
||||||
|
'rateWelchP', 0.070, 'rateMWUp', 0.050, ...
|
||||||
|
'countWelchP', 0.041, 'countMWUp', 0.050), ...
|
||||||
|
'mergeA2_d0_10', struct( ...
|
||||||
|
'scenario', 'mergeA2_d0_10', 'nB', 4, 'nA', 4, ...
|
||||||
|
'rateWelchP', 0.069, 'rateMWUp', 0.057, ...
|
||||||
|
'countWelchP', 0.023, 'countMWUp', 0.014), ...
|
||||||
|
'mergeB2_d0_10', struct( ...
|
||||||
|
'scenario', 'mergeB2_d0_10', 'nB', 5, 'nA', 3, ...
|
||||||
|
'rateWelchP', 0.094, 'rateMWUp', 0.071, ...
|
||||||
|
'countWelchP', 0.023, 'countMWUp', 0.018) ...
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
methods (Test)
|
||||||
|
function testPerAnimalGoldenValues(testCase, perAnimalCase)
|
||||||
|
cfg = tdcs_config();
|
||||||
|
[S, meta] = tdcs_scenario_data(perAnimalCase.scenario);
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
|
||||||
|
testCase.verifyEqual(R.perAnimal.nB, perAnimalCase.nB);
|
||||||
|
testCase.verifyEqual(R.perAnimal.nA, perAnimalCase.nA);
|
||||||
|
testCase.verifyEqual(R.perAnimal.rateWelchP, ...
|
||||||
|
perAnimalCase.rateWelchP, 'AbsTol', 0.005);
|
||||||
|
testCase.verifyEqual(R.perAnimal.rateMWUp, ...
|
||||||
|
perAnimalCase.rateMWUp, 'AbsTol', 0.005);
|
||||||
|
testCase.verifyEqual(R.perAnimal.countWelchP, ...
|
||||||
|
perAnimalCase.countWelchP, 'AbsTol', 0.005);
|
||||||
|
testCase.verifyEqual(R.perAnimal.countMWUp, ...
|
||||||
|
perAnimalCase.countMWUp, 'AbsTol', 0.005);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testGLMEObjectsReturned(testCase)
|
||||||
|
cfg = tdcs_config();
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeB2_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
|
||||||
|
testCase.verifyClass(R.countGLME, 'GeneralizedLinearMixedModel');
|
||||||
|
testCase.verifyClass(R.rateGLME, 'GeneralizedLinearMixedModel');
|
||||||
|
testCase.verifyClass(R.learnGLME, 'GeneralizedLinearMixedModel');
|
||||||
|
end
|
||||||
|
|
||||||
|
function testH2SanityMergeB2Full(testCase)
|
||||||
|
cfg = tdcs_config();
|
||||||
|
[S, meta] = tdcs_scenario_data('mergeB2_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
|
||||||
|
testCase.verifyGreaterThan(R.h2.rateOR, 1.3);
|
||||||
|
testCase.verifyLessThan(R.h2.rateOR, 1.7);
|
||||||
|
testCase.verifyLessThan(R.h2.rateOneSidedP, 0.05);
|
||||||
|
testCase.verifyGreaterThan(R.h2.countIRR, 1);
|
||||||
|
end
|
||||||
|
|
||||||
|
function testClassifyPresentForUnmergedOnly(testCase)
|
||||||
|
cfg = tdcs_config();
|
||||||
|
|
||||||
|
[S, meta] = tdcs_scenario_data('unmerged_full');
|
||||||
|
R = tdcs_models(S, meta, cfg);
|
||||||
|
testCase.verifyTrue(meta.isUnmerged);
|
||||||
|
for i = 1:numel(cfg.unknowns)
|
||||||
|
fieldName = matlab.lang.makeValidName(cfg.unknowns{i});
|
||||||
|
testCase.verifyTrue(isfield(R.classify, fieldName), ...
|
||||||
|
sprintf('R.classify missing entry for unknown "%s".', ...
|
||||||
|
cfg.unknowns{i}));
|
||||||
|
entry = R.classify.(fieldName);
|
||||||
|
testCase.verifyTrue(isfield(entry, 'count'));
|
||||||
|
testCase.verifyTrue(isfield(entry, 'rate'));
|
||||||
|
testCase.verifyTrue(isfield(entry.count, 'label'));
|
||||||
|
testCase.verifyTrue(isfield(entry.count.vsA2, 'ratio'));
|
||||||
|
testCase.verifyTrue(isfield(entry.count.vsA2, 'p'));
|
||||||
|
testCase.verifyTrue(isfield(entry.count.vsB2, 'ratio'));
|
||||||
|
testCase.verifyTrue(isfield(entry.count.vsB2, 'p'));
|
||||||
|
end
|
||||||
|
|
||||||
|
[Sm, metam] = tdcs_scenario_data('mergeB2_full');
|
||||||
|
Rm = tdcs_models(Sm, metam, cfg);
|
||||||
|
testCase.verifyFalse(metam.isUnmerged);
|
||||||
|
testCase.verifyEmpty(fieldnames(Rm.classify));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@@ -0,0 +1,51 @@
|
|||||||
|
classdef tReport < matlab.unittest.TestCase
|
||||||
|
%TREPORT Tests for tdcs_glm / tdcs_report (Task 4: entry point + report).
|
||||||
|
%
|
||||||
|
% MATLAB R2025b + Statistics and Machine Learning Toolbox IS licensed
|
||||||
|
% in this environment, so this suite is actually run. It does not
|
||||||
|
% re-check golden values (tModels already does that); it checks the
|
||||||
|
% wiring: tdcs_glm(scenario) runs end-to-end without error and writes
|
||||||
|
% results/<scenario>.txt, and the otherwise-branch errors on an
|
||||||
|
% unrecognized scenario name.
|
||||||
|
|
||||||
|
methods (Test)
|
||||||
|
function testRunsScenarioAndWritesResultFile(testCase)
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
resultsDir = fullfile(fileparts(thisDir), 'results');
|
||||||
|
outFile = fullfile(resultsDir, 'mergeB2_full.txt');
|
||||||
|
if exist(outFile, 'file')
|
||||||
|
delete(outFile);
|
||||||
|
end
|
||||||
|
|
||||||
|
tdcs_glm('mergeB2_full');
|
||||||
|
|
||||||
|
testCase.verifyTrue(exist(outFile, 'file') == 2, ...
|
||||||
|
'tdcs_glm did not write results/mergeB2_full.txt.');
|
||||||
|
txt = fileread(outFile);
|
||||||
|
testCase.verifyNotEmpty(txt);
|
||||||
|
testCase.verifyTrue(contains(txt, 'DESCRIPTIVES'));
|
||||||
|
testCase.verifyTrue(contains(txt, 'INTERPRETATION'));
|
||||||
|
testCase.verifyTrue(contains(txt, 'CAVEATS'));
|
||||||
|
% Command-Window-only markup must be stripped from the report.
|
||||||
|
testCase.verifyFalse(contains(txt, '<strong>'));
|
||||||
|
end
|
||||||
|
|
||||||
|
function testUnmergedScenarioIncludesClassification(testCase)
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
resultsDir = fullfile(fileparts(thisDir), 'results');
|
||||||
|
outFile = fullfile(resultsDir, 'unmerged_full.txt');
|
||||||
|
|
||||||
|
tdcs_glm('unmerged_full');
|
||||||
|
|
||||||
|
txt = fileread(outFile);
|
||||||
|
testCase.verifyTrue(contains(txt, 'Classification'));
|
||||||
|
testCase.verifyTrue(contains(txt, 'Electrode-Box-A'));
|
||||||
|
testCase.verifyTrue(contains(txt, 'Right-Electrode'));
|
||||||
|
end
|
||||||
|
|
||||||
|
function testBadScenarioErrors(testCase)
|
||||||
|
testCase.verifyError(@() tdcs_glm('not_a_real_scenario'), ...
|
||||||
|
'tdcs_glm:badScenario');
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@@ -0,0 +1,135 @@
|
|||||||
|
classdef tScenarioData < matlab.unittest.TestCase
|
||||||
|
%TSCENARIODATA Tests for tdcs_scenario_data (Task 2).
|
||||||
|
%
|
||||||
|
% Authored per TDD: written before tdcs_scenario_data.m existed, so
|
||||||
|
% it is expected to fail until the function is implemented. As of
|
||||||
|
% this commit MATLAB is not licensed in this environment, so this
|
||||||
|
% suite has been authored and reviewed by hand but has NOT been
|
||||||
|
% executed. Expected per-scenario subject-per-group counts were
|
||||||
|
% hand-derived from the Task 1 curated counts (Naive=4,
|
||||||
|
% Electrode-Box-A2=3, Electrode-Box-B2=3, Electrode-Box-A=1,
|
||||||
|
% Right-Electrode=1) plus each scenario's merge map, and cross-
|
||||||
|
% checked against task-2-brief.md's table (see task-2-report.md).
|
||||||
|
|
||||||
|
properties (TestParameter)
|
||||||
|
scenarioCase = struct( ...
|
||||||
|
'unmerged_full', struct( ...
|
||||||
|
'scenario', 'unmerged_full', ...
|
||||||
|
'maxDay', 26, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A', 'Electrode-Box-A2', ...
|
||||||
|
'Electrode-Box-B2', 'Naive', 'Right-Electrode'}, ...
|
||||||
|
{1, 3, 3, 4, 1})), ...
|
||||||
|
'unmerged_d0_10', struct( ...
|
||||||
|
'scenario', 'unmerged_d0_10', ...
|
||||||
|
'maxDay', 10, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A', 'Electrode-Box-A2', ...
|
||||||
|
'Electrode-Box-B2', 'Naive', 'Right-Electrode'}, ...
|
||||||
|
{1, 3, 3, 4, 1})), ...
|
||||||
|
'mergeA2_full', struct( ...
|
||||||
|
'scenario', 'mergeA2_full', ...
|
||||||
|
'maxDay', 26, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A2', 'Electrode-Box-B2', 'Naive'}, ...
|
||||||
|
{4, 4, 4})), ...
|
||||||
|
'mergeA2_d0_10', struct( ...
|
||||||
|
'scenario', 'mergeA2_d0_10', ...
|
||||||
|
'maxDay', 10, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A2', 'Electrode-Box-B2', 'Naive'}, ...
|
||||||
|
{4, 4, 4})), ...
|
||||||
|
'mergeB2_full', struct( ...
|
||||||
|
'scenario', 'mergeB2_full', ...
|
||||||
|
'maxDay', 26, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A2', 'Electrode-Box-B2', 'Naive'}, ...
|
||||||
|
{3, 5, 4})), ...
|
||||||
|
'mergeB2_d0_10', struct( ...
|
||||||
|
'scenario', 'mergeB2_d0_10', ...
|
||||||
|
'maxDay', 10, ...
|
||||||
|
'counts', containers.Map( ...
|
||||||
|
{'Electrode-Box-A2', 'Electrode-Box-B2', 'Naive'}, ...
|
||||||
|
{3, 5, 4})) ...
|
||||||
|
)
|
||||||
|
end
|
||||||
|
|
||||||
|
methods (Test)
|
||||||
|
function testGroupCountsWindowAndReference(testCase, scenarioCase)
|
||||||
|
[S, meta] = tdcs_scenario_data(scenarioCase.scenario);
|
||||||
|
|
||||||
|
% Day window.
|
||||||
|
testCase.verifyEqual(meta.maxDay, scenarioCase.maxDay);
|
||||||
|
testCase.verifyLessThanOrEqual(max(S.day), scenarioCase.maxDay);
|
||||||
|
|
||||||
|
% meta fields derived from the scenario name.
|
||||||
|
expectedMergeKey = regexprep(scenarioCase.scenario, '_full$|_d0_10$', '');
|
||||||
|
testCase.verifyEqual(meta.mergeKey, expectedMergeKey);
|
||||||
|
testCase.verifyEqual(meta.isUnmerged, strcmp(expectedMergeKey, 'unmerged'));
|
||||||
|
|
||||||
|
% Reference (first) category of S.group is Electrode-Box-B2.
|
||||||
|
allCats = categories(S.group);
|
||||||
|
testCase.verifyEqual(allCats{1}, 'Electrode-Box-B2');
|
||||||
|
|
||||||
|
% Per-group subject counts.
|
||||||
|
actualCounts = tScenarioData.localSubjectCountsByGroup(S);
|
||||||
|
expectedCounts = scenarioCase.counts;
|
||||||
|
|
||||||
|
testCase.verifyEqual(actualCounts.Count, expectedCounts.Count);
|
||||||
|
actualNames = keys(actualCounts);
|
||||||
|
for i = 1:numel(actualNames)
|
||||||
|
name = actualNames{i};
|
||||||
|
testCase.verifyTrue(isKey(expectedCounts, name), ...
|
||||||
|
sprintf('Unexpected group "%s" in scenario "%s".', ...
|
||||||
|
name, scenarioCase.scenario));
|
||||||
|
testCase.verifyEqual(actualCounts(name), expectedCounts(name), ...
|
||||||
|
sprintf('Subject count mismatch for group "%s" in scenario "%s".', ...
|
||||||
|
name, scenarioCase.scenario));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
function testMergeB2FullCsvIsWrittenAndReloads(testCase)
|
||||||
|
[S, ~] = tdcs_scenario_data('mergeB2_full');
|
||||||
|
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
matlabDir = fileparts(thisDir);
|
||||||
|
csvFile = fullfile(matlabDir, 'derived', 'mergeB2_full.csv');
|
||||||
|
|
||||||
|
testCase.verifyTrue(isfile(csvFile), ...
|
||||||
|
'derived/mergeB2_full.csv was not written.');
|
||||||
|
|
||||||
|
Sreload = readtable(csvFile);
|
||||||
|
Sreload.group = categorical(Sreload.group);
|
||||||
|
|
||||||
|
countsOrig = tScenarioData.localSubjectCountsByGroup(S);
|
||||||
|
countsReload = tScenarioData.localSubjectCountsByGroup(Sreload);
|
||||||
|
|
||||||
|
testCase.verifyEqual(countsReload.Count, countsOrig.Count);
|
||||||
|
namesOrig = keys(countsOrig);
|
||||||
|
for i = 1:numel(namesOrig)
|
||||||
|
name = namesOrig{i};
|
||||||
|
testCase.verifyTrue(isKey(countsReload, name), ...
|
||||||
|
sprintf('Group "%s" missing after CSV reload.', name));
|
||||||
|
testCase.verifyEqual(countsReload(name), countsOrig(name), ...
|
||||||
|
sprintf('Subject count mismatch for group "%s" after CSV reload.', name));
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
methods (Static)
|
||||||
|
function m = localSubjectCountsByGroup(T)
|
||||||
|
%LOCALSUBJECTCOUNTSBYGROUP containers.Map group name -> #unique subjects.
|
||||||
|
subjGroup = unique(T(:, {'subject', 'group'}));
|
||||||
|
group = subjGroup.group;
|
||||||
|
if ~iscategorical(group)
|
||||||
|
group = categorical(group);
|
||||||
|
end
|
||||||
|
counts = countcats(group);
|
||||||
|
catNames = categories(group);
|
||||||
|
m = containers.Map('KeyType', 'char', 'ValueType', 'double');
|
||||||
|
for i = 1:numel(catNames)
|
||||||
|
m(catNames{i}) = counts(i);
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
@@ -0,0 +1,151 @@
|
|||||||
|
# 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.
|
||||||
@@ -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 0–26 |
|
||||||
|
| `unmerged_d0_10` | 5 groups | days 0–10 |
|
||||||
|
| `mergeA2_full` | Box-A→A2, Right→B2 | 0–26 |
|
||||||
|
| `mergeA2_d0_10` | Box-A→A2, Right→B2 | 0–10 |
|
||||||
|
| `mergeB2_full` | Box-A→B2, Right→B2 | 0–26 |
|
||||||
|
| `mergeB2_d0_10` | Box-A→B2, Right→B2 | 0–10 |
|
||||||
|
|
||||||
|
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=3–4/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.4–1.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).
|
||||||
Reference in New Issue
Block a user