From b1057f61056c97fa3456f7271029d42a1a8ef4a2 Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Mon, 20 Jul 2026 08:40:47 -0400 Subject: [PATCH] fix(matlab): report learning rate via per-animal slope test, not anticonservative GLMM F-test fitglme coefTest only offers observation-level DF (df2~=178), which overstates the group x day interaction significance for this few-subject design (p~=0 in every scenario). Replace the reported learning-rate result with a cluster-honest per-animal slope test (per-subject OLS slope of success vs day, Welch + MWU, B2 vs A2) -- now correctly parallel in all 6 scenarios, matching the Python GEE. The GLMM F-test is kept as a flagged reference only. Adds a regression test. Co-Authored-By: Claude Opus 4.8 --- analysis/matlab/tdcs_models.m | 45 ++++++++++++++++++++++++++++++--- analysis/matlab/tdcs_report.m | 18 +++++++++---- analysis/matlab/tests/tModels.m | 16 ++++++++++++ 3 files changed, 70 insertions(+), 9 deletions(-) diff --git a/analysis/matlab/tdcs_models.m b/analysis/matlab/tdcs_models.m index 11f00b8..c7b17eb 100644 --- a/analysis/matlab/tdcs_models.m +++ b/analysis/matlab/tdcs_models.m @@ -12,9 +12,18 @@ function R = tdcs_models(S, meta, cfg) % R.learnGLME GeneralizedLinearMixedModel, Poisson learning-rate % model with a group x day_c interaction: % success ~ group*day_c + day_c2 + (1|subject) -% R.learn struct with the joint Wald test (via COEFTEST) that -% all group x day_c interaction terms are zero: -% .jointP, .jointF, .jointDF1, .jointDF2. +% R.learn learning-rate comparison, Box-B2 vs Box-A2. The reported +% result is a PER-ANIMAL slope test (cluster-honest, like +% R.perAnimal): each subject's OLS slope of success on day, +% then Welch t / two-sided Mann-Whitney across groups: +% .slopeWelchP, .slopeMWUp, .nB, .nA, .meanSlopeB, +% .meanSlopeA. The GLMM group x day_c joint F-test is also +% kept for reference (.glmeJointP, .glmeJointF, .glmeJointDF1, +% .glmeJointDF2) but is ANTICONSERVATIVE here: fitglme's +% coefTest offers only observation-level DF, which overstates +% significance for this few-subject design (it reads p~=0 +% even when the per-animal slopes are clearly parallel), so +% it is NOT the basis for the conclusion. % R.perAnimal per-subject Box-B2 vs Box-A2 comparison (pure stats, % no GLME involved): .countWelchP, .countMWUp, % .rateWelchP, .rateMWUp, .nB, .nA. "count" = per-subject @@ -51,10 +60,10 @@ R.rateGLME = fitglme(Sr, countFormula, 'Distribution', 'Binomial', ... learnFormula = 'success ~ group*day_c + day_c2 + (1|subject)'; R.learnGLME = fitglme(S, learnFormula, 'Distribution', 'Poisson'); -R.learn = localJointInteractionTest(R.learnGLME); R.perAnimal = localPerAnimal(S, cfg); R.h2 = localH2(R.countGLME, R.rateGLME, cfg); +R.learn = localLearn(R.learnGLME, S, cfg); if meta.isUnmerged R.classify = localClassifyAll(R.countGLME, R.rateGLME, cfg); @@ -200,6 +209,34 @@ end end % -------------------------------------------------------------- learning +function out = localLearn(learnGLME, S, cfg) +%LOCALLEARN Learning-rate result: per-animal slope test (primary, cluster- +% honest) plus the GLMM interaction joint F-test (reference, anticonservative). +bSlopes = localSubjectSlopes(S(S.group == cfg.anchorHigh, :)); +aSlopes = localSubjectSlopes(S(S.group == cfg.anchorLow, :)); +[~, slopeWelchP] = ttest2(bSlopes, aSlopes, 'Vartype', 'unequal'); +slopeMWUp = ranksum(bSlopes, aSlopes); % two-sided: slope-diff direction unknown + +j = localJointInteractionTest(learnGLME); +out = struct( ... + 'slopeWelchP', slopeWelchP, 'slopeMWUp', slopeMWUp, ... + 'nB', numel(bSlopes), 'nA', numel(aSlopes), ... + 'meanSlopeB', mean(bSlopes), 'meanSlopeA', mean(aSlopes), ... + 'glmeJointP', j.jointP, 'glmeJointF', j.jointF, ... + 'glmeJointDF1', j.jointDF1, 'glmeJointDF2', j.jointDF2); +end + +function slopes = localSubjectSlopes(T) +%LOCALSUBJECTSLOPES Per-subject OLS slope of success on day. +subj = unique(T.subject); +slopes = zeros(numel(subj), 1); +for i = 1:numel(subj) + rows = T.subject == subj(i); + coeffs = polyfit(T.day(rows), T.success(rows), 1); + slopes(i) = coeffs(1); +end +end + function out = localJointInteractionTest(learnGLME) %LOCALJOINTINTERACTIONTEST Joint Wald test that all group x day_c % interaction fixed effects are zero (do groups learn at different diff --git a/analysis/matlab/tdcs_report.m b/analysis/matlab/tdcs_report.m index 124523b..066864d 100644 --- a/analysis/matlab/tdcs_report.m +++ b/analysis/matlab/tdcs_report.m @@ -87,11 +87,19 @@ s = [s localCleanDisp(R.countGLME) sprintf('\n')]; s = [s localHeader('(B) LEVEL / RATE -- Binomial GLMM (subject random intercept)')]; s = [s localCleanDisp(R.rateGLME) sprintf('\n')]; -s = [s localHeader('(C) LEARNING RATE -- joint test (all group x day_c interactions = 0)')]; -s = [s sprintf('F=%.3f, df1=%d, df2=%.2f, p=%.4f\n', ... - R.learn.jointF, R.learn.jointDF1, R.learn.jointDF2, R.learn.jointP)]; -s = [s sprintf([' -> small p = groups improve at DIFFERENT rates; ' ... - 'large p = parallel learning.\n\n'])]; +s = [s localHeader('(C) LEARNING RATE -- per-animal slope test (Box-B2 vs Box-A2)')]; +if R.learn.slopeMWUp < 0.05 + verdict = 'DIFFERENT learning rates'; +else + verdict = 'parallel learning (no slope difference detected)'; +end +s = [s sprintf(['Per-subject OLS slope of success vs day: Box-B2 mean=%.2f, Box-A2 mean=%.2f\n' ... + ' Welch two-sided p=%.3f, Mann-Whitney p=%.3f (nB=%d, nA=%d) -> %s\n'], ... + R.learn.meanSlopeB, R.learn.meanSlopeA, R.learn.slopeWelchP, R.learn.slopeMWUp, ... + R.learn.nB, R.learn.nA, verdict)]; +s = [s sprintf(['(Reference only: the GLMM group x day_c joint F-test gives p=%.4f, but with just\n' ... + ' observation-level DF (df2=%.0f) it is ANTICONSERVATIVE for this few-subject design and is\n' ... + ' NOT the basis for the conclusion above.)\n\n'], R.learn.glmeJointP, R.learn.glmeJointDF2)]; end function s = localCleanDisp(model) diff --git a/analysis/matlab/tests/tModels.m b/analysis/matlab/tests/tModels.m index f421107..46842e6 100644 --- a/analysis/matlab/tests/tModels.m +++ b/analysis/matlab/tests/tModels.m @@ -117,5 +117,21 @@ classdef tModels < matlab.unittest.TestCase testCase.verifyFalse(metam.isUnmerged); testCase.verifyEmpty(fieldnames(Rm.classify)); end + + function testLearningPerAnimalSlopesParallel(testCase) + % The reported learning-rate result is a per-animal slope test and + % must show parallel learning (non-significant) -- matching the + % Python's cluster-robust conclusion. Guards against regressing to + % the anticonservative GLMM interaction F-test (which reads p~=0). + cfg = tdcs_config(); + for sc = {'unmerged_full', 'mergeA2_full', 'mergeB2_full'} + [S, meta] = tdcs_scenario_data(sc{1}); + R = tdcs_models(S, meta, cfg); + testCase.verifyGreaterThan(R.learn.slopeMWUp, 0.05, ... + sprintf('%s: per-animal learning slopes should be parallel (MWU p>0.05)', sc{1})); + testCase.verifyGreaterThan(R.learn.slopeWelchP, 0.05, ... + sprintf('%s: per-animal learning slopes should be parallel (Welch p>0.05)', sc{1})); + end + end end end