feat(matlab): verbatim paper LME (behavior ~ stim+day+stim:day+(1|rat)) + power sim on it

Adds tdcs_paper_lme replicating the paper's exact formula/variable names (paper_*
switch cases): mergeA2 reproduces the interaction F(1)=7.09 vs 7.12, p=0.009 vs
0.008, with a coverage-confound warning. Rewrites tdcs_power_sim to fit that same
formula (interaction term canonicalized to day:stim). Adds a replication test.
Suite 37/37.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

analysis/matlab/README.md

@@ -142,3 +142,16 @@ term, and MATLAB's native `fitglme` already *is* the random-intercept model.
   ```
   Produces `curated_all.csv`, `scenarios/<merge>_<window>.csv` (9), `anchors/<merge>_B2vsA2.csv`
   (3, with the binary tDCS factor), and `phases/<merge>_phase_<lo>-<hi>.csv` (9), plus a MANIFEST.
+
+## Paper replication (verbatim formula)
+
+The paper's model, `fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)')`, is
+replicated verbatim (behavior = successful reaches, stim = tDCS [Box-B2=1 vs Box-A2=0],
+rat = subject) over the full training range:
+```
+matlab -batch "tdcs_glm('paper_mergeA2')"   # or paper_unmerged / paper_mergeB2
+```
+`paper_mergeA2` reproduces the paper's interaction (F(1)=7.09 vs 7.12, p=0.009 vs 0.008); it
+also warns that the full-range interaction is confounded by unequal day coverage (Box-A2 ends
+~day 13) — see the _d0_13 and phased analyses. The Monte-Carlo power analysis (`run_power`,
+`tdcs_power_sim`) uses this same `behavior ~ stim + day + stim:day + (1|rat)` model.

analysis/matlab/run_all.m

@@ -9,7 +9,8 @@ scenarios = {'unmerged_full', 'unmerged_d0_10', ...
              'lme_mergeA2_full', 'lme_mergeA2_d0_10', ...
              'lme_mergeB2_full', 'lme_mergeB2_d0_10', ...
              'lme_unmerged_d0_13', 'lme_mergeA2_d0_13', 'lme_mergeB2_d0_13', ...
-             'phase_unmerged', 'phase_mergeA2', 'phase_mergeB2'};
+             'phase_unmerged', 'phase_mergeA2', 'phase_mergeB2', ...
+             'paper_unmerged', 'paper_mergeA2', 'paper_mergeB2'};
 
 for i = 1:numel(scenarios)
     fprintf('\n\n### Running scenario: %s ###\n', scenarios{i});

analysis/matlab/tdcs_glm.m

@@ -107,6 +107,16 @@ switch scenario
     case 'phase_mergeB2'
         tdcs_phase_lme('mergeB2', cfg);
 
+    % --- Paper replication: behavior ~ stim + day + stim:day + (1|rat) ---
+    case 'paper_unmerged'
+        tdcs_paper_lme('unmerged', cfg);
+
+    case 'paper_mergeA2'
+        tdcs_paper_lme('mergeA2', cfg);
+
+    case 'paper_mergeB2'
+        tdcs_paper_lme('mergeB2', cfg);
+
     otherwise
         error('tdcs_glm:badScenario', ...
             ['Unrecognized scenario "%s". Valid scenarios are: ' ...

analysis/matlab/tdcs_paper_lme.m

@@ -0,0 +1,91 @@
+function L = tdcs_paper_lme(mergeKey, cfg)
+%TDCS_PAPER_LME Replicate the paper's linear mixed model, verbatim formula:
+%       behavior ~ stim + day + stim:day + (1|rat)
+%   fit with fitlme on the Box-B2 (stim = 1, tDCS) vs Box-A2 (stim = 0, control)
+%   subset of the MERGEKEY grouping ('unmerged'|'mergeA2'|'mergeB2'), over the
+%   full training range. `behavior` = successful reaches, `day` = training day
+%   (raw; day 0 = the paper's "Day 1"), `rat` = subject. Reports the three
+%   effects (stim x day interaction, day, stim) as t(df) / F(1) / p with the
+%   interaction 95% CI, alongside the paper's reference values, and writes
+%   results/paper_<MERGEKEY>.txt.
+%
+%   (This is the same model as tdcs_lme -- success ~ day*tDCS + (1|subject) --
+%   written with the paper's exact term order and variable names.)
+%
+%   L fields: .model .nRats .nObs and .stim/.day/.interaction effect structs
+%   (.estimate .se .t .df .p from Coefficients; .F .df1 .df2 .Fp from ANOVA;
+%   .ci = coefficient 95% CI).
+
+Sfull = tdcs_scenario_data([mergeKey '_full']);
+A = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
+
+tbl = table();
+tbl.behavior = A.success;
+tbl.day = A.day;
+tbl.stim = double(A.group == cfg.anchorHigh);   % Box-B2 = 1, Box-A2 = 0
+tbl.rat = A.subject;
+
+model = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
+C = model.Coefficients; An = anova(model); ci = coefCI(model);
+
+L.model = model;
+L.nRats = numel(unique(tbl.rat));
+L.nObs = height(tbl);
+L.maxDayCtrl = max(tbl.day(tbl.stim == 0));
+L.maxDayStim = max(tbl.day(tbl.stim == 1));
+L.stim = localTerm(C, An, ci, 'stim');
+L.day = localTerm(C, An, ci, 'day');
+L.interaction = localTerm(C, An, ci, 'day:stim');   % MATLAB canonicalizes stim:day -> day:stim
+
+localReport(L, mergeKey, cfg);
+
+end
+
+function e = localTerm(C, An, ci, name)
+i = strcmp(C.Name, name);
+if ~any(i)
+    error('tdcs_paper_lme:missingTerm', 'No "%s" coefficient (have: %s).', ...
+        name, strjoin(C.Name, ', '));
+end
+ai = strcmp(An.Term, name);
+e = struct('estimate', C.Estimate(i), 'se', C.SE(i), 't', C.tStat(i), ...
+    'df', C.DF(i), 'p', C.pValue(i), 'F', An.FStat(ai), 'df1', An.DF1(ai), ...
+    'df2', An.DF2(ai), 'Fp', An.pValue(ai), 'ci', ci(i, :));
+end
+
+function localReport(L, mergeKey, cfg)
+bar = repmat('=', 1, 78);
+s = sprintf('%s\n', bar);
+s = [s sprintf('PAPER LME REPLICATION -- %s\n', mergeKey)];
+s = [s sprintf('%s\n', bar)];
+s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat)   [stim: %s=1 vs %s=0]\n', ...
+    cfg.anchorHigh, cfg.anchorLow)];
+s = [s sprintf('N = %d rats, %d sessions   (day raw; day 0 = paper "Day 1")\n', L.nRats, L.nObs)];
+s = [s sprintf('day coverage: stim(B2) 0..%d, control(A2) 0..%d\n', L.maxDayStim, L.maxDayCtrl)];
+if abs(L.maxDayStim - L.maxDayCtrl) > 2
+    s = [s sprintf(['** WARNING: unequal day coverage -- the full-range stim:day interaction\n' ...
+        '   extrapolates the control group''s line and is CONFOUNDED here (the paper''s\n' ...
+        '   groups had equal coverage). See the _d0_13 fair-window and phased analyses. **\n'])];
+end
+s = [s sprintf('\n%-26s %-13s %-13s %s\n', 'effect', 't (df)', 'F (df1)', 'p')];
+s = [s sprintf('%s\n', repmat('-', 1, 66))];
+s = [s localRow('stim x day (interaction)', L.interaction)];
+s = [s localRow('day (learning)', L.day)];
+s = [s localRow('stim (main, Day 1)', L.stim)];
+s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', L.interaction.ci(1), L.interaction.ci(2))];
+s = [s sprintf(['\nPaper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008; ' ...
+    'day t(227)=9.64,\n  F(1)=267.64, p=1.2e-18; stim t(227)=0.23, F(1)=0.053, p=0.81.\n'])];
+
+fprintf('%s', s);
+thisDir = fileparts(mfilename('fullpath'));
+resDir = fullfile(thisDir, 'results');
+if ~exist(resDir, 'dir'); mkdir(resDir); end
+fid = fopen(fullfile(resDir, ['paper_' mergeKey '.txt']), 'w');
+if fid < 0; error('tdcs_paper_lme:fopen', 'Cannot open results file.'); end
+cleanup = onCleanup(@() fclose(fid)); %#ok<NASGU>
+fprintf(fid, '%s', s);
+end
+
+function r = localRow(name, e)
+r = sprintf('%-26s t(%d)=%6.2f  F(%d)=%8.3f  p=%.4g\n', name, e.df, e.t, e.df1, e.F, e.p);
+end

analysis/matlab/tdcs_power_sim.m

@@ -1,60 +1,64 @@
-function PW = tdcs_power_sim(mergeKey, phase, Ns, effMuls, nrep, cfg)
+function PW = tdcs_power_sim(mergeKey, window, Ns, effMuls, nrep, cfg)
-%TDCS_POWER_SIM Monte-Carlo power for the phased days x tDCS interaction.
+%TDCS_POWER_SIM Monte-Carlo power for the paper's stim x day interaction.
-%   PW = TDCS_POWER_SIM(MERGEKEY, PHASE, NS, EFFMULS, NREP, CFG) estimates the
+%   Uses the paper's exact model -- fitlme(tbl, 'behavior ~ stim + day +
-%   power to detect the day x tDCS interaction of the phased LME
+%   stim:day + (1|rat)') -- throughout. The fitted model for MERGEKEY over the
-%   (success ~ dayp*tDCS + (1|subject)). The fitted model for MERGEKEY over the
+%   WINDOW = [lo hi] day range (Box-B2 = stim = 1 vs Box-A2 = stim = 0) is the
-%   PHASE = [lo hi] window is used as ground truth (its fixed effects, subject
+%   ground truth (its fixed effects, rat random-intercept SD, residual SD);
-%   random-intercept SD, and residual SD); NREP datasets are simulated at each
+%   NREP datasets are simulated at each rats-per-group in NS, for each
-%   subjects-per-group in NS, for each true-effect multiplier in EFFMULS (e.g.
+%   true-effect multiplier in EFFMULS (e.g. [1 0.5] = observed and half the
-%   [1 0.5] = observed and half the observed interaction). Each dataset is
+%   observed stim:day slope). Each dataset is scored at alpha = 0.05 two ways:
-%   scored two ways at alpha = 0.05:
+%     - per-animal (cluster-honest): per-rat behavior~day slope, Welch t (stim)
-%     - per-animal (cluster-honest): per-subject slope, Welch t (Box-B2 vs A2)
+%     - LME: the fitlme stim:day interaction p (observation-level DF)
-%     - LME: the fitlme dayp:tDCS interaction p (observation-level DF)
+%   PW is a table (effMul, N, powerPerAnimal, powerLME); also printed. A fixed
-%   PW is a table (effMul, N, powerPerAnimal, powerLME); it is also printed.
+%   RNG seed makes the estimate reproducible.
-%   A fixed RNG seed makes the estimate reproducible.
 %
-%   Defaults: PHASE=[0 5], NS=[3 5 8 12 16 24 30], EFFMULS=[1 0.5], NREP=200.
+%   Defaults: WINDOW=[0 5] (early phase), NS=[3 5 8 12 16 24 30],
+%   EFFMULS=[1 0.5], NREP=200.
 
 if nargin < 1 || isempty(mergeKey); mergeKey = 'mergeA2'; end
-if nargin < 2 || isempty(phase);    phase = [0 5]; end
+if nargin < 2 || isempty(window);   window = [0 5]; end
 if nargin < 3 || isempty(Ns);       Ns = [3 5 8 12 16 24 30]; end
 if nargin < 4 || isempty(effMuls);  effMuls = [1 0.5]; end
 if nargin < 5 || isempty(nrep);     nrep = 200; end
 if nargin < 6 || isempty(cfg);      cfg = tdcs_config(); end
 
 warnState = warning('off', 'all'); cleanupW = onCleanup(@() warning(warnState)); %#ok<NASGU>
-rng(1);   % reproducible
+rng(1);
+FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
 
-% Ground truth = fitted phased LME.
+% Ground truth = fitted paper model over the window.
 Sfull = tdcs_scenario_data([mergeKey '_full']);
-T = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
+A = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
-Tp = T(T.day >= phase(1) & T.day <= phase(2), :);
+A = A(A.day >= window(1) & A.day <= window(2), :);
-Tp.dayp = Tp.day - phase(1);
+tbl0 = table();
-Tp.tDCS = double(Tp.group == cfg.anchorHigh);
+tbl0.behavior = A.success;
-lme = fitlme(Tp, 'success ~ dayp*tDCS + (1|subject)');
+tbl0.day = A.day - window(1);      % 0-based within the window
+tbl0.stim = double(A.group == cfg.anchorHigh);
+tbl0.rat = A.subject;
+lme = fitlme(tbl0, FORMULA);
 cn = lme.CoefficientNames; be = lme.fixedEffects;
-b0 = be(strcmp(cn,'(Intercept)')); bDay = be(strcmp(cn,'dayp'));
+b0 = be(strcmp(cn,'(Intercept)')); bStim = be(strcmp(cn,'stim'));
-bT = be(strcmp(cn,'tDCS'));        bInt = be(strcmp(cn,'dayp:tDCS'));
+bDay = be(strcmp(cn,'day'));       bInt = be(strcmp(cn,'day:stim'));
-psi = covarianceParameters(lme); sSub = sqrt(psi{1}); sRes = sqrt(lme.MSE);
+psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
-days = (phase(1):phase(2))' - phase(1);
+days = (window(1):window(2))' - window(1);
 
-fprintf('Power sim: truth=%s phase %d-%d | interaction=%.2f subjSD=%.2f resSD=%.2f | nrep=%d\n', ...
+fprintf('Power sim (paper formula): truth=%s window %d-%d | stim:day=%.2f ratSD=%.2f resSD=%.2f | nrep=%d\n', ...
-    mergeKey, phase(1), phase(2), bInt, sSub, sRes, nrep);
+    mergeKey, window(1), window(2), bInt, sRat, sRes, nrep);
 
 rows = {};
 for eMul = effMuls
     bI = bInt * eMul;
-    fprintf('\n true interaction = %+.2f (%.0f%% of observed)\n', bI, eMul*100);
+    fprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul*100);
     fprintf('  %-8s | per-animal power | LME power\n', 'N/group');
     for N = Ns
         sigPA = 0; sigL = 0;
         for r = 1:nrep
-            tbl = localSim(N, days, b0, bDay, bT, bI, sSub, sRes);
+            tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
             sigPA = sigPA + (localPerAnimalP(tbl) < 0.05);
             try
-                m = fitlme(tbl, 'success ~ dayp*tDCS + (1|subject)');
+                m = fitlme(tbl, FORMULA);
                 Cm = m.Coefficients;
-                sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'dayp:tDCS')) < 0.05);
+                sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'day:stim')) < 0.05);
             catch
             end
         end
@@ -67,31 +71,31 @@ PW = cell2table(rows, 'VariableNames', {'effMul','N','powerPerAnimal','powerLME'
 
 end
 
-function tbl = localSim(N, days, b0, bDay, bT, bI, sSub, sRes)
+function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
 nd = numel(days); rows = 2*N*nd;
-subj = strings(rows,1); dayp = zeros(rows,1); tDCS = zeros(rows,1); success = zeros(rows,1);
+rat = strings(rows,1); day = zeros(rows,1); stim = zeros(rows,1); behavior = zeros(rows,1);
 k = 0;
 for g = 0:1
     for sIdx = 1:N
-        re = sSub * randn;
+        re = sRat * randn;
-        sid = sprintf('g%d_s%d', g, sIdx);
+        rid = sprintf('g%d_r%d', g, sIdx);
         for d = 1:nd
             k = k+1;
-            subj(k) = sid; dayp(k) = days(d); tDCS(k) = g;
+            rat(k) = rid; day(k) = days(d); stim(k) = g;
-            success(k) = b0 + bDay*days(d) + bT*g + bI*days(d)*g + re + sRes*randn;
+            behavior(k) = b0 + bStim*g + bDay*days(d) + bI*days(d)*g + re + sRes*randn;
         end
     end
 end
-tbl = table(categorical(subj), dayp, tDCS, success, ...
+tbl = table(categorical(rat), day, stim, behavior, ...
-    'VariableNames', {'subject','dayp','tDCS','success'});
+    'VariableNames', {'rat','day','stim','behavior'});
 end
 
 function p = localPerAnimalP(tbl)
-subs = unique(tbl.subject); sl = zeros(numel(subs),1); gr = zeros(numel(subs),1);
+rats = unique(tbl.rat); sl = zeros(numel(rats),1); gr = zeros(numel(rats),1);
-for i = 1:numel(subs)
+for i = 1:numel(rats)
-    r = tbl.subject == subs(i);
+    r = tbl.rat == rats(i);
-    c = polyfit(tbl.dayp(r), tbl.success(r), 1); sl(i) = c(1);
+    c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
-    gr(i) = tbl.tDCS(find(r,1));
+    gr(i) = tbl.stim(find(r,1));
 end
 [~, p] = ttest2(sl(gr==1), sl(gr==0), 'Vartype', 'unequal');
 end

analysis/matlab/tests/tLme.m

@@ -59,6 +59,18 @@ classdef tLme < matlab.unittest.TestCase
             testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario');
         end
 
+        function testPaperFormulaReplicatesInteraction(testCase)
+            % The paper's exact formula (behavior ~ stim + day + stim:day +
+            % (1|rat)) on mergeA2 reproduces the paper's interaction
+            % (paper F(1)=7.12, p=0.008) -- same fit as the day*tDCS form.
+            cfg = tdcs_config();
+            evalc('L = tdcs_paper_lme(''mergeA2'', 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.verifyLessThan(L.day.p, 1e-6);   % strong learning
+        end
+
         function testD0_13WindowFairCoverageAndParallel(testCase)
             % The 0-13 window caps both anchor groups at day 13 (equal coverage),
             % unlike the full window where Box-A2 stops at 13 but Box-B2 runs on.