feat(matlab): add 0-13 day window (lme_*_d0_13) + coverage warning
Box-A2 has no data past ~day 13, so the full-range days x tDCS interaction extrapolates A2's slope and is confounded (spuriously significant, wrong sign). Add the _d0_13 window (equal anchor coverage) and a report warning when the two groups' day coverage is unequal. In the fair 0-13 window the interaction is n.s. with a POSITIVE slope diff and the Day-1 tDCS effect is n.s. -- matching the paper's 'equal on Day 1' direction (underpowered at n=8). Suite 32/32. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
@@ -7,7 +7,8 @@ scenarios = {'unmerged_full', 'unmerged_d0_10', ...
|
|||||||
'mergeB2_full', 'mergeB2_d0_10', ...
|
'mergeB2_full', 'mergeB2_d0_10', ...
|
||||||
'lme_unmerged_full', 'lme_unmerged_d0_10', ...
|
'lme_unmerged_full', 'lme_unmerged_d0_10', ...
|
||||||
'lme_mergeA2_full', 'lme_mergeA2_d0_10', ...
|
'lme_mergeA2_full', 'lme_mergeA2_d0_10', ...
|
||||||
'lme_mergeB2_full', 'lme_mergeB2_d0_10'};
|
'lme_mergeB2_full', 'lme_mergeB2_d0_10', ...
|
||||||
|
'lme_unmerged_d0_13', 'lme_mergeA2_d0_13', 'lme_mergeB2_d0_13'};
|
||||||
|
|
||||||
for i = 1:numel(scenarios)
|
for i = 1:numel(scenarios)
|
||||||
fprintf('\n\n### Running scenario: %s ###\n', scenarios{i});
|
fprintf('\n\n### Running scenario: %s ###\n', scenarios{i});
|
||||||
|
|||||||
@@ -10,9 +10,11 @@ function tdcs_glm(scenario)
|
|||||||
% mergeA2_full, mergeA2_d0_10,
|
% mergeA2_full, mergeA2_d0_10,
|
||||||
% mergeB2_full, mergeB2_d0_10
|
% mergeB2_full, mergeB2_d0_10
|
||||||
% or a linear "days x tDCS" mixed model (Box-B2 vs Box-A2, per the paper):
|
% or a linear "days x tDCS" mixed model (Box-B2 vs Box-A2, per the paper):
|
||||||
% lme_unmerged_full, lme_unmerged_d0_10,
|
% lme_unmerged_full, lme_unmerged_d0_10, lme_unmerged_d0_13,
|
||||||
% lme_mergeA2_full, lme_mergeA2_d0_10,
|
% lme_mergeA2_full, lme_mergeA2_d0_10, lme_mergeA2_d0_13,
|
||||||
% lme_mergeB2_full, lme_mergeB2_d0_10
|
% lme_mergeB2_full, lme_mergeB2_d0_10, lme_mergeB2_d0_13
|
||||||
|
% The _d0_13 window ends where Box-A2 runs out of data, giving both anchor
|
||||||
|
% groups equal day coverage (the fair window for the interaction).
|
||||||
%
|
%
|
||||||
% Each case below is a self-contained, independently invocable call of
|
% Each case below is a self-contained, independently invocable call of
|
||||||
% the same three-function pipeline; they are enumerated explicitly
|
% the same three-function pipeline; they are enumerated explicitly
|
||||||
@@ -83,6 +85,18 @@ switch scenario
|
|||||||
L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_10'), cfg);
|
L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_10'), cfg);
|
||||||
tdcs_lme_report(L, 'lme_mergeB2_d0_10', cfg);
|
tdcs_lme_report(L, 'lme_mergeB2_d0_10', cfg);
|
||||||
|
|
||||||
|
case 'lme_unmerged_d0_13'
|
||||||
|
L = tdcs_lme(tdcs_scenario_data('unmerged_d0_13'), cfg);
|
||||||
|
tdcs_lme_report(L, 'lme_unmerged_d0_13', cfg);
|
||||||
|
|
||||||
|
case 'lme_mergeA2_d0_13'
|
||||||
|
L = tdcs_lme(tdcs_scenario_data('mergeA2_d0_13'), cfg);
|
||||||
|
tdcs_lme_report(L, 'lme_mergeA2_d0_13', cfg);
|
||||||
|
|
||||||
|
case 'lme_mergeB2_d0_13'
|
||||||
|
L = tdcs_lme(tdcs_scenario_data('mergeB2_d0_13'), cfg);
|
||||||
|
tdcs_lme_report(L, 'lme_mergeB2_d0_13', cfg);
|
||||||
|
|
||||||
otherwise
|
otherwise
|
||||||
error('tdcs_glm:badScenario', ...
|
error('tdcs_glm:badScenario', ...
|
||||||
['Unrecognized scenario "%s". Valid scenarios are: ' ...
|
['Unrecognized scenario "%s". Valid scenarios are: ' ...
|
||||||
|
|||||||
@@ -36,6 +36,8 @@ A = anova(lme);
|
|||||||
L.lme = lme;
|
L.lme = lme;
|
||||||
L.nSubjects = numel(unique(T.subject));
|
L.nSubjects = numel(unique(T.subject));
|
||||||
L.nObs = height(T);
|
L.nObs = height(T);
|
||||||
|
L.maxDayA = max(T.day(T.group == cfg.anchorLow)); % last day Box-A2 has data
|
||||||
|
L.maxDayB = max(T.day(T.group == cfg.anchorHigh)); % last day Box-B2 has data
|
||||||
L.interaction = localTerm(C, A, 'day:tDCS');
|
L.interaction = localTerm(C, A, 'day:tDCS');
|
||||||
L.day = localTerm(C, A, 'day');
|
L.day = localTerm(C, A, 'day');
|
||||||
L.tDCS = localTerm(C, A, 'tDCS');
|
L.tDCS = localTerm(C, A, 'tDCS');
|
||||||
|
|||||||
@@ -12,8 +12,16 @@ s = [s sprintf('%s\n', bar)];
|
|||||||
s = [s sprintf('model: success ~ day * tDCS + (1|subject) [tDCS: %s = 1 vs %s = 0]\n', ...
|
s = [s sprintf('model: success ~ day * tDCS + (1|subject) [tDCS: %s = 1 vs %s = 0]\n', ...
|
||||||
cfg.anchorHigh, cfg.anchorLow)];
|
cfg.anchorHigh, cfg.anchorLow)];
|
||||||
s = [s sprintf('N = %d subjects, %d sessions\n', L.nSubjects, L.nObs)];
|
s = [s sprintf('N = %d subjects, %d sessions\n', L.nSubjects, L.nObs)];
|
||||||
|
s = [s sprintf('day coverage: Box-B2 0..%d, Box-A2 0..%d\n', L.maxDayB, L.maxDayA)];
|
||||||
s = [s sprintf('(day is raw and 0-indexed: our day 0 = the paper''s "Day 1", so the tDCS\n')];
|
s = [s sprintf('(day is raw and 0-indexed: our day 0 = the paper''s "Day 1", so the tDCS\n')];
|
||||||
s = [s sprintf(' main effect below is the group difference on Day 1 -- comparable to the paper.)\n\n')];
|
s = [s sprintf(' main effect below is the group difference on Day 1 -- comparable to the paper.)\n')];
|
||||||
|
if abs(L.maxDayB - L.maxDayA) > 2
|
||||||
|
s = [s sprintf(['** WARNING: the groups'' day coverage is UNEQUAL (differ by %d days). The\n' ...
|
||||||
|
' interaction/slope over this window EXTRAPOLATES the shorter group''s line and is\n' ...
|
||||||
|
' CONFOUNDED -- prefer the _d0_%d window (both groups have data throughout). **\n'], ...
|
||||||
|
abs(L.maxDayB - L.maxDayA), min(L.maxDayB, L.maxDayA))];
|
||||||
|
end
|
||||||
|
s = [s sprintf('\n')];
|
||||||
|
|
||||||
s = [s sprintf('%-27s %-14s %-13s %s\n', 'effect', 't (df)', 'F (df1)', 'p')];
|
s = [s sprintf('%-27s %-14s %-13s %s\n', 'effect', 't (df)', 'F (df1)', 'p')];
|
||||||
s = [s sprintf('%s\n', repmat('-', 1, 70))];
|
s = [s sprintf('%s\n', repmat('-', 1, 70))];
|
||||||
|
|||||||
@@ -92,7 +92,7 @@ meta.isUnmerged = strcmp(mergeKey, 'unmerged');
|
|||||||
end
|
end
|
||||||
|
|
||||||
function [mergeKey, maxDay] = localParseScenario(scenario)
|
function [mergeKey, maxDay] = localParseScenario(scenario)
|
||||||
%LOCALPARSESCENARIO Split "<mergeKey>_full" / "<mergeKey>_d0_10" apart.
|
%LOCALPARSESCENARIO Split "<mergeKey>_full" / "_d0_10" / "_d0_13" apart.
|
||||||
scenario = char(scenario);
|
scenario = char(scenario);
|
||||||
if endsWith(scenario, '_full')
|
if endsWith(scenario, '_full')
|
||||||
maxDay = 26;
|
maxDay = 26;
|
||||||
@@ -100,9 +100,12 @@ if endsWith(scenario, '_full')
|
|||||||
elseif endsWith(scenario, '_d0_10')
|
elseif endsWith(scenario, '_d0_10')
|
||||||
maxDay = 10;
|
maxDay = 10;
|
||||||
mergeKey = scenario(1:end - numel('_d0_10'));
|
mergeKey = scenario(1:end - numel('_d0_10'));
|
||||||
|
elseif endsWith(scenario, '_d0_13')
|
||||||
|
maxDay = 13; % day 13 = last day Box-A2 has data (fair A2-vs-B2 window)
|
||||||
|
mergeKey = scenario(1:end - numel('_d0_13'));
|
||||||
else
|
else
|
||||||
error('tdcs_scenario_data:badScenario', ...
|
error('tdcs_scenario_data:badScenario', ...
|
||||||
'Scenario "%s" does not end in "_full" or "_d0_10".', scenario);
|
'Scenario "%s" does not end in "_full", "_d0_10", or "_d0_13".', scenario);
|
||||||
end
|
end
|
||||||
|
|
||||||
if ~ismember(mergeKey, {'unmerged', 'mergeA2', 'mergeB2'})
|
if ~ismember(mergeKey, {'unmerged', 'mergeA2', 'mergeB2'})
|
||||||
|
|||||||
@@ -59,5 +59,21 @@ classdef tLme < matlab.unittest.TestCase
|
|||||||
testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario');
|
testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario');
|
||||||
end
|
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
|
||||||
end
|
end
|
||||||
|
|||||||
Reference in New Issue
Block a user