From 08c14c476c3650c307c5d25a172876fefc216d7f Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Mon, 20 Jul 2026 09:43:08 -0400 Subject: [PATCH] feat(matlab): add linear days x tDCS mixed-model scenarios (lme_*) Adds 6 switch cases replicating the paper's linear mixed model success ~ day * tDCS + (1|subject) on the Box-B2(tDCS=1) vs Box-A2(tDCS=0) arm, reporting the days x tDCS interaction, days, and tDCS effects as t/F/p. mergeA2_full reproduces the paper's interaction (F=7.09 vs 7.12, p=0.009 vs 0.008); interpretation is sign-aware and honestly notes our data diverges from the paper (Box-B2 ahead on day 1, gap narrows -- not 'equal on day 1, accumulates'). Adds tLme tests. Suite 31/31. Co-Authored-By: Claude Opus 4.8 --- analysis/matlab/run_all.m | 8 +++- analysis/matlab/tdcs_glm.m | 32 ++++++++++++++- analysis/matlab/tdcs_lme.m | 50 +++++++++++++++++++++++ analysis/matlab/tdcs_lme_report.m | 67 +++++++++++++++++++++++++++++++ analysis/matlab/tests/tLme.m | 63 +++++++++++++++++++++++++++++ 5 files changed, 217 insertions(+), 3 deletions(-) create mode 100644 analysis/matlab/tdcs_lme.m create mode 100644 analysis/matlab/tdcs_lme_report.m create mode 100644 analysis/matlab/tests/tLme.m diff --git a/analysis/matlab/run_all.m b/analysis/matlab/run_all.m index 9a4ba47..9cd6d11 100644 --- a/analysis/matlab/run_all.m +++ b/analysis/matlab/run_all.m @@ -1,9 +1,13 @@ -%RUN_ALL Run all 6 tDCS GLM scenarios and write results/.txt for each. +%RUN_ALL Run all tDCS scenarios and write results/.txt for each. % Usage: matlab -batch "run_all" +% Covers the 6 GLMM scenarios and their 6 linear "days x tDCS" (lme_*) variants. scenarios = {'unmerged_full', 'unmerged_d0_10', ... 'mergeA2_full', 'mergeA2_d0_10', ... - 'mergeB2_full', 'mergeB2_d0_10'}; + 'mergeB2_full', 'mergeB2_d0_10', ... + 'lme_unmerged_full', 'lme_unmerged_d0_10', ... + 'lme_mergeA2_full', 'lme_mergeA2_d0_10', ... + 'lme_mergeB2_full', 'lme_mergeB2_d0_10'}; for i = 1:numel(scenarios) fprintf('\n\n### Running scenario: %s ###\n', scenarios{i}); diff --git a/analysis/matlab/tdcs_glm.m b/analysis/matlab/tdcs_glm.m index 2dbcd63..762366d 100644 --- a/analysis/matlab/tdcs_glm.m +++ b/analysis/matlab/tdcs_glm.m @@ -9,6 +9,10 @@ function tdcs_glm(scenario) % unmerged_full, unmerged_d0_10, % mergeA2_full, mergeA2_d0_10, % mergeB2_full, mergeB2_d0_10 +% or a linear "days x tDCS" mixed model (Box-B2 vs Box-A2, per the paper): +% lme_unmerged_full, lme_unmerged_d0_10, +% lme_mergeA2_full, lme_mergeA2_d0_10, +% lme_mergeB2_full, lme_mergeB2_d0_10 % % Each case below is a self-contained, independently invocable call of % the same three-function pipeline; they are enumerated explicitly @@ -54,11 +58,37 @@ switch scenario R = tdcs_models(S, meta, cfg); tdcs_report(S, R, meta, cfg, 'mergeB2_d0_10'); + % --- Linear "days x tDCS" mixed model (Box-B2 vs Box-A2), per the paper --- + case 'lme_unmerged_full' + L = tdcs_lme(tdcs_scenario_data('unmerged_full'), cfg); + tdcs_lme_report(L, 'lme_unmerged_full', cfg); + + case 'lme_unmerged_d0_10' + L = tdcs_lme(tdcs_scenario_data('unmerged_d0_10'), cfg); + tdcs_lme_report(L, 'lme_unmerged_d0_10', cfg); + + case 'lme_mergeA2_full' + L = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg); + tdcs_lme_report(L, 'lme_mergeA2_full', cfg); + + case 'lme_mergeA2_d0_10' + L = tdcs_lme(tdcs_scenario_data('mergeA2_d0_10'), cfg); + tdcs_lme_report(L, 'lme_mergeA2_d0_10', cfg); + + case 'lme_mergeB2_full' + L = tdcs_lme(tdcs_scenario_data('mergeB2_full'), cfg); + tdcs_lme_report(L, 'lme_mergeB2_full', cfg); + + case 'lme_mergeB2_d0_10' + L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_10'), cfg); + tdcs_lme_report(L, 'lme_mergeB2_d0_10', cfg); + 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); + 'mergeB2_full, mergeB2_d0_10, and their lme_* variants ' ... + '(lme_unmerged_full, ..., lme_mergeB2_d0_10).'], scenario); end end diff --git a/analysis/matlab/tdcs_lme.m b/analysis/matlab/tdcs_lme.m new file mode 100644 index 0000000..1ba3b6d --- /dev/null +++ b/analysis/matlab/tdcs_lme.m @@ -0,0 +1,50 @@ +function L = tdcs_lme(S, cfg) +%TDCS_LME Linear mixed-effects "days x tDCS" model on the Box-B2 vs Box-A2 arm. +% L = TDCS_LME(S, CFG) subsets one scenario table S to the two anchor groups, +% codes a binary tDCS factor (Box-B2 = 1, Box-A2 = 0), and fits the LINEAR +% mixed model +% +% success ~ day * tDCS + (1|subject) +% +% replicating the published formulation (a linear mixed-effects model for the +% number of successful reaches with a days x tDCS interaction, a main effect +% of days, and a main effect of tDCS). `day` is used raw (not centered), so +% the tDCS main effect is the group difference on the first day ("Day 1"). +% +% L fields: +% .lme the LinearMixedModel object +% .nSubjects number of subjects (Box-A2 + Box-B2) +% .nObs number of sessions +% .interaction / .day / .tDCS effect structs, each with .estimate .se +% .t .df .p (from Coefficients) and .F .df1 .df2 .Fp (from +% ANOVA). For these single-df terms F == t^2 and Fp == p. +% +% NOTE: this is a LINEAR mixed model (Gaussian on the raw count), matching +% the paper's method, unlike the Poisson/Binomial GLMMs in tdcs_models. + +T = S(ismember(S.group, {cfg.anchorLow, cfg.anchorHigh}), :); +T.group = removecats(T.group); +T.tDCS = double(T.group == cfg.anchorHigh); % Box-B2 = 1, Box-A2 = 0 + +lme = fitlme(T, 'success ~ day*tDCS + (1|subject)'); +C = lme.Coefficients; +A = anova(lme); + +L.lme = lme; +L.nSubjects = numel(unique(T.subject)); +L.nObs = height(T); +L.interaction = localTerm(C, A, 'day:tDCS'); +L.day = localTerm(C, A, 'day'); +L.tDCS = localTerm(C, A, 'tDCS'); + +end + +function e = localTerm(C, A, name) +%LOCALTERM Pull one term's coefficient (t/df/p) and ANOVA (F/df/p) stats. +ci = strcmp(C.Name, name); +ai = strcmp(A.Term, name); +e = struct( ... + 'estimate', C.Estimate(ci), 'se', C.SE(ci), ... + 't', C.tStat(ci), 'df', C.DF(ci), 'p', C.pValue(ci), ... + 'F', A.FStat(ai), 'df1', A.DF1(ai), 'df2', A.DF2(ai), 'Fp', A.pValue(ai)); +end diff --git a/analysis/matlab/tdcs_lme_report.m b/analysis/matlab/tdcs_lme_report.m new file mode 100644 index 0000000..e399025 --- /dev/null +++ b/analysis/matlab/tdcs_lme_report.m @@ -0,0 +1,67 @@ +function tdcs_lme_report(L, scenario, cfg) +%TDCS_LME_REPORT Print + save the days x tDCS linear-mixed-model report. +% TDCS_LME_REPORT(L, SCENARIO, CFG) formats the result of TDCS_LME to the +% console and to analysis/matlab/results/.txt, in the style of the +% published result (interaction, days, tDCS effects) plus interpretation and +% the paper's reference numbers. + +bar = repmat('=', 1, 78); +s = sprintf('%s\n', bar); +s = [s sprintf('LINEAR MIXED MODEL (days x tDCS) -- scenario: %s\n', scenario)]; +s = [s sprintf('%s\n', bar)]; +s = [s sprintf('model: success ~ day * tDCS + (1|subject) [tDCS: %s = 1 vs %s = 0]\n', ... + cfg.anchorHigh, cfg.anchorLow)]; +s = [s sprintf('N = %d subjects, %d sessions (day used raw: tDCS main = first-day difference)\n\n', ... + L.nSubjects, L.nObs)]; + +s = [s sprintf('%-27s %-14s %-13s %s\n', 'effect', 't (df)', 'F (df1)', 'p')]; +s = [s sprintf('%s\n', repmat('-', 1, 70))]; +s = [s localRow('days x tDCS (interaction)', L.interaction)]; +s = [s localRow('days (learning)', L.day)]; +s = [s localRow('tDCS (main, at day 0)', L.tDCS)]; + +s = [s sprintf('\nINTERPRETATION\n')]; +if L.interaction.p >= 0.05 + interTxt = 'slopes are parallel (no differential change over training)'; +elseif L.interaction.estimate > 0 + interTxt = 'the tDCS (Box-B2) group improves FASTER -- benefit ACCUMULATES over training'; +else + interTxt = 'the tDCS (Box-B2) group improves SLOWER -- groups CONVERGE (Box-B2 is ahead early, the gap narrows)'; +end +s = [s sprintf(' - days x tDCS interaction: %s (p=%.4f, slope diff=%.2f) -> %s.\n', ... + localSig(L.interaction.p), L.interaction.p, L.interaction.estimate, interTxt)]; +s = [s sprintf(' - days (learning): %s (p=%.2g) -> performance improves with training.\n', ... + localSig(L.day.p), L.day.p)]; +s = [s sprintf(' - tDCS main effect at day 0: %s (p=%.4f) -> the groups %s on the first day.\n', ... + localSig(L.tDCS.p), L.tDCS.p, ... + localPick(L.tDCS.p < 0.05, 'already DIFFER', 'are comparable'))]; + +s = [s sprintf('\nPaper reference (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008;\n')]; +s = [s sprintf(' days t(227)=9.64, F(1)=267.64, p=1.2e-18; tDCS t(227)=0.23, F(1)=0.053, p=0.81.\n')]; +s = [s sprintf(' (Our N and exact statistics differ; this replicates the MODEL FORM on our data.)\n')]; + +fprintf('%s', s); + +thisDir = fileparts(mfilename('fullpath')); +resDir = fullfile(thisDir, 'results'); +if ~exist(resDir, 'dir'); mkdir(resDir); end +fid = fopen(fullfile(resDir, [scenario '.txt']), 'w'); +if fid < 0 + error('tdcs_lme_report:fopen', 'Cannot open results file for "%s".', scenario); +end +cleanup = onCleanup(@() fclose(fid)); %#ok +fprintf(fid, '%s', s); + +end + +function r = localRow(name, e) +r = sprintf('%-27s t(%d)=%6.2f F(%d)=%8.3f p=%.4g\n', name, e.df, e.t, e.df1, e.F, e.p); +end + +function t = localSig(p) +if p < 0.05; t = 'SIGNIFICANT'; else; t = 'n.s.'; end +end + +function t = localPick(b, yes, no) +if b; t = yes; else; t = no; end +end diff --git a/analysis/matlab/tests/tLme.m b/analysis/matlab/tests/tLme.m new file mode 100644 index 0000000..686fc1f --- /dev/null +++ b/analysis/matlab/tests/tLme.m @@ -0,0 +1,63 @@ +classdef tLme < matlab.unittest.TestCase + %TLME Tests for the linear "days x tDCS" mixed model (tdcs_lme) and its + % switch/case entry. Golden values computed from this dataset; the + % mergeA2_full interaction closely reproduces the paper's F=7.12/p=0.008. + + methods (Test) + + function testSubjectCountsAndDayEffect(testCase) + % N (Box-A2 + Box-B2 subjects) per scenario, and the always-present + % strong learning (day) effect. + cfg = tdcs_config(); + expectedN = struct('unmerged_full', 6, 'mergeA2_full', 8, 'mergeB2_full', 8); + for sc = fieldnames(expectedN)' + L = tdcs_lme(tdcs_scenario_data(sc{1}), cfg); + testCase.verifyEqual(L.nSubjects, expectedN.(sc{1}), ... + sprintf('%s: wrong subject count', sc{1})); + testCase.verifyEqual(L.day.df1, 1); % F(1) term + testCase.verifyLessThan(L.day.p, 1e-6, ... + sprintf('%s: day (learning) effect should be highly significant', sc{1})); + end + end + + function testMergeA2ReplicatesPaperInteraction(testCase) + % mergeA2_full reproduces the paper's days x tDCS interaction + % (paper: F(1)=7.12, p=0.008; |t|=2.68). + cfg = tdcs_config(); + L = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg); + testCase.verifyEqual(L.interaction.df1, 1); + testCase.verifyEqual(L.interaction.F, 7.09, 'AbsTol', 0.1); + testCase.verifyEqual(L.interaction.p, 0.009, 'AbsTol', 0.002); + testCase.verifyEqual(abs(L.interaction.t), 2.66, 'AbsTol', 0.05); + % single-df term: F == t^2 + testCase.verifyEqual(L.interaction.F, L.interaction.t^2, 'RelTol', 1e-6); + end + + function testEffectStructShape(testCase) + cfg = tdcs_config(); + L = tdcs_lme(tdcs_scenario_data('unmerged_full'), cfg); + for f = {'interaction', 'day', 'tDCS'} + e = L.(f{1}); + for g = {'estimate', 'se', 't', 'df', 'p', 'F', 'df1', 'df2', 'Fp'} + testCase.verifyTrue(isfield(e, g{1}), ... + sprintf('%s missing field %s', f{1}, g{1})); + end + end + end + + function testSwitchRunsAndWritesResults(testCase) + % The lme_* switch cases run end-to-end and write results/.txt. + here = fileparts(fileparts(mfilename('fullpath'))); % analysis/matlab + outFile = fullfile(here, 'results', 'lme_mergeA2_full.txt'); + if exist(outFile, 'file'); delete(outFile); end + evalc("tdcs_glm('lme_mergeA2_full')"); + testCase.verifyTrue(exist(outFile, 'file') == 2, ... + 'lme_mergeA2_full.txt was not written'); + end + + function testBadScenarioErrors(testCase) + testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario'); + end + + end +end