Experiments DB Dev
b1057f6105
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>
138 lines
6.8 KiB
Matlab
138 lines
6.8 KiB
Matlab
classdef tModels < matlab.unittest.TestCase
|
|
%TMODELS Tests for tdcs_models (Task 3: GLMM + per-animal + classify).
|
|
%
|
|
% Authored per TDD: written before tdcs_models.m existed, so it is
|
|
% expected to fail until the model-fitting function is implemented.
|
|
% MATLAB R2025b + Statistics and Machine Learning Toolbox IS licensed
|
|
% in this environment (unlike Tasks 1-2), so this suite is actually
|
|
% run: RED before tdcs_models.m existed, GREEN after.
|
|
%
|
|
% Per-animal golden values (Box-B2 vs Box-A2, pure stats: per-subject
|
|
% mean success for "count", per-subject pooled sum(success)/sum(total)
|
|
% for "rate"; Welch two-sided TTEST2, one-sided 'right' RANKSUM) were
|
|
% independently reproduced with a throwaway Python script against
|
|
% analysis/tdcs_reach_data.csv (equal_var=False t-test,
|
|
% mannwhitneyu(..., alternative='greater')), confirming the plan's
|
|
% Task 3 golden-value table to better than 1e-3 in every cell (see
|
|
% task-3-report.md). They must match to 3 decimals (AbsTol 0.005) --
|
|
% do not loosen this tolerance.
|
|
%
|
|
% GLME-derived checks (H2 odds/incidence-rate ratios, classification)
|
|
% are checked loosely (direction/magnitude), since GLMM (subject
|
|
% random intercept) is not numerically identical to the Python's GEE
|
|
% (population-average, robust SEs) -- see the plan's Global
|
|
% Constraints and analysis/METHODS.md.
|
|
|
|
properties (TestParameter)
|
|
perAnimalCase = struct( ...
|
|
'unmerged_full', struct( ...
|
|
'scenario', 'unmerged_full', 'nB', 3, 'nA', 3, ...
|
|
'rateWelchP', 0.072, 'rateMWUp', 0.050, ...
|
|
'countWelchP', 0.055, 'countMWUp', 0.050), ...
|
|
'mergeA2_full', struct( ...
|
|
'scenario', 'mergeA2_full', 'nB', 4, 'nA', 4, ...
|
|
'rateWelchP', 0.056, 'rateMWUp', 0.029, ...
|
|
'countWelchP', 0.035, 'countMWUp', 0.029), ...
|
|
'mergeB2_full', struct( ...
|
|
'scenario', 'mergeB2_full', 'nB', 5, 'nA', 3, ...
|
|
'rateWelchP', 0.044, 'rateMWUp', 0.018, ...
|
|
'countWelchP', 0.020, 'countMWUp', 0.018), ...
|
|
'unmerged_d0_10', struct( ...
|
|
'scenario', 'unmerged_d0_10', 'nB', 3, 'nA', 3, ...
|
|
'rateWelchP', 0.070, 'rateMWUp', 0.050, ...
|
|
'countWelchP', 0.041, 'countMWUp', 0.050), ...
|
|
'mergeA2_d0_10', struct( ...
|
|
'scenario', 'mergeA2_d0_10', 'nB', 4, 'nA', 4, ...
|
|
'rateWelchP', 0.069, 'rateMWUp', 0.057, ...
|
|
'countWelchP', 0.023, 'countMWUp', 0.014), ...
|
|
'mergeB2_d0_10', struct( ...
|
|
'scenario', 'mergeB2_d0_10', 'nB', 5, 'nA', 3, ...
|
|
'rateWelchP', 0.094, 'rateMWUp', 0.071, ...
|
|
'countWelchP', 0.023, 'countMWUp', 0.018) ...
|
|
)
|
|
end
|
|
|
|
methods (Test)
|
|
function testPerAnimalGoldenValues(testCase, perAnimalCase)
|
|
cfg = tdcs_config();
|
|
[S, meta] = tdcs_scenario_data(perAnimalCase.scenario);
|
|
R = tdcs_models(S, meta, cfg);
|
|
|
|
testCase.verifyEqual(R.perAnimal.nB, perAnimalCase.nB);
|
|
testCase.verifyEqual(R.perAnimal.nA, perAnimalCase.nA);
|
|
testCase.verifyEqual(R.perAnimal.rateWelchP, ...
|
|
perAnimalCase.rateWelchP, 'AbsTol', 0.005);
|
|
testCase.verifyEqual(R.perAnimal.rateMWUp, ...
|
|
perAnimalCase.rateMWUp, 'AbsTol', 0.005);
|
|
testCase.verifyEqual(R.perAnimal.countWelchP, ...
|
|
perAnimalCase.countWelchP, 'AbsTol', 0.005);
|
|
testCase.verifyEqual(R.perAnimal.countMWUp, ...
|
|
perAnimalCase.countMWUp, 'AbsTol', 0.005);
|
|
end
|
|
|
|
function testGLMEObjectsReturned(testCase)
|
|
cfg = tdcs_config();
|
|
[S, meta] = tdcs_scenario_data('mergeB2_full');
|
|
R = tdcs_models(S, meta, cfg);
|
|
|
|
testCase.verifyClass(R.countGLME, 'GeneralizedLinearMixedModel');
|
|
testCase.verifyClass(R.rateGLME, 'GeneralizedLinearMixedModel');
|
|
testCase.verifyClass(R.learnGLME, 'GeneralizedLinearMixedModel');
|
|
end
|
|
|
|
function testH2SanityMergeB2Full(testCase)
|
|
cfg = tdcs_config();
|
|
[S, meta] = tdcs_scenario_data('mergeB2_full');
|
|
R = tdcs_models(S, meta, cfg);
|
|
|
|
testCase.verifyGreaterThan(R.h2.rateOR, 1.3);
|
|
testCase.verifyLessThan(R.h2.rateOR, 1.7);
|
|
testCase.verifyLessThan(R.h2.rateOneSidedP, 0.05);
|
|
testCase.verifyGreaterThan(R.h2.countIRR, 1);
|
|
end
|
|
|
|
function testClassifyPresentForUnmergedOnly(testCase)
|
|
cfg = tdcs_config();
|
|
|
|
[S, meta] = tdcs_scenario_data('unmerged_full');
|
|
R = tdcs_models(S, meta, cfg);
|
|
testCase.verifyTrue(meta.isUnmerged);
|
|
for i = 1:numel(cfg.unknowns)
|
|
fieldName = matlab.lang.makeValidName(cfg.unknowns{i});
|
|
testCase.verifyTrue(isfield(R.classify, fieldName), ...
|
|
sprintf('R.classify missing entry for unknown "%s".', ...
|
|
cfg.unknowns{i}));
|
|
entry = R.classify.(fieldName);
|
|
testCase.verifyTrue(isfield(entry, 'count'));
|
|
testCase.verifyTrue(isfield(entry, 'rate'));
|
|
testCase.verifyTrue(isfield(entry.count, 'label'));
|
|
testCase.verifyTrue(isfield(entry.count.vsA2, 'ratio'));
|
|
testCase.verifyTrue(isfield(entry.count.vsA2, 'p'));
|
|
testCase.verifyTrue(isfield(entry.count.vsB2, 'ratio'));
|
|
testCase.verifyTrue(isfield(entry.count.vsB2, 'p'));
|
|
end
|
|
|
|
[Sm, metam] = tdcs_scenario_data('mergeB2_full');
|
|
Rm = tdcs_models(Sm, metam, cfg);
|
|
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
|