Experiments DB Dev
8a18c894dd
Every fitlme-based report (lme_*, paper_*, phase_*, and the variations' analyze.m) now shows, per effect: residual-DF p, Satterthwaite-DF p, and -- for the interaction -- an HONEST test from a per-animal random-SLOPE model (day|rat), whose Satterthwaite DF collapses toward the animal count. New: tdcs_random_slope_interaction.m (shared helper). Wired into tdcs_lme, tdcs_paper_lme, tdcs_phase_lme, variation_analyze; SUMMARY.csv gains interaction_p_satt / interaction_p_rs. Regenerated all results/, variations/, matched-effort outputs. Key point this surfaces: Satterthwaite ~= residual on the random-INTERCEPT model (the slope's error is at session level), so it does NOT fix pseudoreplication; the random-slope model does. Effect: full-range mergeA2 interaction 0.009 -> 0.75 (collapses); unmerge_d0_5 0.015 -> 0.13 (n.s.); the pooled-control early windows survive honestly (naive_a2_d0_5 0.001 -> 0.028; naive_boxa_d0_5 0.003 -> 0.036). Suite 42/42. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
177 lines
9.2 KiB
Matlab
177 lines
9.2 KiB
Matlab
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
|
|
|
|
function testPaperMergeNaiveWindowedFairCoverage(testCase)
|
|
% The windowed paper_mergeNaive_* runs cap both anchor arms at the
|
|
% window end (equal coverage -> no confound warning), write their own
|
|
% results file, and embed the full fitlme summary. Unlike mergeA2,
|
|
% the pooled-naive control keeps the interaction significant in the
|
|
% fair d0_13 window (~ the paper's F=7.12, p=0.008).
|
|
here = fileparts(fileparts(mfilename('fullpath'))); % analysis/matlab
|
|
resDir = fullfile(here, 'results');
|
|
for w = {'d0_10', 'd0_13'}
|
|
name = ['paper_mergeNaive_' w{1}];
|
|
evalc(sprintf('tdcs_glm(''%s'')', name));
|
|
txt = fileread(fullfile(resDir, [name '.txt']));
|
|
testCase.verifyTrue(contains(txt, 'Equal day coverage'), ...
|
|
sprintf('%s: expected equal (fair) day coverage', name));
|
|
testCase.verifyFalse(contains(txt, '** WARNING'), ...
|
|
sprintf('%s: fair window should not warn', name));
|
|
testCase.verifyEqual(count(txt, 'Fixed effects coefficients (95% CIs):'), 1);
|
|
end
|
|
% Fair-window d0_13 interaction reproduces the paper (~F=7.1, p~0.009).
|
|
L = tdcs_paper_lme('mergeNaive', tdcs_config(), 'd0_13');
|
|
testCase.verifyEqual(L.maxDayStim, L.maxDayCtrl); % equal coverage
|
|
testCase.verifyLessThan(L.interaction.p, 0.05);
|
|
testCase.verifyGreaterThan(L.interaction.estimate, 0);
|
|
end
|
|
|
|
function testPaperLmeBadWindowErrors(testCase)
|
|
testCase.verifyError(@() tdcs_paper_lme('mergeNaive', tdcs_config(), 'bogus'), ...
|
|
'tdcs_paper_lme:badWindow');
|
|
end
|
|
|
|
function testReportsIncludeFullModelSummary(testCase)
|
|
% Every scenario prints the full native fitlme summary. The lme_*,
|
|
% paper_*, and phase_* reports each embed disp(model), whose
|
|
% 'Fixed effects coefficients' header is the marker (phase_* has
|
|
% one per phase -> 3). Markup must be stripped.
|
|
here = fileparts(fileparts(mfilename('fullpath'))); % analysis/matlab
|
|
resDir = fullfile(here, 'results');
|
|
cases = struct( ...
|
|
'lme_mergeA2_full', 1, ...
|
|
'paper_mergeA2', 1, ...
|
|
'phase_mergeA2', 3);
|
|
for sc = fieldnames(cases)'
|
|
name = sc{1};
|
|
evalc(sprintf('tdcs_glm(''%s'')', name));
|
|
txt = fileread(fullfile(resDir, [name '.txt']));
|
|
testCase.verifyEqual(count(txt, 'Fixed effects coefficients (95% CIs):'), ...
|
|
cases.(name), sprintf('%s: missing full model summary', name));
|
|
testCase.verifyFalse(contains(txt, '<strong>'), ...
|
|
sprintf('%s: markup not stripped', name));
|
|
end
|
|
end
|
|
|
|
function testMergeNaiveGrouping(testCase)
|
|
% mergeNaive: control = Box-A2 + Naive (7), tDCS = Box-B2 + Right (4),
|
|
% Box-A unchanged (1); Naive is relabeled away.
|
|
[S, ~] = tdcs_scenario_data('mergeNaive_full');
|
|
sg = unique(S(:, {'subject','group'}));
|
|
n = @(g) sum(sg.group == g);
|
|
testCase.verifyEqual(n('Electrode-Box-A2'), 7);
|
|
testCase.verifyEqual(n('Electrode-Box-B2'), 4);
|
|
testCase.verifyEqual(n('Electrode-Box-A'), 1);
|
|
testCase.verifyFalse(any(sg.group == 'Naive'));
|
|
end
|
|
|
|
function testLmeReportsHaveSatterthwaiteAndHonest(testCase)
|
|
% Every fitlme report adds Satterthwaite DF and the honest
|
|
% random-slope interaction test.
|
|
cfg = tdcs_config();
|
|
L = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg);
|
|
testCase.verifyTrue(isfield(L.interaction, 'pSatt'));
|
|
testCase.verifyTrue(isfield(L.interaction, 'dfSatt'));
|
|
testCase.verifyTrue(isfield(L, 'interRS') && islogical(L.interRS.ok));
|
|
% The full-range mergeA2 interaction is significant under obs-level
|
|
% DF but NOT under the honest per-animal random slope (it collapses).
|
|
testCase.verifyLessThan(L.interaction.p, 0.05);
|
|
testCase.verifyTrue(L.interRS.ok);
|
|
testCase.verifyGreaterThan(L.interRS.p, 0.05);
|
|
% Report text carries the Satterthwaite column and honest block.
|
|
here = fileparts(fileparts(mfilename('fullpath')));
|
|
evalc("tdcs_glm('lme_mergeA2_full')");
|
|
txt = fileread(fullfile(here, 'results', 'lme_mergeA2_full.txt'));
|
|
testCase.verifyTrue(contains(txt, 'Satterthwaite'));
|
|
testCase.verifyTrue(contains(txt, 'HONEST LME'));
|
|
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.
|
|
cfg = tdcs_config();
|
|
Lfull = tdcs_lme(tdcs_scenario_data('mergeA2_full'), cfg);
|
|
L13 = tdcs_lme(tdcs_scenario_data('mergeA2_d0_13'), cfg);
|
|
% Full window: unequal coverage (A2 ends ~13, B2 runs later).
|
|
testCase.verifyGreaterThan(abs(Lfull.maxDayB - Lfull.maxDayA), 2);
|
|
% 0-13 window: equal coverage and the (truncation-driven) interaction
|
|
% is no longer significant -- the fair-window result.
|
|
testCase.verifyLessThanOrEqual(L13.maxDayA, 13);
|
|
testCase.verifyLessThanOrEqual(L13.maxDayB, 13);
|
|
testCase.verifyGreaterThan(L13.interaction.p, 0.05);
|
|
testCase.verifyLessThan(L13.day.p, 1e-6); % learning still strong
|
|
end
|
|
|
|
end
|
|
end
|