Experiments DB Dev
173 lines
7.5 KiB
Matlab
173 lines
7.5 KiB
Matlab
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
|