From 7b5c27df8c8e68c894788985b52282c2703c1f77 Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Mon, 20 Jul 2026 12:40:26 -0400 Subject: [PATCH] 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 --- analysis/matlab/README.md | 13 +++++ analysis/matlab/run_all.m | 3 +- analysis/matlab/tdcs_glm.m | 10 ++++ analysis/matlab/tdcs_paper_lme.m | 91 ++++++++++++++++++++++++++++++ analysis/matlab/tdcs_power_sim.m | 96 +++++++++++++++++--------------- analysis/matlab/tests/tLme.m | 12 ++++ 6 files changed, 178 insertions(+), 47 deletions(-) create mode 100644 analysis/matlab/tdcs_paper_lme.m diff --git a/analysis/matlab/README.md b/analysis/matlab/README.md index 5efeb22..c0290da 100644 --- a/analysis/matlab/README.md +++ b/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/_.csv` (9), `anchors/_B2vsA2.csv` (3, with the binary tDCS factor), and `phases/_phase_-.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. diff --git a/analysis/matlab/run_all.m b/analysis/matlab/run_all.m index b013bee..5441a4f 100644 --- a/analysis/matlab/run_all.m +++ b/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}); diff --git a/analysis/matlab/tdcs_glm.m b/analysis/matlab/tdcs_glm.m index 91a8b2d..eb39bda 100644 --- a/analysis/matlab/tdcs_glm.m +++ b/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: ' ... diff --git a/analysis/matlab/tdcs_paper_lme.m b/analysis/matlab/tdcs_paper_lme.m new file mode 100644 index 0000000..da5bbd8 --- /dev/null +++ b/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_.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 +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 diff --git a/analysis/matlab/tdcs_power_sim.m b/analysis/matlab/tdcs_power_sim.m index a72f7ef..486ccba 100644 --- a/analysis/matlab/tdcs_power_sim.m +++ b/analysis/matlab/tdcs_power_sim.m @@ -1,60 +1,64 @@ -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. +function PW = tdcs_power_sim(mergeKey, window, Ns, effMuls, nrep, cfg) +%TDCS_POWER_SIM Monte-Carlo power for the paper's stim x day interaction. +% Uses the paper's exact model -- fitlme(tbl, 'behavior ~ stim + day + +% stim:day + (1|rat)') -- throughout. The fitted model for MERGEKEY over the +% WINDOW = [lo hi] day range (Box-B2 = stim = 1 vs Box-A2 = stim = 0) is the +% ground truth (its fixed effects, rat random-intercept SD, residual SD); +% NREP datasets are simulated at each rats-per-group in NS, for each +% true-effect multiplier in EFFMULS (e.g. [1 0.5] = observed and half the +% observed stim:day slope). Each dataset is scored at alpha = 0.05 two ways: +% - per-animal (cluster-honest): per-rat behavior~day slope, Welch t (stim) +% - LME: the fitlme stim:day interaction p (observation-level DF) +% PW is a table (effMul, N, powerPerAnimal, powerLME); 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. +% 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 -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}), :); -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)'); +A = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :); +A = A(A.day >= window(1) & A.day <= window(2), :); +tbl0 = table(); +tbl0.behavior = A.success; +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')); -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); +b0 = be(strcmp(cn,'(Intercept)')); bStim = be(strcmp(cn,'stim')); +bDay = be(strcmp(cn,'day')); bInt = be(strcmp(cn,'day:stim')); +psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); +days = (window(1):window(2))' - window(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); +fprintf('Power sim (paper formula): truth=%s window %d-%d | stim:day=%.2f ratSD=%.2f resSD=%.2f | nrep=%d\n', ... + 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; - sid = sprintf('g%d_s%d', g, sIdx); + re = sRat * randn; + 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; - success(k) = b0 + bDay*days(d) + bT*g + bI*days(d)*g + re + sRes*randn; + rat(k) = rid; day(k) = days(d); stim(k) = g; + 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, ... - 'VariableNames', {'subject','dayp','tDCS','success'}); +tbl = table(categorical(rat), day, stim, behavior, ... + 'VariableNames', {'rat','day','stim','behavior'}); 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)); +rats = unique(tbl.rat); sl = zeros(numel(rats),1); gr = zeros(numel(rats),1); +for i = 1:numel(rats) + r = tbl.rat == rats(i); + c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1); + gr(i) = tbl.stim(find(r,1)); end [~, p] = ttest2(sl(gr==1), sl(gr==0), 'Vartype', 'unequal'); end diff --git a/analysis/matlab/tests/tLme.m b/analysis/matlab/tests/tLme.m index 16928d0..4c671fd 100644 --- a/analysis/matlab/tests/tLme.m +++ b/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.