analysis(matlab): formalize cross-study (current vs _f) rate comparison
Add make_crossstudy_compare.m producing crossstudy/result.txt + crossstudy_summary.csv: current control (A2, A2+Naive) and anodal (B2) arms vs the previous study's a2_f/b2_f, on per-animal success RATE, at full range (0-9 both) and attempt-matched (current 0-3 vs previous full). Finding: current > previous for both arms at full range (B2 0.54 vs 0.41 p=0.024; A2+Naive 0.42 vs 0.33 p=0.014) but the gap vanishes/reverses at matched effort (B2 0.33 vs 0.41 n.s.; A2 0.30 vs 0.33 n.s.) -- a practice/ attempts artifact (~2.7x more attempts/session), not a stronger cohort. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
@@ -0,0 +1,104 @@
|
||||
function make_crossstudy_compare()
|
||||
%MAKE_CROSSSTUDY_COMPARE Cross-study current-vs-previous(_f) rate comparisons.
|
||||
% MAKE_CROSSSTUDY_COMPARE() compares the current study's control (Box-A2,
|
||||
% and Box-A2+Naive) and anodal (Box-B2) arms against the previous (Forouzan)
|
||||
% study's control (a2_f) and anodal (b2_f) arms, on the per-animal success
|
||||
% RATE (success/attempts) -- the fair cross-study metric, since the current
|
||||
% protocol runs ~2.7x more attempts per session.
|
||||
%
|
||||
% Two windows are reported for each arm:
|
||||
% full both studies over days 0-9 (equal CALENDAR time)
|
||||
% matched current days 0-3 vs previous full 0-9 (equal cumulative ATTEMPTS;
|
||||
% the previous study's whole span ~= the current study's first
|
||||
% ~3 days by attempt count)
|
||||
%
|
||||
% Each row reports per-animal rate, mean cumulative attempts (the effort being
|
||||
% matched), and honest two-group tests (Welch t, Mann-Whitney U). Writes
|
||||
% analysis/matlab/crossstudy/result.txt and crossstudy_summary.csv.
|
||||
%
|
||||
% Key finding: the current study out-performs the previous one at full range
|
||||
% for BOTH arms, but that disappears (slightly reverses) at matched effort --
|
||||
% a practice/attempts artifact, not a higher-performing cohort.
|
||||
|
||||
thisDir = fileparts(mfilename('fullpath'));
|
||||
outDir = fullfile(thisDir, 'crossstudy');
|
||||
if ~exist(outDir, 'dir'); mkdir(outDir); end
|
||||
|
||||
Tc = tdcs_load_data();
|
||||
Tf = forouzan_load_data();
|
||||
|
||||
% arm | curLabel | curGroups | curWin | prevLabel | prevWin | mode
|
||||
specs = {
|
||||
'Control', 'A2', {'Electrode-Box-A2'}, [0 9], 'a2_f', [0 9], 'full'
|
||||
'Control', 'A2', {'Electrode-Box-A2'}, [0 3], 'a2_f', [0 9], 'matched'
|
||||
'Control', 'A2+Naive', {'Electrode-Box-A2', 'Naive'}, [0 9], 'a2_f', [0 9], 'full'
|
||||
'Control', 'A2+Naive', {'Electrode-Box-A2', 'Naive'}, [0 3], 'a2_f', [0 9], 'matched'
|
||||
'Anodal', 'B2', {'Electrode-Box-B2'}, [0 9], 'b2_f', [0 9], 'full'
|
||||
'Anodal', 'B2', {'Electrode-Box-B2'}, [0 3], 'b2_f', [0 9], 'matched'
|
||||
};
|
||||
|
||||
bar = repmat('=', 1, 92);
|
||||
s = sprintf('%s\nCROSS-STUDY COMPARISON: current vs previous (Forouzan, _f) -- per-animal success RATE\n%s\n', bar, bar);
|
||||
s = [s sprintf(['Metric: per-animal rate = sum(success)/sum(attempts). Rate is used because the\n' ...
|
||||
'current protocol runs ~2.7x more attempts/session, so raw counts are not comparable.\n' ...
|
||||
'Anodal = Box-B2 (contralateral tDCS); control = Box-A2 (ipsilateral sham) +/- Naive.\n' ...
|
||||
'Previous study: b2_f = Anodal, a2_f = Control (validated: b2_f reproduces the paper''s\n' ...
|
||||
'positive tDCS x day interaction). "matched" = current days 0-3 vs previous full 0-9,\n' ...
|
||||
'equalising cumulative attempts (the previous study''s whole span ~= current''s first 3 days).\n\n'])];
|
||||
|
||||
hdr = sprintf('%-8s %-10s %-8s %5s %6s %8s %6s %6s %8s %8s %8s', ...
|
||||
'arm', 'current', 'window', 'nCur', 'rate', 'cumAtt', 'prev', 'rate', 'cumAtt', 'Welch p', 'MWU p');
|
||||
s = [s hdr sprintf('\n') repmat('-', 1, length(hdr)) sprintf('\n')];
|
||||
|
||||
rows = {};
|
||||
for i = 1:size(specs, 1)
|
||||
arm = specs{i, 1}; curLab = specs{i, 2}; curGrp = specs{i, 3};
|
||||
curWin = specs{i, 4}; prevLab = specs{i, 5}; prevWin = specs{i, 6}; mode = specs{i, 7};
|
||||
|
||||
[cR, cA] = localRate(Tc, curGrp, curWin);
|
||||
[pR, pA] = localRate(Tf, {prevLab}, prevWin);
|
||||
[~, welch] = ttest2(cR, pR, 'Vartype', 'unequal');
|
||||
mwu = ranksum(cR, pR);
|
||||
|
||||
winStr = sprintf('%d-%d', curWin(1), curWin(2));
|
||||
s = [s sprintf('%-8s %-10s %-8s %5d %6.3f %8.0f %6s %6.3f %8.0f %8.3f %8.3f\n', ...
|
||||
arm, curLab, winStr, numel(cR), mean(cR), mean(cA), ...
|
||||
prevLab, mean(pR), mean(pA), welch, mwu)]; %#ok<AGROW>
|
||||
rows(end + 1, :) = {arm, curLab, mode, winStr, numel(cR), mean(cR), mean(cA), ...
|
||||
prevLab, mean(pR), mean(pA), welch, mwu}; %#ok<AGROW>
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nINTERPRETATION\n' ...
|
||||
' - FULL range (0-9 both): current > previous for BOTH arms (control and anodal),\n' ...
|
||||
' e.g. B2 0.54 vs b2_f 0.41 (p=0.024); A2+Naive 0.42 vs a2_f 0.33 (p=0.014).\n' ...
|
||||
' - MATCHED effort (current 0-3 vs previous full): the gap disappears and slightly\n' ...
|
||||
' reverses -- B2 0.33 vs b2_f 0.41 (n.s.); A2 0.30 vs a2_f 0.33 (n.s.).\n' ...
|
||||
' => The current study''s apparent superiority is a PRACTICE/ATTEMPTS artifact: its\n' ...
|
||||
' protocol packs ~2.7x more attempts/session, pushing every arm further along the\n' ...
|
||||
' learning curve by a given day. At equal effort the two cohorts perform the same.\n' ...
|
||||
' Note: A2+Naive 0-3 is slightly under-matched on attempts (294 vs ~384); A2-alone\n' ...
|
||||
' 0-3 (361 vs 384) is the cleaner apples-to-apples match.\n'])];
|
||||
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(outDir, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
T = cell2table(rows, 'VariableNames', {'arm', 'current', 'mode', 'window', ...
|
||||
'nCur', 'rateCur', 'cumAttCur', 'prev', 'ratePrev', 'cumAttPrev', 'welchP', 'mwuP'});
|
||||
writetable(T, fullfile(outDir, 'crossstudy_summary.csv'));
|
||||
fprintf('\nWrote %s and crossstudy_summary.csv\n', fullfile(outDir, 'result.txt'));
|
||||
|
||||
end
|
||||
|
||||
% per-animal pooled rate and cumulative attempts over a window
|
||||
function [rate, cumAtt] = localRate(T, groups, win)
|
||||
T = T(ismember(cellstr(T.group), groups) & T.day >= win(1) & T.day <= win(2), :);
|
||||
subj = unique(T.subject);
|
||||
rate = zeros(numel(subj), 1); cumAtt = rate;
|
||||
for i = 1:numel(subj)
|
||||
r = T.subject == subj(i);
|
||||
rate(i) = sum(T.success(r)) / sum(T.total(r));
|
||||
cumAtt(i) = sum(T.total(r));
|
||||
end
|
||||
end
|
||||
Reference in New Issue
Block a user