analysis(matlab): convert previous (Forouzan) study to our long format

Add forouzan_load_data.m (joins raw/Data_Forouzan.csv + raw/Total_Success.csv)
and make_forouzan_data.m, producing analysis/forouzan_reach_data.csv in our
curated column format (subject, group, day, success, total). Group is mapped
Anodal -> b2_f (tDCS), Control -> a2_f (control); Hand ignored; day 0-indexed
(Session 1 = day 0 = paper 'Day 1'); early-stopped sessions (total=0) dropped.

Conversion validated by refitting the paper's own LME on the converted data:
interaction t(227)=2.66, F(1)=7.09, p=0.0083 -- reproduces the paper's
reported t(227)=2.68, F(1)=7.12, p=0.008, N=24 (12 Anodal / 12 Control).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
2026-07-22 14:18:47 -04:00
parent a4dff31772
commit d9b65a3279
3 changed files with 348 additions and 0 deletions
+38
View File
@@ -0,0 +1,38 @@
function make_forouzan_data()
%MAKE_FOROUZAN_DATA Convert the previous (Forouzan) study to our long format.
% MAKE_FOROUZAN_DATA() calls FOROUZAN_LOAD_DATA and writes the curated long
% table to analysis/forouzan_reach_data.csv (same columns as our
% analysis/tdcs_reach_data.csv: subject, group, day, success, total), with
% group in {a2_f (Control), b2_f (Anodal)}. It then prints a descriptives
% summary and, as a conversion sanity-check, fits the paper's LME
% (behavior ~ stim + day + stim:day + (1|rat)) and compares the interaction
% to the paper's reported values (t(227)=2.68, F(1)=7.12, p=0.008, N=24).
thisDir = fileparts(mfilename('fullpath'));
T = forouzan_load_data();
outFile = fullfile(fileparts(thisDir), 'forouzan_reach_data.csv');
writetable(T, outFile);
fprintf('Wrote %s (%d rows)\n\n', outFile, height(T));
fprintf('%-6s %6s %10s %12s %10s %8s\n', ...
'group', 'nRats', 'nSessions', 'mean_success', 'mean_rate', 'max_day');
for g = categories(T.group)'
r = T.group == g{1};
fprintf('%-6s %6d %10d %12.1f %10.3f %8d\n', g{1}, ...
numel(unique(T.subject(r))), sum(r), mean(T.success(r)), ...
sum(T.success(r)) / sum(T.total(r)), max(T.day(r)));
end
% --- conversion sanity-check: the paper's own LME on the converted data ---
tbl = table(T.success, T.day, double(T.group == 'b2_f'), categorical(T.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m);
ii = strcmp(C.Name, 'day:stim'); ai = strcmp(A.Term, 'day:stim');
fprintf(['\nPaper-LME check (behavior ~ stim + day + stim:day + (1|rat), N=%d):\n' ...
' stim x day interaction: t(%d)=%.2f, F(1)=%.2f, p=%.4g\n' ...
' paper reported: t(227)=2.68, F(1)=7.12, p=0.008\n'], ...
numel(unique(tbl.rat)), C.DF(ii), C.tStat(ii), A.FStat(ai), C.pValue(ii));
end