analysis(matlab): rate + count outputs, variation batch 2

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
2026-07-24 02:32:52 -04:00
parent 0a438a4649
commit 6ba21a1a35
72 changed files with 1874 additions and 896 deletions
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 36 KiB

After

Width:  |  Height:  |  Size: 35 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 38 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- matched_current_d0_3 LOG-DAY MODEL + COHEN'S f + POWER -- matched_current_d0_3 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05 observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,20)=7.043 p(resid)=0.01524 p(Satt)=0.0139 (df
honest per-animal random slope (log-day): F(1,18.5)=7.24 p=0.01475 honest per-animal random slope (log-day): F(1,18.5)=7.24 p=0.01475
Cohen's f (interaction, partial eta^2=0.102) = 0.336 (medium; f: .10 small, .25 medium, .40 large) Cohen's f (interaction, partial eta^2=0.102) = 0.336 (medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+20.79, ratSD=0.00, resSD=9.99) --- --- power simulation (log-day ground truth: stim:logday=+20.79, ratSD=0, resSD=9.991) ---
true stim:log(day) = +20.79 (100% of observed) true stim:log(day) = +20.79 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.102) = 0.336 (medium; f: .10 small, .25
16 | 1.00 | 1.00 16 | 1.00 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:log(day) = +10.40 (50% of observed) true stim:log(day) = +10.4 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.11 | 0.23 <- observed 3 | 0.11 | 0.23 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.102) = 0.336 (medium; f: .10 small, .25
16 | 0.79 | 0.84 16 | 0.79 | 0.84
24 | 0.95 | 0.97 24 | 0.95 | 0.97
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- matched_current_d0_3 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,20)=7.974 p(resid)=0.01049 p(Satt)=0.01125 (df=18)
honest per-animal random slope (log-day): F(1,18.3)=7.99 p=0.01106
Cohen's f (interaction, partial eta^2=0.157) = 0.431 (large; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+0.1785, ratSD=0.01747, resSD=0.08062) ---
true stim:log(day) = +0.1785 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.42 | 0.82 <- observed
5 | 0.88 | 0.95
8 | 1.00 | 0.99
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +0.08926 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.26 <- observed
5 | 0.39 | 0.50
8 | 0.53 | 0.63
12 | 0.82 | 0.82
16 | 0.87 | 0.90
24 | 0.97 | 0.98
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- matched_current_d0_3 POWER SIMULATION -- matched_current_d0_3 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05 observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
ground truth: stim:day=+10.20/day, rat SD=2.92, residual SD=8.54, days=4 ground truth: stim:day=+10.2/day, rat SD=2.916, residual SD=8.54, days=4
true stim:day interaction = +10.20 (100% of observed) true stim:day interaction = +10.2 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.55 | 0.91 <- observed 3 | 0.55 | 0.91 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=+10.20/day, rat SD=2.92, residual SD=8.54, days=4
16 | 1.00 | 1.00 16 | 1.00 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:day interaction = +5.10 (50% of observed) true stim:day interaction = +5.1 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.11 | 0.33 <- observed 3 | 0.11 | 0.33 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=+10.20/day, rat SD=2.92, residual SD=8.54, days=4
16 | 0.97 | 0.97 16 | 0.97 | 0.97
24 | 1.00 | 1.00 24 | 1.00 | 1.00
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- matched_current_d0_3 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
ground truth: stim:day=+0.07723/day, rat SD=0.01983, residual SD=0.07841, days=4
true stim:day interaction = +0.07723 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.39 | 0.78 <- observed
5 | 0.86 | 0.94
8 | 0.98 | 0.99
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:day interaction = +0.03861 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.21 <- observed
5 | 0.35 | 0.47
8 | 0.53 | 0.62
12 | 0.78 | 0.79
16 | 0.82 | 0.87
24 | 0.97 | 0.97
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: matched_current_d0_3 VARIATION: matched_current_d0_3 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(20)= 3.27 F(1)= 10.700 p=0.003822 p=0.004244 (df=18) stim x day (interaction) t(20)= 3.27 F(1)= 10.700 p=0.003822 p=0.004244 (df=18)
day (learning) t(20)= 2.89 F(1)= 8.337 p=0.009105 p=0.009807 (df=18) day (learning) t(20)= 2.89 F(1)= 8.337 p=0.009105 p=0.009807 (df=18)
stim (main, window start) t(20)= -1.42 F(1)= 2.025 p=0.1701 p=0.1703 (df=20) stim (main, window start) t(20)= -1.42 F(1)= 2.025 p=0.1701 p=0.1703 (df=20)
interaction 95% CI: [+3.70, +16.70] interaction 95% CI: [+3.695, +16.7]
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.8)=9.76, p=0.01107 HONEST LME (per-animal random slope, day|rat): interaction F(1,9.8)=9.76, p=0.01107
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) model above is the honest learning-rate test -- DF collapses toward the animal count.)
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.003822, slope diff=+10.20) INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.003822, slope diff=+10.2)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: matched_current_d0_3 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2
N = 6 rats, 24 sessions raw day coverage: treat 0..3, control 0..3
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 24
Fixed effects coefficients 4
Random effects coefficients 6
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-40.724 -33.656 26.362 -52.724
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 0.23952 0.039567 6.0535 20 6.443e-06
{'day' } 0.027914 0.020245 1.3788 20 0.18317
{'stim' } -0.089781 0.055956 -1.6045 20 0.12428
{'day:stim' } 0.077228 0.02863 2.6974 20 0.013857
Lower Upper
0.15698 0.32205
-0.014316 0.070144
-0.2065 0.026941
0.017506 0.13695
Random effects covariance parameters (95% CIs):
Group: rat (6 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.019829
Lower Upper
0.00093276 0.42155
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.078408 0.056558 0.1087
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(20)= 2.70 F(1)= 7.276 p=0.01386 p=0.01473 (df=18)
day (learning) t(20)= 1.38 F(1)= 1.901 p=0.1832 p=0.1848 (df=18)
stim (main, window start) t(20)= -1.60 F(1)= 2.574 p=0.1243 p=0.1235 (df=21)
interaction 95% CI: [+0.01751, +0.1369]
HONEST LME (per-animal random slope, day|rat): interaction F(1,18.0)=7.28, p=0.01473
(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 SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.01386, slope diff=+0.07723)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 41 KiB

After

Width:  |  Height:  |  Size: 41 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 42 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_10 LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_10 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
@@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,100)=4.633 p(resid)=0.03378 p(Satt)=0.03391 (
honest per-animal random slope (log-day): F(1,7.2)=2.24 p=0.177 honest per-animal random slope (log-day): F(1,7.2)=2.24 p=0.177
Cohen's f (interaction, partial eta^2=0.010) = 0.101 (small-medium; f: .10 small, .25 medium, .40 large) Cohen's f (interaction, partial eta^2=0.010) = 0.101 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+9.59, ratSD=8.17, resSD=14.84) --- --- power simulation (log-day ground truth: stim:logday=+9.586, ratSD=8.17, resSD=14.84) ---
true stim:log(day) = +9.59 (100% of observed) true stim:log(day) = +9.586 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.23 | 0.44 <- observed 3 | 0.23 | 0.44 <- observed
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.010) = 0.101 (small-medium; f: .10 smal
16 | 0.99 | 1.00 16 | 0.99 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:log(day) = +4.79 (50% of observed) true stim:log(day) = +4.793 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.07 | 0.10 <- observed 3 | 0.07 | 0.10 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.010) = 0.101 (small-medium; f: .10 smal
16 | 0.62 | 0.64 16 | 0.62 | 0.64
24 | 0.75 | 0.78 24 | 0.75 | 0.78
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_10 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,98)=6.255 p(resid)=0.01404 p(Satt)=0.01413 (df=93)
honest per-animal random slope (log-day): F(1,6.5)=3.64 p=0.1012
Cohen's f (interaction, partial eta^2=0.019) = 0.138 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+0.07796, ratSD=0.04949, resSD=0.1018) ---
true stim:log(day) = +0.07796 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.25 | 0.57 <- observed
5 | 0.71 | 0.84
8 | 0.92 | 0.95
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +0.03898 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.09 | 0.15 <- observed
5 | 0.14 | 0.26
8 | 0.41 | 0.46
12 | 0.52 | 0.56
16 | 0.78 | 0.78
24 | 0.88 | 0.89
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- naive_a2_d0_10 POWER SIMULATION -- naive_a2_d0_10 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+1.43/day, rat SD=8.75, residual SD=13.78, days=11 ground truth: stim:day=+1.434/day, rat SD=8.753, residual SD=13.78, days=11
true stim:day interaction = +1.43 (100% of observed) true stim:day interaction = +1.434 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.09 | 0.23 <- observed 3 | 0.09 | 0.23 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.43/day, rat SD=8.75, residual SD=13.78, days=11
16 | 0.84 | 0.84 16 | 0.84 | 0.84
24 | 0.95 | 0.94 24 | 0.95 | 0.94
true stim:day interaction = +0.72 (50% of observed) true stim:day interaction = +0.7169 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.07 | 0.09 <- observed 3 | 0.07 | 0.09 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.43/day, rat SD=8.75, residual SD=13.78, days=11
16 | 0.38 | 0.36 16 | 0.38 | 0.36
24 | 0.47 | 0.47 24 | 0.47 | 0.47
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- naive_a2_d0_10 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+0.01123/day, rat SD=0.05096, residual SD=0.098, days=11
true stim:day interaction = +0.01123 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.28 <- observed
5 | 0.39 | 0.47
8 | 0.54 | 0.59
12 | 0.84 | 0.84
16 | 0.89 | 0.93
24 | 0.97 | 0.98
true stim:day interaction = +0.005617 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.08 | 0.10 <- observed
5 | 0.07 | 0.12
8 | 0.18 | 0.19
12 | 0.28 | 0.28
16 | 0.43 | 0.40
24 | 0.56 | 0.56
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: naive_a2_d0_10 VARIATION: naive_a2_d0_10 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(100)= 1.54 F(1)= 2.382 p=0.1259 p=0.1261 (df=95) stim x day (interaction) t(100)= 1.54 F(1)= 2.382 p=0.1259 p=0.1261 (df=95)
day (learning) t(100)= 15.09 F(1)=227.710 p=1.598e-27 p=4.059e-27 (df=96) day (learning) t(100)= 15.09 F(1)=227.710 p=1.598e-27 p=4.059e-27 (df=96)
stim (main, window start) t(100)= 1.48 F(1)= 2.192 p=0.1418 p=0.1543 (df=20) stim (main, window start) t(100)= 1.48 F(1)= 2.192 p=0.1418 p=0.1543 (df=20)
interaction 95% CI: [-0.41, +3.28] interaction 95% CI: [-0.4093, +3.277]
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.0)=1.16, p=0.3093 HONEST LME (per-animal random slope, day|rat): interaction F(1,9.0)=1.16, p=0.3093
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) 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.1259, slope diff=+1.43) INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.1259, slope diff=+1.434)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: naive_a2_d0_10 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2, Naive
N = 10 rats, 102 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 102
Fixed effects coefficients 4
Random effects coefficients 10
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-159.21 -143.46 85.607 -171.21
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 0.19588 0.02981 6.571 98 2.4247e-09
{'day' } 0.04169 0.0039441 10.57 98 7.0294e-18
{'stim' } 0.05903 0.05266 1.121 98 0.26504
{'day:stim' } 0.011235 0.0066828 1.6812 98 0.095917
Lower Upper
0.13672 0.25504
0.033863 0.049517
-0.045471 0.16353
-0.002027 0.024497
Random effects covariance parameters (95% CIs):
Group: rat (10 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.050962
Lower Upper
0.027158 0.095633
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.098003 0.084761 0.11331
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(98)= 1.68 F(1)= 2.826 p=0.09592 p=0.09609 (df=93)
day (learning) t(98)= 10.57 F(1)=111.730 p=7.029e-18 p=9.281e-18 (df=96)
stim (main, window start) t(98)= 1.12 F(1)= 1.257 p=0.265 p=0.2735 (df=24)
interaction 95% CI: [-0.002027, +0.0245]
HONEST LME (per-animal random slope, day|rat): interaction F(1,8.4)=2.15, p=0.1786
(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.09592, slope diff=+0.01123)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 42 KiB

After

Width:  |  Height:  |  Size: 42 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 40 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_13 LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_13 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
@@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,120)=6.786 p(resid)=0.01035 p(Satt)=0.01039 (
honest per-animal random slope (log-day): F(1,5.5)=5.06 p=0.06966 honest per-animal random slope (log-day): F(1,5.5)=5.06 p=0.06966
Cohen's f (interaction, partial eta^2=0.011) = 0.104 (small-medium; f: .10 small, .25 medium, .40 large) Cohen's f (interaction, partial eta^2=0.011) = 0.104 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+9.94, ratSD=8.36, resSD=14.40) --- --- power simulation (log-day ground truth: stim:logday=+9.939, ratSD=8.359, resSD=14.4) ---
true stim:log(day) = +9.94 (100% of observed) true stim:log(day) = +9.939 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.27 | 0.71 <- observed 3 | 0.27 | 0.71 <- observed
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.011) = 0.104 (small-medium; f: .10 smal
16 | 1.00 | 1.00 16 | 1.00 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:log(day) = +4.97 (50% of observed) true stim:log(day) = +4.969 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.12 | 0.25 <- observed 3 | 0.12 | 0.25 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.011) = 0.104 (small-medium; f: .10 smal
16 | 0.80 | 0.85 16 | 0.80 | 0.85
24 | 0.92 | 0.93 24 | 0.92 | 0.93
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_13 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,118)=8.542 p(resid)=0.00416 p(Satt)=0.00418 (df=115)
honest per-animal random slope (log-day): F(1,5.1)=7.00 p=0.04487
Cohen's f (interaction, partial eta^2=0.019) = 0.141 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+0.07789, ratSD=0.04775, resSD=0.09872) ---
true stim:log(day) = +0.07789 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.38 | 0.80 <- observed
5 | 0.84 | 0.93
8 | 1.00 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +0.03895 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.13 | 0.31 <- observed
5 | 0.31 | 0.38
8 | 0.59 | 0.63
12 | 0.75 | 0.78
16 | 0.92 | 0.92
24 | 0.95 | 0.95
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- naive_a2_d0_13 POWER SIMULATION -- naive_a2_d0_13 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+1.61/day, rat SD=8.46, residual SD=15.26, days=14 ground truth: stim:day=+1.614/day, rat SD=8.46, residual SD=15.26, days=14
true stim:day interaction = +1.61 (100% of observed) true stim:day interaction = +1.614 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.22 | 0.53 <- observed 3 | 0.22 | 0.53 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.61/day, rat SD=8.46, residual SD=15.26, days=14
16 | 0.98 | 0.98 16 | 0.98 | 0.98
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:day interaction = +0.81 (50% of observed) true stim:day interaction = +0.8068 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.12 | 0.21 <- observed 3 | 0.12 | 0.21 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.61/day, rat SD=8.46, residual SD=15.26, days=14
16 | 0.64 | 0.72 16 | 0.64 | 0.72
24 | 0.78 | 0.75 24 | 0.78 | 0.75
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- naive_a2_d0_13 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+0.01209/day, rat SD=0.0442, residual SD=0.1045, days=14
true stim:day interaction = +0.01209 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.27 | 0.64 <- observed
5 | 0.62 | 0.75
8 | 0.88 | 0.93
12 | 0.98 | 1.00
16 | 0.99 | 0.99
24 | 1.00 | 1.00
true stim:day interaction = +0.006044 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.13 | 0.24 <- observed
5 | 0.21 | 0.24
8 | 0.46 | 0.47
12 | 0.53 | 0.55
16 | 0.70 | 0.77
24 | 0.82 | 0.81
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: naive_a2_d0_13 VARIATION: naive_a2_d0_13 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(120)= 2.08 F(1)= 4.320 p=0.0398 p=0.03986 (df=117) stim x day (interaction) t(120)= 2.08 F(1)= 4.320 p=0.0398 p=0.03986 (df=117)
day (learning) t(120)= 14.57 F(1)=212.424 p=2.547e-28 p=3.246e-28 (df=119) day (learning) t(120)= 14.57 F(1)=212.424 p=2.547e-28 p=3.246e-28 (df=119)
stim (main, window start) t(120)= 1.44 F(1)= 2.062 p=0.1536 p=0.1655 (df=21) stim (main, window start) t(120)= 1.44 F(1)= 2.062 p=0.1536 p=0.1655 (df=21)
interaction 95% CI: [+0.08, +3.15] interaction 95% CI: [+0.07651, +3.151]
HONEST LME (per-animal random slope, day|rat): interaction F(1,41.7)=3.81, p=0.05781 HONEST LME (per-animal random slope, day|rat): interaction F(1,41.7)=3.81, p=0.05781
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) model above is the honest learning-rate test -- DF collapses toward the animal count.)
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.0398, slope diff=+1.61) INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.0398, slope diff=+1.614)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: naive_a2_d0_13 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2, Naive
N = 10 rats, 122 sessions raw day coverage: treat 0..13, control 0..13
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 122
Fixed effects coefficients 4
Random effects coefficients 10
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-181.33 -164.5 96.663 -193.33
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 0.23325 0.027564 8.4623 118 8.2776e-14
{'day' } 0.031488 0.0030917 10.184 118 7.3141e-18
{'stim' } 0.055421 0.048873 1.134 118 0.2591
{'day:stim' } 0.012089 0.0053578 2.2563 118 0.025894
Lower Upper
0.17867 0.28783
0.025365 0.03761
-0.041361 0.1522
0.0014788 0.022699
Random effects covariance parameters (95% CIs):
Group: rat (10 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.044199
Lower Upper
0.022679 0.086139
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.10453 0.091674 0.1192
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(118)= 2.26 F(1)= 5.091 p=0.02589 p=0.02593 (df=116)
day (learning) t(118)= 10.18 F(1)=103.724 p=7.314e-18 p=6.719e-18 (df=119)
stim (main, window start) t(118)= 1.13 F(1)= 1.286 p=0.2591 p=0.2666 (df=27)
interaction 95% CI: [+0.001479, +0.0227]
HONEST LME (per-animal random slope, day|rat): interaction F(1,97.3)=5.06, p=0.02667
(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 SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.02589, slope diff=+0.01209)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 35 KiB

After

Width:  |  Height:  |  Size: 35 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 36 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_5 LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_5 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,55)=9.364 p(resid)=0.003417 p(Satt)=0.003585
honest per-animal random slope (log-day): F(1,13.5)=7.67 p=0.01548 honest per-animal random slope (log-day): F(1,13.5)=7.67 p=0.01548
Cohen's f (interaction, partial eta^2=0.057) = 0.245 (small-medium; f: .10 small, .25 medium, .40 large) Cohen's f (interaction, partial eta^2=0.057) = 0.245 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+19.33, ratSD=9.47, resSD=13.37) --- --- power simulation (log-day ground truth: stim:logday=+19.33, ratSD=9.472, resSD=13.37) ---
true stim:log(day) = +19.33 (100% of observed) true stim:log(day) = +19.33 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.057) = 0.245 (small-medium; f: .10 smal
16 | 1.00 | 1.00 16 | 1.00 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:log(day) = +9.67 (50% of observed) true stim:log(day) = +9.665 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.12 | 0.28 <- observed 3 | 0.12 | 0.28 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.057) = 0.245 (small-medium; f: .10 smal
16 | 0.83 | 0.84 16 | 0.83 | 0.84
24 | 0.95 | 0.96 24 | 0.95 | 0.96
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_5 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,53)=14.105 p(resid)=0.0004314 p(Satt)=0.0004753 (df=47)
honest per-animal random slope (log-day): F(1,11.9)=12.05 p=0.004674
Cohen's f (interaction, partial eta^2=0.124) = 0.376 (medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+0.1804, ratSD=0.05407, resSD=0.09942) ---
true stim:log(day) = +0.1804 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.53 | 0.89 <- observed
5 | 0.92 | 0.97
8 | 1.00 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +0.09019 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.13 | 0.38 <- observed
5 | 0.44 | 0.54
8 | 0.72 | 0.73
12 | 0.78 | 0.88
16 | 0.94 | 0.97
24 | 1.00 | 1.00
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- naive_a2_d0_5 POWER SIMULATION -- naive_a2_d0_5 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05 observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+6.96/day, rat SD=9.64, residual SD=12.11, days=6 ground truth: stim:day=+6.962/day, rat SD=9.643, residual SD=12.11, days=6
true stim:day interaction = +6.96 (100% of observed) true stim:day interaction = +6.962 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.43 | 0.85 <- observed 3 | 0.43 | 0.85 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=+6.96/day, rat SD=9.64, residual SD=12.11, days=6
16 | 1.00 | 1.00 16 | 1.00 | 1.00
24 | 1.00 | 1.00 24 | 1.00 | 1.00
true stim:day interaction = +3.48 (50% of observed) true stim:day interaction = +3.481 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.12 | 0.31 <- observed 3 | 0.12 | 0.31 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=+6.96/day, rat SD=9.64, residual SD=12.11, days=6
16 | 0.90 | 0.89 16 | 0.90 | 0.89
24 | 0.97 | 0.99 24 | 0.97 | 0.99
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- naive_a2_d0_5 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
ground truth: stim:day=+0.06206/day, rat SD=0.05458, residual SD=0.09694, days=6
true stim:day interaction = +0.06206 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.53 | 0.90 <- observed
5 | 0.93 | 0.98
8 | 1.00 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:day interaction = +0.03103 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.35 <- observed
5 | 0.44 | 0.52
8 | 0.66 | 0.72
12 | 0.83 | 0.88
16 | 0.94 | 0.95
24 | 0.99 | 0.99
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: naive_a2_d0_5 VARIATION: naive_a2_d0_5 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(55)= 3.45 F(1)= 11.924 p=0.001073 p=0.001152 (df=49) stim x day (interaction) t(55)= 3.45 F(1)= 11.924 p=0.001073 p=0.001152 (df=49)
day (learning) t(55)= 7.32 F(1)= 53.591 p=1.125e-09 p=2.092e-09 (df=49) day (learning) t(55)= 7.32 F(1)= 53.591 p=1.125e-09 p=2.092e-09 (df=49)
stim (main, window start) t(55)= 0.05 F(1)= 0.002 p=0.9606 p=0.9609 (df=20) stim (main, window start) t(55)= 0.05 F(1)= 0.002 p=0.9606 p=0.9609 (df=20)
interaction 95% CI: [+2.92, +11.00] interaction 95% CI: [+2.922, +11]
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.7)=6.65, p=0.02815 HONEST LME (per-animal random slope, day|rat): interaction F(1,9.7)=6.65, p=0.02815
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) model above is the honest learning-rate test -- DF collapses toward the animal count.)
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.001073, slope diff=+6.96) INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.001073, slope diff=+6.962)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: naive_a2_d0_5 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2, Naive
N = 10 rats, 57 sessions raw day coverage: treat 0..5, control 0..5
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 57
Fixed effects coefficients 4
Random effects coefficients 10
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-81.991 -69.732 46.995 -93.991
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF
{'(Intercept)'} 0.21873 0.03639 6.0107 53
{'day' } 0.030326 0.0094955 3.1937 53
{'stim' } -0.053886 0.062912 -0.85653 53
{'day:stim' } 0.062061 0.016406 3.7829 53
pValue Lower Upper
1.7409e-07 0.14574 0.29172
0.002365 0.01128 0.049372
0.39556 -0.18007 0.072299
0.00039586 0.029156 0.094967
Random effects covariance parameters (95% CIs):
Group: rat (10 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.054581
Lower Upper
0.026757 0.11134
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.096936 0.079103 0.11879
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(53)= 3.78 F(1)= 14.310 p=0.0003959 p=0.0004381 (df=47)
day (learning) t(53)= 3.19 F(1)= 10.200 p=0.002365 p=0.002483 (df=48)
stim (main, window start) t(53)= -0.86 F(1)= 0.734 p=0.3956 p=0.3992 (df=27)
interaction 95% CI: [+0.02916, +0.09497]
HONEST LME (per-animal random slope, day|rat): interaction F(1,11.5)=11.53, p=0.005637
(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 SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.0003959, slope diff=+0.06206)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 36 KiB

After

Width:  |  Height:  |  Size: 36 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 35 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_10 LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_10 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05 observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.006) = 0.076 (small; f: .10 small, .25
16 | 0.34 | 0.36 16 | 0.34 | 0.36
24 | 0.51 | 0.53 24 | 0.51 | 0.53
true stim:log(day) = -8.46 (50% of observed) true stim:log(day) = -8.464 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.03 | 0.07 <- observed 3 | 0.03 | 0.07 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.006) = 0.076 (small; f: .10 small, .25
16 | 0.11 | 0.12 16 | 0.11 | 0.12
24 | 0.13 | 0.14 24 | 0.13 | 0.14
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_10 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,41)=2.375 p(resid)=0.131 p(Satt)=0.1321 (df=36)
honest per-animal random slope (log-day): F(1,9.0)=1.58 p=0.2398
Cohen's f (interaction, partial eta^2=0.023) = 0.153 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=-0.2002, ratSD=0.06161, resSD=0.06563) ---
true stim:log(day) = -0.2002 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.26 <- observed
5 | 0.27 | 0.39
8 | 0.50 | 0.57
12 | 0.78 | 0.82
16 | 0.88 | 0.91
24 | 0.97 | 0.97
true stim:log(day) = -0.1001 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.03 | 0.10 <- observed
5 | 0.07 | 0.12
8 | 0.17 | 0.22
12 | 0.17 | 0.19
16 | 0.33 | 0.39
24 | 0.36 | 0.39
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- naive_a2_d6_10 POWER SIMULATION -- naive_a2_d6_10 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05 observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
ground truth: stim:day=-1.80/day, rat SD=10.99, residual SD=10.55, days=5 ground truth: stim:day=-1.8/day, rat SD=10.99, residual SD=10.55, days=5
true stim:day interaction = -1.80 (100% of observed) true stim:day interaction = -1.8 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.06 | 0.09 <- observed 3 | 0.06 | 0.09 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=-1.80/day, rat SD=10.99, residual SD=10.55, days=5
16 | 0.33 | 0.34 16 | 0.33 | 0.34
24 | 0.45 | 0.49 24 | 0.45 | 0.49
true stim:day interaction = -0.90 (50% of observed) true stim:day interaction = -0.9 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.03 | 0.07 <- observed 3 | 0.03 | 0.07 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=-1.80/day, rat SD=10.99, residual SD=10.55, days=5
16 | 0.11 | 0.11 16 | 0.11 | 0.11
24 | 0.13 | 0.14 24 | 0.13 | 0.14
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- naive_a2_d6_10 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
ground truth: stim:day=-0.02201/day, rat SD=0.06153, residual SD=0.06598, days=5
true stim:day interaction = -0.02201 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.10 | 0.23 <- observed
5 | 0.23 | 0.37
8 | 0.47 | 0.55
12 | 0.75 | 0.78
16 | 0.88 | 0.88
24 | 0.95 | 0.96
true stim:day interaction = -0.01101 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.03 | 0.10 <- observed
5 | 0.07 | 0.11
8 | 0.17 | 0.21
12 | 0.15 | 0.18
16 | 0.33 | 0.38
24 | 0.34 | 0.37
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: naive_a2_d6_10 VARIATION: naive_a2_d6_10 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(41)= -0.76 F(1)= 0.582 p=0.4499 p=0.4505 (df=36) stim x day (interaction) t(41)= -0.76 F(1)= 0.582 p=0.4499 p=0.4505 (df=36)
day (learning) t(41)= 4.38 F(1)= 19.183 p=8.027e-05 p=9.813e-05 (df=36) day (learning) t(41)= 4.38 F(1)= 19.183 p=8.027e-05 p=9.813e-05 (df=36)
stim (main, window start) t(41)= 2.73 F(1)= 7.453 p=0.00929 p=0.01541 (df=15) stim (main, window start) t(41)= 2.73 F(1)= 7.453 p=0.00929 p=0.01541 (df=15)
interaction 95% CI: [-6.57, +2.97] interaction 95% CI: [-6.565, +2.965]
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.0)=0.44, p=0.5223 HONEST LME (per-animal random slope, day|rat): interaction F(1,9.0)=0.44, p=0.5223
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) 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.4499, slope diff=-1.80) INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.4499, slope diff=-1.8)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: naive_a2_d6_10 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2, Naive
N = 9 rats, 45 sessions raw day coverage: treat 6..10, control 6..10
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 45
Fixed effects coefficients 4
Random effects coefficients 9
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-89.86 -79.02 50.93 -101.86
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 0.46152 0.032656 14.133 41 2.355e-17
{'day' } 0.034211 0.008518 4.0163 41 0.00024603
{'stim' } 0.18209 0.056562 3.2193 41 0.0025147
{'day:stim' } -0.022011 0.014754 -1.4919 41 0.14338
Lower Upper
0.39557 0.52747
0.017008 0.051413
0.067859 0.29632
-0.051806 0.0077844
Random effects covariance parameters (95% CIs):
Group: rat (9 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.061534
Lower Upper
0.034776 0.10888
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.06598 0.052372 0.083124
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(41)= -1.49 F(1)= 2.226 p=0.1434 p=0.1444 (df=36)
day (learning) t(41)= 4.02 F(1)= 16.131 p=0.000246 p=0.0002875 (df=36)
stim (main, window start) t(41)= 3.22 F(1)= 10.364 p=0.002515 p=0.005217 (df=16)
interaction 95% CI: [-0.05181, +0.007784]
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.0)=1.61, p=0.2369
(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.1434, slope diff=-0.02201)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -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) % 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')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
C = m.Coefficients; A = anova(m); ci = coefCI(m); 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; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
wst = warning('off', 'all'); wst = warning('off', 'all');
try 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)), ... 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))); 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)); 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)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
if abs(maxT - maxC) > 2 if abs(maxT - maxC) > 2
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
else else
cov = '(equal day coverage over this window)'; cov = '(equal day coverage over this window)';
end 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); bar = repmat('=', 1, 78);
raw = regexprep(evalc('disp(m)'), '</?strong>', ''); 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('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('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('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('stim x day (interaction)', 'day:stim')];
s = [s row('day (learning)', 'day')]; s = [s row('day (learning)', 'day')];
s = [s row('stim (main, window start)', 'stim')]; 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 if rsOk
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
else else
@@ -82,17 +105,18 @@ else
end end
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... 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'])]; ' 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); 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);
Binary file not shown.

Before

Width:  |  Height:  |  Size: 40 KiB

After

Width:  |  Height:  |  Size: 40 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 39 KiB

@@ -1,7 +1,7 @@
============================================================================== ==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_13 LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_13 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
log(day) uses 1-indexed training day (our day 0 = paper "Day 1") log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05 observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
@@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,61)=0.405 p(resid)=0.5267 p(Satt)=0.5268 (df=
honest per-animal random slope (log-day): F(1,8.2)=0.20 p=0.663 honest per-animal random slope (log-day): F(1,8.2)=0.20 p=0.663
Cohen's f (interaction, partial eta^2=0.002) = 0.047 (small; f: .10 small, .25 medium, .40 large) Cohen's f (interaction, partial eta^2=0.002) = 0.047 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+9.35, ratSD=9.78, resSD=11.90) --- --- power simulation (log-day ground truth: stim:logday=+9.352, ratSD=9.782, resSD=11.9) ---
true stim:log(day) = +9.35 (100% of observed) true stim:log(day) = +9.352 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.07 | 0.10 <- observed 3 | 0.07 | 0.10 <- observed
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.002) = 0.047 (small; f: .10 small, .25
16 | 0.30 | 0.32 16 | 0.30 | 0.32
24 | 0.42 | 0.42 24 | 0.42 | 0.42
true stim:log(day) = +4.68 (50% of observed) true stim:log(day) = +4.676 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.05 | 0.09 <- observed 3 | 0.05 | 0.09 <- observed
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.002) = 0.047 (small; f: .10 small, .25
16 | 0.07 | 0.07 16 | 0.07 | 0.07
24 | 0.17 | 0.17 24 | 0.17 | 0.17
Read the per-animal column as the honest power; the LME column matches the Read per-animal as the honest power; LME matches the paper's power code (optimistic).
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,35 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_13 [metric: success RATE]
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,61)=0.040 p(resid)=0.843 p(Satt)=0.843 (df=58)
honest per-animal random slope (log-day): F(1,7.8)=0.00 p=0.9656
Cohen's f (interaction, partial eta^2=0.000) = 0.010 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+0.01851, ratSD=0.05356, resSD=0.07547) ---
true stim:log(day) = +0.01851 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.05 <- observed
5 | 0.09 | 0.08
8 | 0.05 | 0.03
12 | 0.04 | 0.06
16 | 0.07 | 0.08
24 | 0.03 | 0.07
true stim:log(day) = +0.009253 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.07 <- observed
5 | 0.05 | 0.06
8 | 0.07 | 0.04
12 | 0.04 | 0.05
16 | 0.04 | 0.05
24 | 0.07 | 0.07
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
@@ -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). % 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.) % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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('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)]; 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 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 end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA); full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); 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); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part)); cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); 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', ... s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)]; eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects; cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); 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)]; bInt, sRat, sRes)];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -102,19 +102,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -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 % "Lines indicate mean (and SEM) across animals in the anodal (red) and
% control (blue) groups." % 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" % Run: matlab -batch "plotcurve"
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) % (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, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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]; red = [0.85 0.10 0.10];
blue = [0.10 0.30 0.85]; 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]); fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
hold on hold on
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); 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 hold off
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
'Location', 'northwest', 'Box', 'off'); 'Location', 'northwest', 'Box', 'off');
xlabel('training day'); xlabel('training day'); ylabel(ylab);
ylabel('# successes');
title(vname, 'Interpreter', 'none'); title(vname, 'Interpreter', 'none');
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); 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); exportgraphics(fig, outFile, 'Resolution', 150);
close(fig); 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)); subs = unique(D.subject(D.stim == stimVal));
n = numel(subs); n = numel(subs);
X = nan(numel(days), n); X = nan(numel(days), n);
for j = 1:n for j = 1:n
for i = 1:numel(days) for i = 1:numel(days)
r = D.subject == subs(j) & D.day == days(i); 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
end end
M = mean(X, 2, 'omitnan'); M = mean(X, 2, 'omitnan');
@@ -1,11 +1,11 @@
============================================================================== ==============================================================================
POWER SIMULATION -- naive_a2_d6_13 POWER SIMULATION -- naive_a2_d6_13 [metric: # successes (count)]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05 observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
ground truth: stim:day=+1.16/day, rat SD=9.76, residual SD=12.03, days=8 ground truth: stim:day=+1.155/day, rat SD=9.759, residual SD=12.03, days=8
true stim:day interaction = +1.16 (100% of observed) true stim:day interaction = +1.155 (100% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.08 | 0.15 <- observed 3 | 0.08 | 0.15 <- observed
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.16/day, rat SD=9.76, residual SD=12.03, days=8
16 | 0.42 | 0.47 16 | 0.42 | 0.47
24 | 0.58 | 0.59 24 | 0.58 | 0.59
true stim:day interaction = +0.58 (50% of observed) true stim:day interaction = +0.5776 (50% of observed)
N/group | per-animal power | LME power N/group | per-animal power | LME power
------------------------------------------ ------------------------------------------
3 | 0.06 | 0.09 <- observed 3 | 0.06 | 0.09 <- observed
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.16/day, rat SD=9.76, residual SD=12.03, days=8
16 | 0.12 | 0.12 16 | 0.12 | 0.12
24 | 0.21 | 0.23 24 | 0.21 | 0.23
Read the per-animal column as the honest power. At the observed N this study Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
@@ -0,0 +1,28 @@
==============================================================================
POWER SIMULATION -- naive_a2_d6_13 [metric: success RATE]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
ground truth: stim:day=+0.003468/day, rat SD=0.05337, residual SD=0.07609, days=8
true stim:day interaction = +0.003468 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.05 <- observed
5 | 0.12 | 0.13
8 | 0.09 | 0.11
12 | 0.09 | 0.08
16 | 0.15 | 0.17
24 | 0.14 | 0.15
true stim:day interaction = +0.001734 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.07 <- observed
5 | 0.04 | 0.07
8 | 0.07 | 0.07
12 | 0.05 | 0.05
16 | 0.06 | 0.07
24 | 0.08 | 0.08
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
@@ -1,45 +1,49 @@
% Variation power simulation -- Monte-Carlo power for the paper's stim x day % 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 % 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) % 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.) % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
here = fileparts(mfilename('fullpath')); here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', ''); vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); 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)'; 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); 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'}); 'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1))); nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0))); nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78); 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)]; 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 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 end
lme = fitlme(tbl0, FORMULA); 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); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
days = (0:max(tbl0.day))'; 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))]; bInt, sRat, sRes, numel(days))];
for eMul = EFFMULS for eMul = EFFMULS
bI = bInt * eMul; 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 for N = NS
sigPA = 0; sigL = 0; sigPA = 0; sigL = 0;
for r = 1:NREP for r = 1:NREP
@@ -70,20 +73,18 @@ for eMul = EFFMULS
end end
star = ''; star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end 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
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); warning(warnState);
end
% ---------------------------------------------------------------- helpers % ---------------------------------------------------------------- helpers
function localFinish(s, here) function localFinish(s, here, suffix)
fprintf('%s', s); 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); fprintf(fid, '%s', s); fclose(fid);
end end
@@ -1,5 +1,5 @@
============================================================================== ==============================================================================
VARIATION: naive_a2_d6_13 VARIATION: naive_a2_d6_13 [metric: success COUNT]
============================================================================== ==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT) model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
day = training day within window (0 = first analyzed day) day = training day within window (0 = first analyzed day)
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
stim x day (interaction) t(61)= 0.78 F(1)= 0.605 p=0.4398 p=0.4399 (df=58) stim x day (interaction) t(61)= 0.78 F(1)= 0.605 p=0.4398 p=0.4399 (df=58)
day (learning) t(61)= 3.39 F(1)= 11.485 p=0.001235 p=0.001274 (df=57) day (learning) t(61)= 3.39 F(1)= 11.485 p=0.001235 p=0.001274 (df=57)
stim (main, window start) t(61)= 2.46 F(1)= 6.072 p=0.01657 p=0.02473 (df=17) stim (main, window start) t(61)= 2.46 F(1)= 6.072 p=0.01657 p=0.02473 (df=17)
interaction 95% CI: [-1.82, +4.13] interaction 95% CI: [-1.815, +4.126]
HONEST LME (per-animal random slope, day|rat): interaction F(1,7.7)=0.36, p=0.5635 HONEST LME (per-animal random slope, day|rat): interaction F(1,7.7)=0.36, p=0.5635
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope (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.) 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.4398, slope diff=+1.16) INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.4398, slope diff=+1.155)
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
@@ -0,0 +1,68 @@
==============================================================================
VARIATION: naive_a2_d6_13 [metric: success RATE (success/attempts)]
==============================================================================
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
day = training day within window (0 = first analyzed day)
treatment (stim=1): Electrode-Box-B2
control (stim=0): Electrode-Box-A2, Naive
N = 9 rats, 65 sessions raw day coverage: treat 6..13, control 6..13
(equal day coverage over this window)
==============================================================================
FULL MODEL SUMMARY -- fitlme
==============================================================================
Linear mixed-effects model fit by ML
Model information:
Number of observations 65
Fixed effects coefficients 4
Random effects coefficients 9
Covariance parameters 2
Formula:
behavior ~ 1 + day*stim + (1 | rat)
Model fit statistics:
AIC BIC LogLikelihood Deviance
-124.85 -111.8 68.424 -136.85
Fixed effects coefficients (95% CIs):
Name Estimate SE tStat DF pValue
{'(Intercept)'} 0.49381 0.029854 16.541 61 3.1921e-24
{'day' } 0.012224 0.0052666 2.321 61 0.023647
{'stim' } 0.14414 0.051782 2.7835 61 0.007149
{'day:stim' } 0.0034682 0.0093748 0.36995 61 0.7127
Lower Upper
0.43411 0.5535
0.0016924 0.022755
0.040593 0.24768
-0.015278 0.022214
Random effects covariance parameters (95% CIs):
Group: rat (9 Levels)
Name1 Name2 Type Estimate
{'(Intercept)'} {'(Intercept)'} {'std'} 0.053368
Lower Upper
0.029753 0.095724
Group: Error
Name Estimate Lower Upper
{'Res Std'} 0.076088 0.063265 0.091511
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
----------------------------------------------------------------------------
stim x day (interaction) t(61)= 0.37 F(1)= 0.137 p=0.7127 p=0.7128 (df=58)
day (learning) t(61)= 2.32 F(1)= 5.387 p=0.02365 p=0.02384 (df=58)
stim (main, window start) t(61)= 2.78 F(1)= 7.748 p=0.007149 p=0.01177 (df=19)
interaction 95% CI: [-0.01528, +0.02221]
HONEST LME (per-animal random slope, day|rat): interaction F(1,7.4)=0.04, p=0.8393
(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.7127, slope diff=+0.003468)
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.