From e730aa00209d517d9a684ee7df6b5bd0a8c45e19 Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Mon, 20 Jul 2026 11:06:51 -0400 Subject: [PATCH] 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 --- analysis/matlab/README.md | 21 +++++++ analysis/matlab/run_all.m | 3 +- analysis/matlab/run_power.m | 8 +++ analysis/matlab/tdcs_glm.m | 15 ++++- analysis/matlab/tdcs_phase_lme.m | 81 ++++++++++++++++++++++++ analysis/matlab/tdcs_power_sim.m | 97 +++++++++++++++++++++++++++++ analysis/matlab/tests/tPhasePower.m | 41 ++++++++++++ 7 files changed, 263 insertions(+), 3 deletions(-) create mode 100644 analysis/matlab/run_power.m create mode 100644 analysis/matlab/tdcs_phase_lme.m create mode 100644 analysis/matlab/tdcs_power_sim.m create mode 100644 analysis/matlab/tests/tPhasePower.m diff --git a/analysis/matlab/README.md b/analysis/matlab/README.md index 32af9d0..c2aec5e 100644 --- a/analysis/matlab/README.md +++ b/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. diff --git a/analysis/matlab/run_all.m b/analysis/matlab/run_all.m index 88f5873..b013bee 100644 --- a/analysis/matlab/run_all.m +++ b/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}); diff --git a/analysis/matlab/run_power.m b/analysis/matlab/run_power.m new file mode 100644 index 0000000..ffa2668 --- /dev/null +++ b/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]); diff --git a/analysis/matlab/tdcs_glm.m b/analysis/matlab/tdcs_glm.m index 65e0ca6..91a8b2d 100644 --- a/analysis/matlab/tdcs_glm.m +++ b/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 ' ... - '(lme_unmerged_full, ..., lme_mergeB2_d0_10).'], scenario); + 'mergeB2_full, mergeB2_d0_10, their lme_* variants ' ... + '(lme_unmerged_full, ..., lme_mergeB2_d0_13), and the phased LMEs ' ... + 'phase_unmerged, phase_mergeA2, phase_mergeB2.'], scenario); end end diff --git a/analysis/matlab/tdcs_phase_lme.m b/analysis/matlab/tdcs_phase_lme.m new file mode 100644 index 0000000..4e08910 --- /dev/null +++ b/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_.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 +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 +fprintf(fid, '%s', s); +end diff --git a/analysis/matlab/tdcs_power_sim.m b/analysis/matlab/tdcs_power_sim.m new file mode 100644 index 0000000..a72f7ef --- /dev/null +++ b/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 +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 + 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 diff --git a/analysis/matlab/tests/tPhasePower.m b/analysis/matlab/tests/tPhasePower.m new file mode 100644 index 0000000..94e431c --- /dev/null +++ b/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