analysis(matlab): Box-A pooling variations (b2 vs a2+a, b2+a vs a2)

Add make_boxa_variations.m producing 6 variation folders (data.csv + analyze.m
+ result.txt) for the paper LME over windows 0-5, 6-10, 0-10:
  boxa_a2  B2         vs A2 + Box-A   (b2 vs a2+a)
  boxa_b2  B2 + Box-A vs A2           (b2+a vs a2)
plus variations/boxa_summary.csv (residual / Satterthwaite / random-slope
interaction p). Early window (0-5) is obs-level significant (res p~0.03-0.04)
but n.s. under the honest random-slope test (rs p~0.13); later windows n.s.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
2026-07-23 13:39:01 -04:00
parent b03549982e
commit e1af1af7f4
20 changed files with 1370 additions and 0 deletions
@@ -0,0 +1,98 @@
% Variation analysis -- the paper's linear mixed model on the successful-reach
% COUNT, fit on this folder's curated data subset.
%
% model: behavior ~ stim + day + stim:day + (1|rat)
% behavior = successful reaches (count per session)
% stim = 1 for the treatment group(s), 0 for the control group(s)
% day = training day within this window (0 = first analyzed day)
% rat = subject (random intercept)
%
% Self-contained: reads data.csv beside this script and writes result.txt.
% Run headless from this folder with: matlab -batch "analyze"
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m);
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
% Honest test: refit with a per-animal random SLOPE so the interaction DF
% collapses toward the animal count (guarded -- may not converge in short windows).
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all');
try
mr = fitlme(tbl, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite');
ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
warning(wst);
gi = @(t) find(strcmp(C.Name, t), 1);
ga = @(t) find(strcmp(A.Term, t), 1);
gs = @(t) find(strcmp(As.Term, t), 1);
row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f)\n', nm, ...
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
As.pValue(gs(t)), As.DF2(gs(t)));
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2
cov = '** WARNING: unequal day coverage -- interaction may be confounded. **';
else
cov = '(equal day coverage over this window)';
end
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
if pI >= 0.05
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
elseif eI > 0
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
else
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
end
bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
s = [s sprintf('N = %d rats, %d sessions raw day coverage: treat %d..%d, control %d..%d\n', ...
numel(unique(D.subject)), height(D), minT, maxT, minC, maxC)];
s = [s sprintf('%s\n\n%s\nFULL MODEL SUMMARY -- fitlme\n%s\n%s\n', cov, bar, bar, raw)];
s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (resid)', 'Satterthwaite: p (df)', repmat('-', 1, 76))];
s = [s row('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')];
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): did not converge for this window.\n')];
end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
fprintf('%s', s);
fid = fopen(fullfile(here, 'result.txt'), 'w');
fprintf(fid, '%s', s);
fclose(fid);
% Machine-readable handoff for the summary table (see make_variations.m).
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
'covEqual', abs(maxT - maxC) <= 2);
@@ -0,0 +1,73 @@
subject,group,day,success,total,stim
Banh-mi-1,Electrode-Box-B2,0,13,84,1
Banh-mi-1,Electrode-Box-B2,1,35,86,1
Banh-mi-1,Electrode-Box-B2,2,47,110,1
Banh-mi-1,Electrode-Box-B2,3,65,140,1
Banh-mi-1,Electrode-Box-B2,4,70,127,1
Banh-mi-1,Electrode-Box-B2,5,102,142,1
Banh-mi-1,Electrode-Box-B2,6,90,131,1
Banh-mi-1,Electrode-Box-B2,7,109,148,1
Banh-mi-1,Electrode-Box-B2,8,104,137,1
Banh-mi-1,Electrode-Box-B2,9,119,150,1
Banh-mi-1,Electrode-Box-B2,10,121,158,1
Banh-mi-2,Electrode-Box-A2,0,25,97,0
Banh-mi-2,Electrode-Box-A2,1,22,101,0
Banh-mi-2,Electrode-Box-A2,2,22,119,0
Banh-mi-2,Electrode-Box-A2,3,29,118,0
Banh-mi-2,Electrode-Box-A2,4,27,136,0
Banh-mi-2,Electrode-Box-A2,5,43,146,0
Banh-mi-2,Electrode-Box-A2,6,74,146,0
Banh-mi-2,Electrode-Box-A2,7,70,148,0
Banh-mi-2,Electrode-Box-A2,8,65,130,0
Banh-mi-2,Electrode-Box-A2,9,79,151,0
Banh-mi-2,Electrode-Box-A2,10,93,152,0
Egg-tart-1,Electrode-Box-B2,0,7,56,1
Egg-tart-1,Electrode-Box-B2,1,16,78,1
Egg-tart-1,Electrode-Box-B2,2,23,103,1
Egg-tart-1,Electrode-Box-B2,3,63,120,1
Egg-tart-1,Electrode-Box-B2,4,69,132,1
Egg-tart-1,Electrode-Box-B2,5,83,136,1
Egg-tart-1,Electrode-Box-B2,6,71,142,1
Egg-tart-1,Electrode-Box-B2,7,79,138,1
Egg-tart-1,Electrode-Box-B2,8,98,142,1
Egg-tart-1,Electrode-Box-B2,9,89,139,1
Egg-tart-1,Electrode-Box-B2,10,96,143,1
Egg-tart-2,Electrode-Box-A2,0,9,32,0
Egg-tart-2,Electrode-Box-A2,1,2,38,0
Egg-tart-2,Electrode-Box-A2,2,31,93,0
Egg-tart-2,Electrode-Box-A2,3,44,101,0
Egg-tart-2,Electrode-Box-A2,4,54,131,0
Egg-tart-2,Electrode-Box-A2,5,84,139,0
Egg-tart-2,Electrode-Box-A2,6,85,145,0
Egg-tart-2,Electrode-Box-A2,7,79,143,0
Egg-tart-2,Electrode-Box-A2,8,76,131,0
Egg-tart-2,Electrode-Box-A2,9,88,149,0
Egg-tart-2,Electrode-Box-A2,10,81,151,0
Khoai-tay-1,Electrode-Box-A,0,14,62,1
Khoai-tay-1,Electrode-Box-A,1,22,83,1
Khoai-tay-1,Electrode-Box-A,2,6,79,1
Khoai-tay-1,Electrode-Box-A,3,28,97,1
Khoai-tay-1,Electrode-Box-A,4,60,134,1
Khoai-tay-1,Electrode-Box-A,5,75,138,1
Khoai-tay-1,Electrode-Box-A,6,81,137,1
Khoai-tay-1,Electrode-Box-A,7,78,147,1
Khoai-tay-1,Electrode-Box-A,8,91,132,1
Khoai-tay-1,Electrode-Box-A,9,99,146,1
Khoai-tay-1,Electrode-Box-A,10,94,143,1
Root-beer-1,Electrode-Box-B2,0,11,85,1
Root-beer-1,Electrode-Box-B2,1,18,76,1
Root-beer-1,Electrode-Box-B2,2,40,105,1
Root-beer-1,Electrode-Box-B2,3,55,134,1
Root-beer-1,Electrode-Box-B2,4,75,136,1
Root-beer-1,Electrode-Box-B2,5,64,133,1
Root-beer-1,Electrode-Box-B2,6,104,139,1
Root-beer-1,Electrode-Box-B2,7,98,148,1
Root-beer-1,Electrode-Box-B2,8,81,145,1
Root-beer-1,Electrode-Box-B2,9,89,156,1
Root-beer-1,Electrode-Box-B2,10,105,158,1
Root-beer-2,Electrode-Box-A2,0,22,74,0
Root-beer-2,Electrode-Box-A2,1,31,87,0
Root-beer-2,Electrode-Box-A2,2,49,134,0
Root-beer-2,Electrode-Box-A2,3,31,89,0
Root-beer-2,Electrode-Box-A2,4,60,140,0
Root-beer-2,Electrode-Box-A2,5,84,147,0
1 subject group day success total stim
2 Banh-mi-1 Electrode-Box-B2 0 13 84 1
3 Banh-mi-1 Electrode-Box-B2 1 35 86 1
4 Banh-mi-1 Electrode-Box-B2 2 47 110 1
5 Banh-mi-1 Electrode-Box-B2 3 65 140 1
6 Banh-mi-1 Electrode-Box-B2 4 70 127 1
7 Banh-mi-1 Electrode-Box-B2 5 102 142 1
8 Banh-mi-1 Electrode-Box-B2 6 90 131 1
9 Banh-mi-1 Electrode-Box-B2 7 109 148 1
10 Banh-mi-1 Electrode-Box-B2 8 104 137 1
11 Banh-mi-1 Electrode-Box-B2 9 119 150 1
12 Banh-mi-1 Electrode-Box-B2 10 121 158 1
13 Banh-mi-2 Electrode-Box-A2 0 25 97 0
14 Banh-mi-2 Electrode-Box-A2 1 22 101 0
15 Banh-mi-2 Electrode-Box-A2 2 22 119 0
16 Banh-mi-2 Electrode-Box-A2 3 29 118 0
17 Banh-mi-2 Electrode-Box-A2 4 27 136 0
18 Banh-mi-2 Electrode-Box-A2 5 43 146 0
19 Banh-mi-2 Electrode-Box-A2 6 74 146 0
20 Banh-mi-2 Electrode-Box-A2 7 70 148 0
21 Banh-mi-2 Electrode-Box-A2 8 65 130 0
22 Banh-mi-2 Electrode-Box-A2 9 79 151 0
23 Banh-mi-2 Electrode-Box-A2 10 93 152 0
24 Egg-tart-1 Electrode-Box-B2 0 7 56 1
25 Egg-tart-1 Electrode-Box-B2 1 16 78 1
26 Egg-tart-1 Electrode-Box-B2 2 23 103 1
27 Egg-tart-1 Electrode-Box-B2 3 63 120 1
28 Egg-tart-1 Electrode-Box-B2 4 69 132 1
29 Egg-tart-1 Electrode-Box-B2 5 83 136 1
30 Egg-tart-1 Electrode-Box-B2 6 71 142 1
31 Egg-tart-1 Electrode-Box-B2 7 79 138 1
32 Egg-tart-1 Electrode-Box-B2 8 98 142 1
33 Egg-tart-1 Electrode-Box-B2 9 89 139 1
34 Egg-tart-1 Electrode-Box-B2 10 96 143 1
35 Egg-tart-2 Electrode-Box-A2 0 9 32 0
36 Egg-tart-2 Electrode-Box-A2 1 2 38 0
37 Egg-tart-2 Electrode-Box-A2 2 31 93 0
38 Egg-tart-2 Electrode-Box-A2 3 44 101 0
39 Egg-tart-2 Electrode-Box-A2 4 54 131 0
40 Egg-tart-2 Electrode-Box-A2 5 84 139 0
41 Egg-tart-2 Electrode-Box-A2 6 85 145 0
42 Egg-tart-2 Electrode-Box-A2 7 79 143 0
43 Egg-tart-2 Electrode-Box-A2 8 76 131 0
44 Egg-tart-2 Electrode-Box-A2 9 88 149 0
45 Egg-tart-2 Electrode-Box-A2 10 81 151 0
46 Khoai-tay-1 Electrode-Box-A 0 14 62 1
47 Khoai-tay-1 Electrode-Box-A 1 22 83 1
48 Khoai-tay-1 Electrode-Box-A 2 6 79 1
49 Khoai-tay-1 Electrode-Box-A 3 28 97 1
50 Khoai-tay-1 Electrode-Box-A 4 60 134 1
51 Khoai-tay-1 Electrode-Box-A 5 75 138 1
52 Khoai-tay-1 Electrode-Box-A 6 81 137 1
53 Khoai-tay-1 Electrode-Box-A 7 78 147 1
54 Khoai-tay-1 Electrode-Box-A 8 91 132 1
55 Khoai-tay-1 Electrode-Box-A 9 99 146 1
56 Khoai-tay-1 Electrode-Box-A 10 94 143 1
57 Root-beer-1 Electrode-Box-B2 0 11 85 1
58 Root-beer-1 Electrode-Box-B2 1 18 76 1
59 Root-beer-1 Electrode-Box-B2 2 40 105 1
60 Root-beer-1 Electrode-Box-B2 3 55 134 1
61 Root-beer-1 Electrode-Box-B2 4 75 136 1
62 Root-beer-1 Electrode-Box-B2 5 64 133 1
63 Root-beer-1 Electrode-Box-B2 6 104 139 1
64 Root-beer-1 Electrode-Box-B2 7 98 148 1
65 Root-beer-1 Electrode-Box-B2 8 81 145 1
66 Root-beer-1 Electrode-Box-B2 9 89 156 1
67 Root-beer-1 Electrode-Box-B2 10 105 158 1
68 Root-beer-2 Electrode-Box-A2 0 22 74 0
69 Root-beer-2 Electrode-Box-A2 1 31 87 0
70 Root-beer-2 Electrode-Box-A2 2 49 134 0
71 Root-beer-2 Electrode-Box-A2 3 31 89 0
72 Root-beer-2 Electrode-Box-A2 4 60 140 0
73 Root-beer-2 Electrode-Box-A2 5 84 147 0
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: boxa_b2_d0_10
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
control (stim=0): Electrode-Box-A2
N = 7 rats, 72 sessions raw day coverage: treat 0..10, control 0..10
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 72
Fixed effects coefficients 4
Random effects coefficients 7
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
587.74 601.4 -287.87 575.74
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 17.308 5.4187 3.1941 68 0.0021263
{'day' } 8.0028 0.78816 10.154 68 2.9164e-15
{'stim' } 1.8623 7.1264 0.26132 68 0.79463
{'day:stim' } 1.604 0.98592 1.6269 68 0.10838
Lower Upper
6.4953 28.121
6.4301 9.5756
-12.358 16.083
-0.36337 3.5714
Random effects covariance parameters (95% CIs):
Group: rat (7 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 6.0474
Lower Upper
2.8281 12.931
Group: Error
Name Estimate Lower Upper
{'Res Std'} 12.425 10.46 14.759
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(68)= 1.63 F(1)= 2.647 p=0.1084 p=0.1084 (df=68)
day (learning) t(68)= 10.15 F(1)=103.100 p=2.916e-15 p=2.368e-15 (df=69)
stim (main, window start) t(68)= 0.26 F(1)= 0.068 p=0.7946 p=0.7968 (df=18)
interaction 95% CI: [-0.36, +3.57]
HONEST LME (per-animal random slope, day|rat): interaction F(1,22.2)=2.05, p=0.1666
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
model above is the honest learning-rate test -- DF collapses toward the animal count.)
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.1084, slope diff=+1.60)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.