feat(matlab): phased days x tDCS LME + Monte-Carlo power analysis

Adds tdcs_phase_lme (same LME refit within learning phases 0-5/6-10/6-13, with
per-phase slopes + interaction 95% CI) as phase_* switch cases, and
tdcs_power_sim (Monte-Carlo power for the early-phase interaction, scored by the
cluster-honest per-animal test and the LME test) + run_power. Findings: the
Box-B2 faster-acquisition signal is in the early phase; at n=3-5/group honest
power is 0.3-0.7 even at the observed effect (need ~8/group if effect is as
observed, ~20/group if half). Adds tPhasePower tests. Suite 35/35.

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

analysis/matlab/README.md

@@ -107,3 +107,24 @@ sensitivity sections — `report_anchor_only` (a two-anchor-only refit) and
 `report_mixed` (a Bayesian MAP random-intercept model) — are intentionally
 **not** ported: the anchor-only refit is closely tracked by the main-model H2
 term, and MATLAB's native `fitglme` already *is* the random-intercept model.
+
+## Phased learning analysis and power
+
+- **Phased days × tDCS LME** — the same `success ~ day*tDCS + (1|subject)` model refit within
+  learning phases (0–5, 6–10, 6–13), Box-B2 vs Box-A2, with per-phase slopes and the interaction
+  95% CI:
+  ```
+  matlab -batch "tdcs_glm('phase_mergeA2')"   # or phase_unmerged / phase_mergeB2
+  ```
+  The Box-B2 faster-acquisition signal sits in the early (0–5) phase; late phases converge. The
+  interaction p/CI use fitlme observation-level DF and are anticonservative at these small n
+  (see power analysis).
+
+- **Monte-Carlo power** — power to detect the early-phase interaction, scored with the
+  cluster-honest per-animal test and the LME test, across subjects-per-group and effect sizes:
+  ```
+  matlab -batch "run_power"          # default: mergeA2, phase 0–5
+  ```
+  Result: at our n (3–5/group) honest power is ~0.3–0.7 even at the observed effect; ~8/group
+  reaches ~80–90% if the effect is as large as observed, ~20/group if it is half that. The LME
+  test is anticonservative only at small n and converges to the honest test by n≈12.

analysis/matlab/run_all.m

@@ -8,7 +8,8 @@ scenarios = {'unmerged_full', 'unmerged_d0_10', ...
              'lme_unmerged_full', 'lme_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'};
+             'lme_unmerged_d0_13', 'lme_mergeA2_d0_13', 'lme_mergeB2_d0_13', ...
+             'phase_unmerged', 'phase_mergeA2', 'phase_mergeB2'};
 
 for i = 1:numel(scenarios)
     fprintf('\n\n### Running scenario: %s ###\n', scenarios{i});

analysis/matlab/run_power.m

@@ -0,0 +1,8 @@
+%RUN_POWER Monte-Carlo power analysis for the early-phase days x tDCS interaction.
+%   Usage: matlab -batch "run_power"
+%   Uses the fitted mergeA2 early-phase (0-5) LME as ground truth and reports
+%   power (cluster-honest per-animal test and the LME test) across a range of
+%   subjects-per-group, at the observed and half-observed interaction sizes.
+%   See tdcs_power_sim for parameters.
+
+tdcs_power_sim('mergeA2', [0 5]);

analysis/matlab/tdcs_glm.m

@@ -97,12 +97,23 @@ switch scenario
         L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_13'), cfg);
         tdcs_lme_report(L, 'lme_mergeB2_d0_13', cfg);
 
+    % --- Phased days x tDCS LME (fast/slow learning phases), Box-B2 vs Box-A2 ---
+    case 'phase_unmerged'
+        tdcs_phase_lme('unmerged', cfg);
+
+    case 'phase_mergeA2'
+        tdcs_phase_lme('mergeA2', cfg);
+
+    case 'phase_mergeB2'
+        tdcs_phase_lme('mergeB2', 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, and their lme_* variants ' ...
+            'mergeB2_full, mergeB2_d0_10, their lme_* variants ' ...
-            '(lme_unmerged_full, ..., lme_mergeB2_d0_10).'], scenario);
+            '(lme_unmerged_full, ..., lme_mergeB2_d0_13), and the phased LMEs ' ...
+            'phase_unmerged, phase_mergeA2, phase_mergeB2.'], scenario);
 end
 
 end

analysis/matlab/tdcs_phase_lme.m

@@ -0,0 +1,81 @@
+function P = tdcs_phase_lme(mergeKey, cfg, phases)
+%TDCS_PHASE_LME Refit the days x tDCS LME within learning phases (Box-B2 vs A2).
+%   P = TDCS_PHASE_LME(MERGEKEY, CFG) fits the same linear mixed model as the
+%   lme_* scenarios -- success ~ dayp*tDCS + (1|subject) -- SEPARATELY within
+%   each learning phase, for the Box-B2 (tDCS=1) vs Box-A2 (tDCS=0) subset of
+%   the MERGEKEY grouping ('unmerged' | 'mergeA2' | 'mergeB2'). `dayp` is `day`
+%   centered at each phase's start, so the tDCS main effect reads as the group
+%   difference on the phase's first day (the interaction and slopes are
+%   invariant to this centering). Prints a per-phase table and writes
+%   results/phase_<MERGEKEY>.txt.
+%
+%   P = TDCS_PHASE_LME(MERGEKEY, CFG, PHASES) uses custom phase windows
+%   (default {[0 5],[6 10],[6 13]}), each a [loDay hiDay] pair.
+%
+%   P.phases{k} has: .phase .nSub .a2Slope .b2Slope .dayP (learning) .tDCSlevelP
+%   (between-group level at phase start) .interP .interEst .interCI (the
+%   between-group learning-rate difference and its 95% CI).
+
+if nargin < 3 || isempty(phases)
+    phases = {[0 5], [6 10], [6 13]};
+end
+
+Sfull = tdcs_scenario_data([mergeKey '_full']);
+T = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
+
+P = struct('mergeKey', mergeKey, 'phases', {cell(1, numel(phases))});
+
+bar = repmat('=', 1, 78);
+s = sprintf('%s\n', bar);
+s = [s sprintf('PHASED days x tDCS LME -- %s (Box-B2 vs Box-A2)\n', mergeKey)];
+s = [s sprintf('%s\n', bar)];
+s = [s sprintf('model per phase: success ~ dayp*tDCS + (1|subject)  [dayp = day - phaseStart]\n\n')];
+s = [s sprintf('%-9s N  A2slope B2slope  day p     tDCS(lvl) p  interaction p  slopeDiff [95%% CI]\n', 'phase')];
+s = [s sprintf('%s\n', repmat('-', 1, 92))];
+
+for k = 1:numel(phases)
+    ph = phases{k};
+    Tp = T(T.day >= ph(1) & T.day <= ph(2), :);
+    Tp.dayp = Tp.day - ph(1);
+    Tp.tDCS = double(Tp.group == cfg.anchorHigh);
+
+    lme = fitlme(Tp, 'success ~ dayp*tDCS + (1|subject)');
+    C = lme.Coefficients; A = anova(lme); ci = coefCI(lme);
+    ii = strcmp(C.Name, 'dayp:tDCS');
+    di = strcmp(C.Name, 'dayp');
+
+    e = struct();
+    e.phase = ph;
+    e.nSub = numel(unique(Tp.subject));
+    e.a2Slope = C.Estimate(di);
+    e.b2Slope = C.Estimate(di) + C.Estimate(ii);
+    e.dayP = A.pValue(strcmp(A.Term, 'dayp'));
+    e.tDCSlevelP = A.pValue(strcmp(A.Term, 'tDCS'));
+    e.interP = A.pValue(strcmp(A.Term, 'dayp:tDCS'));
+    e.interEst = C.Estimate(ii);
+    e.interCI = ci(ii, :);
+    P.phases{k} = e;
+
+    s = [s sprintf('%d-%-6d %d  %6.2f  %6.2f  %-9.2g %-11.3f  %-13.3f %+.2f [%+.2f, %+.2f]\n', ...
+        ph(1), ph(2), e.nSub, e.a2Slope, e.b2Slope, e.dayP, e.tDCSlevelP, ...
+        e.interP, e.interEst, e.interCI(1), e.interCI(2))]; %#ok<AGROW>
+end
+
+s = [s sprintf(['\nNote: the interaction p (and CI) use fitlme observation-level DF and are\n' ...
+    'ANTICONSERVATIVE at these small subject counts (see tdcs_power_sim). The early\n' ...
+    'phase carries the Box-B2 faster-acquisition signal; late phases converge.\n'])];
+
+fprintf('%s', s);
+localWrite(['phase_' mergeKey], s);
+
+end
+
+function localWrite(name, s)
+thisDir = fileparts(mfilename('fullpath'));
+resDir = fullfile(thisDir, 'results');
+if ~exist(resDir, 'dir'); mkdir(resDir); end
+fid = fopen(fullfile(resDir, [name '.txt']), 'w');
+if fid < 0; error('tdcs_phase_lme:fopen', 'Cannot open results file for "%s".', name); end
+cleanup = onCleanup(@() fclose(fid)); %#ok<NASGU>
+fprintf(fid, '%s', s);
+end

analysis/matlab/tdcs_power_sim.m

@@ -0,0 +1,97 @@
+function PW = tdcs_power_sim(mergeKey, phase, Ns, effMuls, nrep, cfg)
+%TDCS_POWER_SIM Monte-Carlo power for the phased days x tDCS interaction.
+%   PW = TDCS_POWER_SIM(MERGEKEY, PHASE, NS, EFFMULS, NREP, CFG) estimates the
+%   power to detect the day x tDCS interaction of the phased LME
+%   (success ~ dayp*tDCS + (1|subject)). The fitted model for MERGEKEY over the
+%   PHASE = [lo hi] window is used as ground truth (its fixed effects, subject
+%   random-intercept SD, and residual SD); NREP datasets are simulated at each
+%   subjects-per-group in NS, for each true-effect multiplier in EFFMULS (e.g.
+%   [1 0.5] = observed and half the observed interaction). Each dataset is
+%   scored two ways at alpha = 0.05:
+%     - per-animal (cluster-honest): per-subject slope, Welch t (Box-B2 vs A2)
+%     - LME: the fitlme dayp:tDCS interaction p (observation-level DF)
+%   PW is a table (effMul, N, powerPerAnimal, powerLME); it is also printed.
+%   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.
+
+if nargin < 1 || isempty(mergeKey); mergeKey = 'mergeA2'; end
+if nargin < 2 || isempty(phase);    phase = [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
+
+% Ground truth = fitted phased LME.
+Sfull = tdcs_scenario_data([mergeKey '_full']);
+T = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
+Tp = T(T.day >= phase(1) & T.day <= phase(2), :);
+Tp.dayp = Tp.day - phase(1);
+Tp.tDCS = double(Tp.group == cfg.anchorHigh);
+lme = fitlme(Tp, 'success ~ dayp*tDCS + (1|subject)');
+cn = lme.CoefficientNames; be = lme.fixedEffects;
+b0 = be(strcmp(cn,'(Intercept)')); bDay = be(strcmp(cn,'dayp'));
+bT = be(strcmp(cn,'tDCS'));        bInt = be(strcmp(cn,'dayp:tDCS'));
+psi = covarianceParameters(lme); sSub = sqrt(psi{1}); sRes = sqrt(lme.MSE);
+days = (phase(1):phase(2))' - phase(1);
+
+fprintf('Power sim: truth=%s phase %d-%d | interaction=%.2f subjSD=%.2f resSD=%.2f | nrep=%d\n', ...
+    mergeKey, phase(1), phase(2), bInt, sSub, sRes, nrep);
+
+rows = {};
+for eMul = effMuls
+    bI = bInt * eMul;
+    fprintf('\n true 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);
+            sigPA = sigPA + (localPerAnimalP(tbl) < 0.05);
+            try
+                m = fitlme(tbl, 'success ~ dayp*tDCS + (1|subject)');
+                Cm = m.Coefficients;
+                sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'dayp:tDCS')) < 0.05);
+            catch
+            end
+        end
+        pPA = sigPA / nrep; pL = sigL / nrep;
+        fprintf('  %-8d | %5.2f            | %5.2f\n', N, pPA, pL);
+        rows(end+1, :) = {eMul, N, pPA, pL}; %#ok<AGROW>
+    end
+end
+PW = cell2table(rows, 'VariableNames', {'effMul','N','powerPerAnimal','powerLME'});
+
+end
+
+function tbl = localSim(N, days, b0, bDay, bT, bI, sSub, sRes)
+nd = numel(days); rows = 2*N*nd;
+subj = strings(rows,1); dayp = zeros(rows,1); tDCS = zeros(rows,1); success = zeros(rows,1);
+k = 0;
+for g = 0:1
+    for sIdx = 1:N
+        re = sSub * randn;
+        sid = sprintf('g%d_s%d', g, sIdx);
+        for d = 1:nd
+            k = k+1;
+            subj(k) = sid; dayp(k) = days(d); tDCS(k) = g;
+            success(k) = b0 + bDay*days(d) + bT*g + bI*days(d)*g + re + sRes*randn;
+        end
+    end
+end
+tbl = table(categorical(subj), dayp, tDCS, success, ...
+    'VariableNames', {'subject','dayp','tDCS','success'});
+end
+
+function p = localPerAnimalP(tbl)
+subs = unique(tbl.subject); sl = zeros(numel(subs),1); gr = zeros(numel(subs),1);
+for i = 1:numel(subs)
+    r = tbl.subject == subs(i);
+    c = polyfit(tbl.dayp(r), tbl.success(r), 1); sl(i) = c(1);
+    gr(i) = tbl.tDCS(find(r,1));
+end
+[~, p] = ttest2(sl(gr==1), sl(gr==0), 'Vartype', 'unequal');
+end

analysis/matlab/tests/tPhasePower.m

@@ -0,0 +1,41 @@
+classdef tPhasePower < matlab.unittest.TestCase
+    %TPHASEPOWER Tests for the phased days x tDCS LME (tdcs_phase_lme) and the
+    %   Monte-Carlo power analysis (tdcs_power_sim).
+
+    methods (Test)
+
+        function testPhaseEarlyFasterAndSlowDown(testCase)
+            % Early phase: Box-B2 slope > Box-A2 slope (faster acquisition) and
+            % strong learning; late phase: shallower slopes (fast->slow).
+            cfg = tdcs_config();
+            evalc('P = tdcs_phase_lme(''mergeA2'', cfg);');   % run in-workspace, suppress print
+            early = P.phases{1};   % [0 5]
+            late = P.phases{2};    % [6 10]
+            testCase.verifyEqual(early.phase, [0 5]);
+            testCase.verifyGreaterThan(early.b2Slope, early.a2Slope);   % B2 faster early
+            testCase.verifyLessThan(early.dayP, 0.01);                  % strong early learning
+            testCase.verifyGreaterThan(early.b2Slope, late.b2Slope);    % slow-down (fast->slow)
+            testCase.verifyGreaterThan(early.interEst, 0);              % positive interaction
+        end
+
+        function testPhaseSwitchWritesResults(testCase)
+            here = fileparts(fileparts(mfilename('fullpath')));  % analysis/matlab
+            outFile = fullfile(here, 'results', 'phase_mergeA2.txt');
+            if exist(outFile, 'file'); delete(outFile); end
+            evalc("tdcs_glm('phase_mergeA2')");
+            testCase.verifyTrue(exist(outFile, 'file') == 2);
+        end
+
+        function testPowerMonotonicAndBounded(testCase)
+            % Small, fast config: power rises with N and stays in [0,1].
+            cfg = tdcs_config();
+            PW = tdcs_power_sim('mergeA2', [0 5], [4 24], 1.0, 30, cfg);
+            testCase.verifyTrue(all(PW.powerPerAnimal >= 0 & PW.powerPerAnimal <= 1));
+            testCase.verifyTrue(all(PW.powerLME >= 0 & PW.powerLME <= 1));
+            pSmall = PW.powerPerAnimal(PW.N == 4);
+            pLarge = PW.powerPerAnimal(PW.N == 24);
+            testCase.verifyGreaterThan(pLarge, pSmall);   % more subjects -> more power
+        end
+
+    end
+end