diff --git a/analysis/matlab/README.md b/analysis/matlab/README.md new file mode 100644 index 0000000..8fbd77c --- /dev/null +++ b/analysis/matlab/README.md @@ -0,0 +1,99 @@ +# 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/.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/.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. diff --git a/analysis/matlab/run_all.m b/analysis/matlab/run_all.m new file mode 100644 index 0000000..9a4ba47 --- /dev/null +++ b/analysis/matlab/run_all.m @@ -0,0 +1,13 @@ +%RUN_ALL Run all 6 tDCS GLM scenarios and write results/.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)); diff --git a/analysis/matlab/run_tests.m b/analysis/matlab/run_tests.m new file mode 100644 index 0000000..7626020 --- /dev/null +++ b/analysis/matlab/run_tests.m @@ -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 diff --git a/analysis/matlab/tdcs_glm.m b/analysis/matlab/tdcs_glm.m new file mode 100644 index 0000000..2dbcd63 --- /dev/null +++ b/analysis/matlab/tdcs_glm.m @@ -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/.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 diff --git a/analysis/matlab/tdcs_report.m b/analysis/matlab/tdcs_report.m new file mode 100644 index 0000000..124523b --- /dev/null +++ b/analysis/matlab/tdcs_report.m @@ -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/.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 +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 +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 +% ... 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, '', ''); +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 + s = [s localClassifyLine('count/level', entry.count)]; %#ok + s = [s localClassifyLine('rate/level ', entry.rate)]; %#ok + 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 diff --git a/analysis/matlab/tests/tReport.m b/analysis/matlab/tests/tReport.m new file mode 100644 index 0000000..51ca31e --- /dev/null +++ b/analysis/matlab/tests/tReport.m @@ -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/.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, '')); + 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