diff --git a/analysis/matlab/run_all.m b/analysis/matlab/run_all.m index 9cd6d11..88f5873 100644 --- a/analysis/matlab/run_all.m +++ b/analysis/matlab/run_all.m @@ -7,7 +7,8 @@ scenarios = {'unmerged_full', 'unmerged_d0_10', ... 'mergeB2_full', 'mergeB2_d0_10', ... 'lme_unmerged_full', 'lme_unmerged_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) fprintf('\n\n### Running scenario: %s ###\n', scenarios{i}); diff --git a/analysis/matlab/tdcs_glm.m b/analysis/matlab/tdcs_glm.m index 762366d..65e0ca6 100644 --- a/analysis/matlab/tdcs_glm.m +++ b/analysis/matlab/tdcs_glm.m @@ -10,9 +10,11 @@ function tdcs_glm(scenario) % mergeA2_full, mergeA2_d0_10, % mergeB2_full, mergeB2_d0_10 % or a linear "days x tDCS" mixed model (Box-B2 vs Box-A2, per the paper): -% lme_unmerged_full, lme_unmerged_d0_10, -% lme_mergeA2_full, lme_mergeA2_d0_10, -% lme_mergeB2_full, lme_mergeB2_d0_10 +% lme_unmerged_full, lme_unmerged_d0_10, lme_unmerged_d0_13, +% lme_mergeA2_full, lme_mergeA2_d0_10, lme_mergeA2_d0_13, +% 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 % 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); 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 error('tdcs_glm:badScenario', ... ['Unrecognized scenario "%s". Valid scenarios are: ' ... diff --git a/analysis/matlab/tdcs_lme.m b/analysis/matlab/tdcs_lme.m index daef865..923ac4b 100644 --- a/analysis/matlab/tdcs_lme.m +++ b/analysis/matlab/tdcs_lme.m @@ -36,6 +36,8 @@ A = anova(lme); L.lme = lme; L.nSubjects = numel(unique(T.subject)); 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.day = localTerm(C, A, 'day'); L.tDCS = localTerm(C, A, 'tDCS'); diff --git a/analysis/matlab/tdcs_lme_report.m b/analysis/matlab/tdcs_lme_report.m index 14401aa..1a0e993 100644 --- a/analysis/matlab/tdcs_lme_report.m +++ b/analysis/matlab/tdcs_lme_report.m @@ -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', ... cfg.anchorHigh, cfg.anchorLow)]; 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(' 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('%s\n', repmat('-', 1, 70))]; diff --git a/analysis/matlab/tdcs_scenario_data.m b/analysis/matlab/tdcs_scenario_data.m index 4e605ad..6daea18 100644 --- a/analysis/matlab/tdcs_scenario_data.m +++ b/analysis/matlab/tdcs_scenario_data.m @@ -92,7 +92,7 @@ meta.isUnmerged = strcmp(mergeKey, 'unmerged'); end function [mergeKey, maxDay] = localParseScenario(scenario) -%LOCALPARSESCENARIO Split "_full" / "_d0_10" apart. +%LOCALPARSESCENARIO Split "_full" / "_d0_10" / "_d0_13" apart. scenario = char(scenario); if endsWith(scenario, '_full') maxDay = 26; @@ -100,9 +100,12 @@ if endsWith(scenario, '_full') elseif endsWith(scenario, '_d0_10') maxDay = 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 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 if ~ismember(mergeKey, {'unmerged', 'mergeA2', 'mergeB2'}) diff --git a/analysis/matlab/tests/tLme.m b/analysis/matlab/tests/tLme.m index 686fc1f..16928d0 100644 --- a/analysis/matlab/tests/tLme.m +++ b/analysis/matlab/tests/tLme.m @@ -59,5 +59,21 @@ classdef tLme < matlab.unittest.TestCase testCase.verifyError(@() tdcs_glm('lme_nonsense'), 'tdcs_glm:badScenario'); 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