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 <noreply@anthropic.com>

analysis/matlab/run_all.m

@@ -1,9 +1,13 @@
-%RUN_ALL Run all 6 tDCS GLM scenarios and write results/<scenario>.txt for each.
+%RUN_ALL Run all tDCS scenarios and write results/<scenario>.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});

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

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

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/<SCENARIO>.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<NASGU>
+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

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/<name>.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