From 27d388530a3a4525d6a57dc1d954ce0fc1b1a116 Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Mon, 20 Jul 2026 07:51:30 -0400 Subject: [PATCH] matlab(tdcs): GLMM + per-animal + classification models --- analysis/matlab/tdcs_models.m | 216 ++++++++++++++++++++++++++++++++ analysis/matlab/tests/tModels.m | 121 ++++++++++++++++++ 2 files changed, 337 insertions(+) create mode 100644 analysis/matlab/tdcs_models.m create mode 100644 analysis/matlab/tests/tModels.m diff --git a/analysis/matlab/tdcs_models.m b/analysis/matlab/tdcs_models.m new file mode 100644 index 0000000..11f00b8 --- /dev/null +++ b/analysis/matlab/tdcs_models.m @@ -0,0 +1,216 @@ +function R = tdcs_models(S, meta, cfg) +%TDCS_MODELS Fit the tDCS GLMMs, per-animal tests, and (unmerged) classification. +% R = TDCS_MODELS(S, META, CFG) consumes one scenario table S and its +% META (from TDCS_SCENARIO_DATA) plus CFG (from TDCS_CONFIG), and +% returns a struct R: +% +% R.countGLME GeneralizedLinearMixedModel, Poisson count-level model: +% success ~ group + day_c + day_c2 + (1|subject) +% R.rateGLME GeneralizedLinearMixedModel, Binomial rate-level model +% (total==0 sessions dropped first): +% success ~ group + day_c + day_c2 + (1|subject) +% 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.perAnimal per-subject Box-B2 vs Box-A2 comparison (pure stats, +% no GLME involved): .countWelchP, .countMWUp, +% .rateWelchP, .rateMWUp, .nB, .nA. "count" = per-subject +% mean success; "rate" = per-subject pooled +% sum(success)/sum(total). Welch = TTEST2(..., +% 'Vartype','unequal') (two-sided); MWU = RANKSUM(..., +% 'tail','right') (one-sided, B2 > A2). +% R.h2 anchor check (Box-A2 vs Box-B2) read off the +% count/rate GLMEs' group_ fixed effect +% (coded vs the cfg.ref = anchorHigh reference level): +% .countIRR, .countOneSidedP, .rateOR, .rateOneSidedP. +% IRR/OR = 1/exp(coef) is the Box-B2-over-Box-A2 +% advantage (so > 1 supports H2); one-sided p = +% NORMCDF(coef/se), testing coef < 0 (B2 above A2). +% R.classify (unmerged scenarios only, i.e. meta.isUnmerged) struct +% with one field per unknown in cfg.unknowns (field name +% via matlab.lang.makeValidName), each holding .count and +% .rate sub-structs with .vsA2 (.ratio, .p), .vsB2 +% (.ratio, .p), and .label (one of 'B2-like', 'A2-like', +% 'AMBIGUOUS (nearer )', 'UNLIKE BOTH'), per +% METHODS.md's classification logic. For merged +% scenarios (meta.isUnmerged == false), R.classify is an +% empty (no-field) struct. +% +% See analysis/tdcs_glm.py (the Python this ports) and +% docs/superpowers/plans/2026-07-20-matlab-tdcs-analysis.md (Task 3). + +countFormula = 'success ~ group + day_c + day_c2 + (1|subject)'; +R.countGLME = fitglme(S, countFormula, 'Distribution', 'Poisson', 'Link', 'log'); + +Sr = S(S.total > 0, :); +R.rateGLME = fitglme(Sr, countFormula, 'Distribution', 'Binomial', ... + 'BinomialSize', Sr.total, 'Link', 'logit'); + +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); + +if meta.isUnmerged + R.classify = localClassifyAll(R.countGLME, R.rateGLME, cfg); +else + R.classify = struct(); +end + +end + +% --------------------------------------------------------------- per-animal +function out = localPerAnimal(S, cfg) +%LOCALPERANIMAL Per-subject Box-B2 vs Box-A2 Welch t-test / MWU (pure stats). +b = S(S.group == cfg.anchorHigh, :); +a = S(S.group == cfg.anchorLow, :); + +[bCount, bRate, nB] = localSubjectSummary(b); +[aCount, aRate, nA] = localSubjectSummary(a); + +[~, countWelchP] = ttest2(bCount, aCount, 'Vartype', 'unequal'); +countMWUp = ranksum(bCount, aCount, 'tail', 'right'); +[~, rateWelchP] = ttest2(bRate, aRate, 'Vartype', 'unequal'); +rateMWUp = ranksum(bRate, aRate, 'tail', 'right'); + +out = struct( ... + 'countWelchP', countWelchP, 'countMWUp', countMWUp, ... + 'rateWelchP', rateWelchP, 'rateMWUp', rateMWUp, ... + 'nB', nB, 'nA', nA); +end + +function [meanSuccess, pooledRate, nSubj] = localSubjectSummary(T) +%LOCALSUBJECTSUMMARY Per-subject mean success and pooled success/total. +subj = unique(T.subject); +nSubj = numel(subj); +meanSuccess = zeros(nSubj, 1); +pooledRate = zeros(nSubj, 1); +for i = 1:nSubj + rows = T.subject == subj(i); + meanSuccess(i) = mean(T.success(rows)); + pooledRate(i) = sum(T.success(rows)) / sum(T.total(rows)); +end +end + +% --------------------------------------------------------------------- H2 +function h2 = localH2(countGLME, rateGLME, cfg) +%LOCALH2 Anchor check: Box-A2 vs Box-B2, read off each GLME's group term. +[countCoef, countSE] = localCoefAndSE(countGLME, cfg.anchorLow, cfg.ref); +[rateCoef, rateSE] = localCoefAndSE(rateGLME, cfg.anchorLow, cfg.ref); + +h2 = struct( ... + 'countIRR', 1 / exp(countCoef), ... + 'countOneSidedP', normcdf(countCoef / countSE), ... + 'rateOR', 1 / exp(rateCoef), ... + 'rateOneSidedP', normcdf(rateCoef / rateSE)); +end + +% -------------------------------------------------------------- classify +function out = localClassifyAll(countGLME, rateGLME, cfg) +%LOCALCLASSIFYALL Classify every unknown in cfg.unknowns vs both anchors. +out = struct(); +for i = 1:numel(cfg.unknowns) + unknown = cfg.unknowns{i}; + fieldName = matlab.lang.makeValidName(unknown); + out.(fieldName) = struct( ... + 'count', localClassifyOne(countGLME, unknown, cfg), ... + 'rate', localClassifyOne(rateGLME, unknown, cfg)); +end +end + +function verdict = localClassifyOne(model, unknown, cfg) +%LOCALCLASSIFYONE Compare UNKNOWN vs both anchors on one GLME; label it. +% Mirrors classify_one() in analysis/tdcs_glm.py / METHODS.md's +% "Classification logic (unknowns)": indistinguishable (p >= 0.05) from +% one anchor and different from the other => that anchor's label; +% indistinguishable from both => AMBIGUOUS (report the numerically +% nearer one); different from both => UNLIKE BOTH. +[coefLo, pLo] = localContrast(model, unknown, cfg.anchorLow, cfg.ref); +[coefHi, pHi] = localContrast(model, unknown, cfg.anchorHigh, cfg.ref); + +vsA2 = struct('ratio', exp(coefLo), 'p', pLo); +vsB2 = struct('ratio', exp(coefHi), 'p', pHi); + +likeLo = pLo >= 0.05; +likeHi = pHi >= 0.05; +if likeHi && ~likeLo + label = 'B2-like'; +elseif likeLo && ~likeHi + label = 'A2-like'; +elseif likeLo && likeHi + if abs(coefHi) < abs(coefLo) + nearer = cfg.anchorHigh; + else + nearer = cfg.anchorLow; + end + label = sprintf('AMBIGUOUS (nearer %s)', nearer); +else + label = 'UNLIKE BOTH'; +end + +verdict = struct('vsA2', vsA2, 'vsB2', vsB2, 'label', label); +end + +% ------------------------------------------------------------- contrasts +function [coef, se] = localCoefAndSE(model, aName, refName) +%LOCALCOEFANDSE Coefficient + SE for group level A_NAME vs the reference. +v = localContrastVector(model, aName, refName, refName); +beta = model.Coefficients.Estimate; +Cov = model.CoefficientCovariance; +coef = v * beta; +se = sqrt(v * Cov * v'); +end + +function [coef, p] = localContrast(model, aName, bName, refName) +%LOCALCONTRAST Linear contrast A_NAME vs B_NAME (both coded vs REF_NAME). +% Two-sided p via COEFTEST (an F-test on one row = the two-sided +% Wald/t-test); works for any pair, including B_NAME == REF_NAME. +v = localContrastVector(model, aName, bName, refName); +beta = model.Coefficients.Estimate; +coef = v * beta; +p = coefTest(model, v); +end + +function v = localContrastVector(model, aName, bName, refName) +%LOCALCONTRASTVECTOR +1 at group_A_NAME (if not the reference), -1 at +% group_B_NAME (if not the reference), over MODEL.CoefficientNames. +cn = model.CoefficientNames; +v = zeros(1, numel(cn)); +if ~strcmp(aName, refName) + v(localCoefIndex(cn, aName)) = v(localCoefIndex(cn, aName)) + 1; +end +if ~strcmp(bName, refName) + v(localCoefIndex(cn, bName)) = v(localCoefIndex(cn, bName)) - 1; +end +end + +function idx = localCoefIndex(coefNames, groupName) +%LOCALCOEFINDEX Index into CoefficientNames for a "group_" term. +idx = find(strcmp(coefNames, ['group_' groupName]), 1); +if isempty(idx) + error('tdcs_models:missingGroupTerm', ... + 'No "group_%s" coefficient in this GLME (CoefficientNames: %s).', ... + groupName, strjoin(coefNames, ', ')); +end +end + +% -------------------------------------------------------------- learning +function out = localJointInteractionTest(learnGLME) +%LOCALJOINTINTERACTIONTEST Joint Wald test that all group x day_c +% interaction fixed effects are zero (do groups learn at different +% rates?), mirroring the GEE wald_test in analysis/tdcs_glm.py. +cn = learnGLME.CoefficientNames; +isInteraction = contains(cn, ':day_c') & ~contains(cn, 'day_c2'); +idx = find(isInteraction); +H = zeros(numel(idx), numel(cn)); +for i = 1:numel(idx) + H(i, idx(i)) = 1; +end +[p, F, df1, df2] = coefTest(learnGLME, H); +out = struct('jointP', p, 'jointF', F, 'jointDF1', df1, 'jointDF2', df2); +end diff --git a/analysis/matlab/tests/tModels.m b/analysis/matlab/tests/tModels.m new file mode 100644 index 0000000..f421107 --- /dev/null +++ b/analysis/matlab/tests/tModels.m @@ -0,0 +1,121 @@ +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 + end +end