matlab(tdcs): report, switch/case entry, runners, README
This commit is contained in:
Experiments DB Dev
@@ -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/<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.
|
||||
@@ -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,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,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,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
|
||||
Reference in New Issue
Block a user