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 <noreply@anthropic.com>

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<AGROW>
+                r.stimP, r.dayP, r.interPrate, r.interPrsRate}; %#ok<AGROW>
             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);

analysis/matlab/variation_analyze.m

@@ -1,30 +1,58 @@
-% Variation analysis -- the paper's linear mixed model on the successful-reach
+% Variation analysis -- the paper's linear mixed model on this folder's data,
-% COUNT, fit on this folder's curated data subset.
+% 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 treatment / 0 control; day = training day within window (0 =
-%     stim     = 1 for the treatment group(s), 0 for the control group(s)
+%     first analyzed day); rat = subject (random intercept).
-%     day      = training day within this window (0 = first analyzed day)
+% For the interaction we report residual DF, Satterthwaite DF, and the honest
-%     rat      = subject (random intercept)
+% per-animal random-slope test. Self-contained: reads data.csv beside this
-%
+% script. Run headless with:  matlab -batch "analyze"
-% Self-contained: reads data.csv beside this script and writes result.txt.
+% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
-% 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
+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
-
+if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
-ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
+elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
-if pI >= 0.05
+else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
-    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 = 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 = success COUNT)\n')];
+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('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\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('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');
+fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
-fprintf(fid, '%s', s);
+fprintf(fid, '%s', s); fclose(fid);
-fclose(fid);
+end
-
-% 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);

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
+% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
-% "Day 1"), so log(day + 1) reproduces their transform exactly.
+% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
-%
+% refits that model, reports the interaction (residual / Satterthwaite / honest
-% This script (a) refits that log-day model on data.csv, (b) reports the
+% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
-% interaction (residual DF, Satterthwaite DF, and the honest per-animal
+% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
-% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
+% Run:  matlab -batch "logpowersim"
-% 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"
 % (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');
+localLogPower(D, 'count', here, vname);
-rng(1);
+localLogPower(D, 'rate', here, vname);
 
-logday = log(D.day + 1);                       % 0-indexed day -> their log(1-indexed day)
+% ---------------------------------------------------------------- per metric
-tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
+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 = 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)   (success COUNT)\n')];
+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')];
+    s = [s sprintf('\nInsufficient data for this analysis.\n')];
-    localFinish(s, here); warning(warnState); return
+    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<AGROW>
+    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))]; %#ok<AGROW>
+    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<AGROW>
+        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' ...
+s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
-    'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
+localFinish(s, here, suffix);
-
-localFinish(s, here);
 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
 

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 /
+% Produces TWO figures from this folder's data.csv:
-% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
+%   learning_curve.png       # successes (count)   per training day
-% this folder's data.csv and saving learning_curve.png. The per-group N is read
+%   learning_curve_rate.png  success rate (success/attempts) per training day
-% from the data (each variation pools different groups), so the legend shows the
+% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
-% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
+% 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
+days = unique(D.day);
-xd = days + 1;                     % plot as 1-indexed training day (paper axis)
+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)
+localPlot(D, days, xd, 'count', '# successes', ...
-[Mc, Sc, nc] = localCurve(D, 0, days);   % control (stim = 0)
+    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');
+xlabel('training day'); ylabel(ylab);
-ylabel('# successes');
 title(vname, 'Interpreter', 'none');
 set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
-
+xtickangle(90);                  % vertical x-axis tick labels
-outFile = fullfile(here, 'learning_curve.png');
 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, metric)
-function [M, S, n] = localCurve(D, stimVal, days)
+%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
-%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
 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');

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,
+% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
-% and residual SD generate NREP synthetic datasets at each rats-per-group N and
+% generate NREP synthetic datasets at each rats-per-group N and each true-effect
-% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
+% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
-% scored at alpha = 0.05 two ways:
+%   per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
-%   per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
-%                honest power, matching the random-slope / per-animal inference)
 %   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];
+if strcmp(metric, 'rate')
-EFFMULS = [1 0.5];
+    D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
-NREP = 120;
+else
-
+    beh = D.success; mlabel = '# successes (count)'; suffix = '';
-warnState = warning('off', 'all');
+end
-rng(1);
+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 = sprintf('%s\nPOWER SIMULATION -- %s   [metric: %s]\n%s\n', bar, vname, mlabel, bar);
-s = [s sprintf('model: %s   (success COUNT; day within-window)\n', FORMULA)];
+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')];
+    s = [s sprintf('\nInsufficient data for a power simulation.\n')];
-    localFinish(s, here); warning(warnState); return
+    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<AGROW>
+    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', 'N/group')]; %#ok<AGROW>
+    s = [s sprintf('  %-8s | per-animal power | LME power\n  %s\n', 'N/group', repmat('-', 1, 42))];
-    s = [s sprintf('  %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
     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<AGROW>
+        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; LME is optimistic (obs-level DF).\n')];
-s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
+localFinish(s, here, suffix);
-    'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
-
-localFinish(s, here);
 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
 

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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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
+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