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 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
@@ -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
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user