From 1009f263e846af01cfedbc8986ecaf9d20348862 Mon Sep 17 00:00:00 2001 From: Experiments DB Dev Date: Fri, 24 Jul 2026 02:32:43 -0400 Subject: [PATCH] analysis(matlab): rate-variant + vertical x-axis (templates, generators, summary) Refactor variation_{analyze,power,logpower,plot}.m to run both count and rate metrics; plot x-axis labels vertical (xtickangle 90). SUMMARY.csv gains interaction_p_rate / interaction_p_rs_rate. Co-Authored-By: Claude Opus 4.8 --- analysis/matlab/make_variations.m | 5 +- analysis/matlab/variation_analyze.m | 102 +++++++++++++++---------- analysis/matlab/variation_logpower.m | 77 +++++++++---------- analysis/matlab/variation_plot.m | 50 +++++++----- analysis/matlab/variation_power.m | 69 ++++++++--------- analysis/matlab/variations/SUMMARY.csv | 42 +++++----- 6 files changed, 190 insertions(+), 155 deletions(-) diff --git a/analysis/matlab/make_variations.m b/analysis/matlab/make_variations.m index 381bf53..cdf39a4 100644 --- a/analysis/matlab/make_variations.m +++ b/analysis/matlab/make_variations.m @@ -75,7 +75,7 @@ for gi = 1:size(groupings, 1) r = localRun(fullfile(folder, 'analyze.m')); rows(end + 1, :) = {vname, gname, wname, r.nRats, r.nObs, ... r.covEqual, r.interP, r.interPsatt, r.interPrs, r.interEst, ... - r.stimP, r.dayP}; %#ok + r.stimP, r.dayP, r.interPrate, r.interPrsRate}; %#ok fprintf(' %-16s N=%2d obs=%4d interaction p=%.4g (%+.2f) %s\n', ... vname, r.nRats, r.nObs, r.interP, r.interEst, ... localTern(r.covEqual, 'equal-cov', 'UNEQUAL-cov')); @@ -89,7 +89,8 @@ end % Top-level index of all variations. S = cell2table(rows, 'VariableNames', {'variation', 'grouping', 'window', ... 'nRats', 'nObs', 'covEqual', 'interaction_p', 'interaction_p_satt', ... - 'interaction_p_rs', 'interaction_est', 'stim_p', 'day_p'}); + 'interaction_p_rs', 'interaction_est', 'stim_p', 'day_p', ... + 'interaction_p_rate', 'interaction_p_rs_rate'}); writetable(S, fullfile(root, 'SUMMARY.csv')); fprintf('\nWrote %d variations under %s\n(index: variations/SUMMARY.csv)\n', ... size(rows, 1), root); diff --git a/analysis/matlab/variation_analyze.m b/analysis/matlab/variation_analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variation_analyze.m +++ b/analysis/matlab/variation_analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % 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.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (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 +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, 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 +As = anova(m, 'DFMethod', 'satterthwaite'); -% 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 @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f 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))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); 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 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 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 +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)'), '', ''); -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 = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; 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))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res 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))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\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 @@ -82,17 +105,18 @@ else 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')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) 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); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variation_logpower.m b/analysis/matlab/variation_logpower.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variation_logpower.m +++ b/analysis/matlab/variation_logpower.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variation_plot.m b/analysis/matlab/variation_plot.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variation_plot.m +++ b/analysis/matlab/variation_plot.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variation_power.m b/analysis/matlab/variation_power.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variation_power.m +++ b/analysis/matlab/variation_power.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/SUMMARY.csv b/analysis/matlab/variations/SUMMARY.csv index d156d26..668b1a3 100644 --- a/analysis/matlab/variations/SUMMARY.csv +++ b/analysis/matlab/variations/SUMMARY.csv @@ -1,21 +1,21 @@ -variation,grouping,window,nRats,nObs,covEqual,interaction_p,interaction_p_satt,interaction_p_rs,interaction_est,stim_p,day_p -unmerge_d0_10,unmerge,d0_10,6,61,1,0.142831730250263,0.142797566084653,0.247946685371652,1.55766564374789,0.497389134925754,2.83166776054496e-14 -unmerge_d0_13,unmerge,d0_13,6,70,1,0.1037612192755,0.103490077847449,0.222686826261486,1.54665472676189,0.343602740986316,3.0346358823857e-13 -unmerge_d0_5,unmerge,d0_5,6,36,1,0.015149263054839,0.0154996899553134,0.133145634512884,5.36190476190475,0.688839435920596,1.57055423419133e-07 -unmerge_d6_10,unmerge,d6_10,5,25,1,0.485980598616448,0.486366669144796,0.554919025308352,1.76666666666666,0.0975215172907211,0.227259314093653 -unmerge_d6_13,unmerge,d6_13,5,34,1,0.0971009786566259,0.0970128901807668,0.13816578257917,2.26696773140328,0.101395665857404,0.117493504965036 -right_only_d0_10,right_only,d0_10,7,72,1,0.072180139606096,0.0721794737454714,0.124574209202325,1.73004431141938,0.715879551990348,5.94289273347193e-16 -right_only_d0_13,right_only,d0_13,7,84,1,0.0480416099228058,0.0478799998667326,0.102857124294903,1.7001712792199,0.528517264593511,1.17078898575779e-14 -right_only_d0_5,right_only,d0_5,7,42,1,0.00589730538358453,0.00612750828678172,0.0876259332262465,5.30714285714286,0.469556667463413,1.45827052706827e-08 -right_only_d6_10,right_only,d6_10,6,30,1,0.288979620167748,0.289798539517805,0.378365014630273,2.425,0.119238243812818,0.200957769624237 -right_only_d6_13,right_only,d6_13,6,42,1,0.0258497485171371,0.0259134063416617,0.0369743190514673,2.75792254325755,0.109557181635453,0.0976080018855909 -naive_a2_d0_10,naive_a2,d0_10,10,104,1,0.125891264174292,0.126071775330177,0.30932733005532,1.43371876298297,0.141836960529428,1.59796354184479e-27 -naive_a2_d0_13,naive_a2,d0_13,10,124,1,0.0397967517409217,0.0398599319703312,0.0578149015961064,1.61366204723722,0.153633906173184,2.5469330031931e-28 -naive_a2_d0_5,naive_a2,d0_5,10,59,1,0.00107331575782455,0.00115218904921896,0.028152310883297,6.96223661591522,0.960550975833492,1.1254222201539e-09 -naive_a2_d6_10,naive_a2,d6_10,9,45,1,0.449922139297732,0.450525162691322,0.522263670801849,-1.8,0.00929038977362001,8.02730473236241e-05 -naive_a2_d6_13,naive_a2,d6_13,9,65,1,0.439759713746994,0.439915577313227,0.563543413702517,1.15527349799111,0.0165675905421826,0.00123453301903328 -naive_boxa_d0_10,naive_boxa,d0_10,11,115,1,0.194116831474294,0.194284662755081,0.383551190839722,1.17792772470162,0.124842156542362,1.7431506758907e-32 -naive_boxa_d0_13,naive_boxa,d0_13,11,138,1,0.0848056366333164,0.0848864932608599,0.140324990567599,1.30637615965008,0.140502409490992,1.21292542612342e-34 -naive_boxa_d0_5,naive_boxa,d0_5,11,65,1,0.00254470765200637,0.00267934065308424,0.0356145839564805,6.40329932104724,0.881419104187055,2.50618640106882e-11 -naive_boxa_d6_10,naive_boxa,d6_10,10,50,1,0.46817073850354,0.468724498567948,0.534407300553877,-1.61904761904762,0.014689379915428,1.87688778586107e-05 -naive_boxa_d6_13,naive_boxa,d6_13,10,73,1,0.530807487326757,0.530942127564996,0.696589180707779,0.873689202983884,0.0303916991304334,6.44631091768978e-05 +variation,grouping,window,nRats,nObs,covEqual,interaction_p,interaction_p_satt,interaction_p_rs,interaction_est,stim_p,day_p,interaction_p_rate,interaction_p_rs_rate +unmerge_d0_10,unmerge,d0_10,6,61,1,0.142831730250263,0.142797566084653,0.247946685371652,1.55766564374789,0.497389134925754,2.83166776054496e-14,0.101273272943887,0.113635949369058 +unmerge_d0_13,unmerge,d0_13,6,70,1,0.1037612192755,0.103490077847449,0.222686826261486,1.54665472676189,0.343602740986316,3.0346358823857e-13,0.0729395844508811,0.0810848229835278 +unmerge_d0_5,unmerge,d0_5,6,36,1,0.015149263054839,0.0154996899553134,0.133145634512884,5.36190476190475,0.688839435920596,1.57055423419133e-07,0.00312728156805176,0.0378990615877136 +unmerge_d6_10,unmerge,d6_10,5,25,1,0.485980598616448,0.486366669144796,0.554919025308352,1.76666666666666,0.0975215172907211,0.227259314093653,0.889539569780549,0.921910870826714 +unmerge_d6_13,unmerge,d6_13,5,34,1,0.0971009786566259,0.0970128901807668,0.13816578257917,2.26696773140328,0.101395665857404,0.117493504965036,0.168410922832656,0.874930263016369 +right_only_d0_10,right_only,d0_10,7,72,1,0.072180139606096,0.0721794737454714,0.124574209202325,1.73004431141938,0.715879551990348,5.94289273347193e-16,0.0406170660292579,0.0396836748954027 +right_only_d0_13,right_only,d0_13,7,84,1,0.0480416099228058,0.0478799998667326,0.102857124294903,1.7001712792199,0.528517264593511,1.17078898575779e-14,0.0182494380376794,0.01709229589027 +right_only_d0_5,right_only,d0_5,7,42,1,0.00589730538358453,0.00612750828678172,0.0876259332262465,5.30714285714286,0.469556667463413,1.45827052706827e-08,0.000673572060461354,0.0180977446376148 +right_only_d6_10,right_only,d6_10,6,30,1,0.288979620167748,0.289798539517805,0.378365014630273,2.425,0.119238243812818,0.200957769624237,0.514255837007347,0.649543797033249 +right_only_d6_13,right_only,d6_13,6,42,1,0.0258497485171371,0.0259134063416617,0.0369743190514673,2.75792254325755,0.109557181635453,0.0976080018855909,0.029980908747733,0.649212848234672 +naive_a2_d0_10,naive_a2,d0_10,10,104,1,0.125891264174292,0.126071775330177,0.30932733005532,1.43371876298297,0.141836960529428,1.59796354184479e-27,0.095917013160988,0.178648474821737 +naive_a2_d0_13,naive_a2,d0_13,10,124,1,0.0397967517409217,0.0398599319703312,0.0578149015961064,1.61366204723722,0.153633906173184,2.5469330031931e-28,0.025894419150261,0.0266660380156418 +naive_a2_d0_5,naive_a2,d0_5,10,59,1,0.00107331575782455,0.00115218904921896,0.028152310883297,6.96223661591522,0.960550975833492,1.1254222201539e-09,0.000395858261274059,0.00563668931877035 +naive_a2_d6_10,naive_a2,d6_10,9,45,1,0.449922139297732,0.450525162691322,0.522263670801849,-1.8,0.00929038977362001,8.02730473236241e-05,0.143377545072213,0.23692983335782 +naive_a2_d6_13,naive_a2,d6_13,9,65,1,0.439759713746994,0.439915577313227,0.563543413702517,1.15527349799111,0.0165675905421826,0.00123453301903328,0.712698437296026,0.839317422332943 +naive_boxa_d0_10,naive_boxa,d0_10,11,115,1,0.194116831474294,0.194284662755081,0.383551190839722,1.17792772470162,0.124842156542362,1.7431506758907e-32,0.174712550461141,0.268147954364157 +naive_boxa_d0_13,naive_boxa,d0_13,11,138,1,0.0848056366333164,0.0848864932608599,0.140324990567599,1.30637615965008,0.140502409490992,1.21292542612342e-34,0.0717859200229646,0.121612217594298 +naive_boxa_d0_5,naive_boxa,d0_5,11,65,1,0.00254470765200637,0.00267934065308424,0.0356145839564805,6.40329932104724,0.881419104187055,2.50618640106882e-11,0.0010002309013673,0.00912728188012594 +naive_boxa_d6_10,naive_boxa,d6_10,10,50,1,0.46817073850354,0.468724498567948,0.534407300553877,-1.61904761904762,0.014689379915428,1.87688778586107e-05,0.141244057488063,0.219772818577224 +naive_boxa_d6_13,naive_boxa,d6_13,10,73,1,0.530807487326757,0.530942127564996,0.696589180707779,0.873689202983884,0.0303916991304334,6.44631091768978e-05,0.852260403603235,0.983649216326651