analysis(matlab): add log(day) model + Cohen's f + power to every variation

Add variation_logpower.m (self-contained) + make_variation_logpower.m,
dropping logpowersim.m + logpower_result.txt into all 28 variation folders.
Matches the paper's power code: models behavior ~ stim + log(day) +
stim:log(day) + (1|rat) (log(day+1), since our day 0 = paper Day 1), reports
Cohen's f (partial eta^2 of the interaction) and the interaction under
residual/Satterthwaite/honest random-slope DF, then runs the Monte-Carlo power
sim on the log-day ground truth (per-animal cluster-honest + LME power).

Notable: under log(day) the accumulating divergence is captured more sharply,
so several full-window scenarios reach honest significance that were n.s. under
raw day (e.g. right_only_d0_13 honest p=0.003, unmerge_d0_13 0.021,
unmerge_d0_10 0.036); Cohen's f is small-medium (~0.10-0.34).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
2026-07-23 14:37:50 -04:00
parent bb5e2d88c0
commit 378ad6da05
58 changed files with 5328 additions and 0 deletions
+28
View File
@@ -0,0 +1,28 @@
function make_variation_logpower()
%MAKE_VARIATION_LOGPOWER Add the log-day model + Cohen's f + power to each variation.
% MAKE_VARIATION_LOGPOWER() copies variation_logpower.m into every
% variations/<name>/ folder with a data.csv (as logpowersim.m) and runs it,
% producing logpower_result.txt. Re-runnable; picks up new folders.
thisDir = fileparts(mfilename('fullpath'));
template = fullfile(thisDir, 'variation_logpower.m');
root = fullfile(thisDir, 'variations');
d = dir(root);
n = 0;
for i = 1:numel(d)
if ~d(i).isdir || ismember(d(i).name, {'.', '..'}); continue; end
folder = fullfile(root, d(i).name);
if ~exist(fullfile(folder, 'data.csv'), 'file'); continue; end
copyfile(template, fullfile(folder, 'logpowersim.m'));
fprintf('\n### %s ###\n', d(i).name);
localRun(fullfile(folder, 'logpowersim.m'));
n = n + 1;
end
fprintf('\nAdded logpowersim.m + logpower_result.txt to %d variation folders.\n', n);
end
function localRun(scriptPath)
run(scriptPath);
end
+148
View File
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,68)=3.821 p(resid)=0.05474 p(Satt)=0.05451 (df=72)
honest per-animal random slope (log-day): F(1,14.6)=3.22 p=0.09348
Cohen's f (interaction, partial eta^2=0.009) = 0.097 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+8.69, ratSD=0.00, resSD=13.42) ---
true stim:log(day) = +8.69 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.23 | 0.47 <- observed
5 | 0.57 | 0.72
8 | 0.83 | 0.87
12 | 0.98 | 0.99
16 | 0.99 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +4.34 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.12 <- observed
5 | 0.11 | 0.19
8 | 0.30 | 0.38
12 | 0.34 | 0.40
16 | 0.62 | 0.64
24 | 0.75 | 0.78
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,38)=4.320 p(resid)=0.04447 p(Satt)=0.04382 (df=42)
honest per-animal random slope (log-day): F(1,10.2)=3.63 p=0.0853
Cohen's f (interaction, partial eta^2=0.032) = 0.181 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+15.16, ratSD=0.00, resSD=14.15) ---
true stim:log(day) = +15.16 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.22 | 0.51 <- observed
5 | 0.55 | 0.68
8 | 0.81 | 0.85
12 | 0.95 | 0.97
16 | 1.00 | 0.99
24 | 1.00 | 1.00
true stim:log(day) = +7.58 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.08 | 0.23 <- observed
5 | 0.17 | 0.23
8 | 0.24 | 0.27
12 | 0.42 | 0.47
16 | 0.55 | 0.58
24 | 0.85 | 0.85
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,26)=0.252 p(resid)=0.6199 p(Satt)=0.6202 (df=24)
honest per-animal random slope (log-day): F(1,6.0)=0.19 p=0.6783
Cohen's f (interaction, partial eta^2=0.004) = 0.065 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+9.52, ratSD=6.32, resSD=8.30) ---
true stim:log(day) = +9.52 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.02 | 0.06 <- observed
5 | 0.09 | 0.13
8 | 0.12 | 0.12
12 | 0.11 | 0.17
16 | 0.14 | 0.16
24 | 0.32 | 0.34
true stim:log(day) = +4.76 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.06 | 0.10 <- observed
5 | 0.07 | 0.07
8 | 0.05 | 0.07
12 | 0.09 | 0.10
16 | 0.09 | 0.12
24 | 0.12 | 0.15
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,68)=5.859 p(resid)=0.01817 p(Satt)=0.0182 (df=67)
honest per-animal random slope (log-day): F(1,21.8)=5.17 p=0.03322
Cohen's f (interaction, partial eta^2=0.013) = 0.115 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+10.40, ratSD=5.52, resSD=12.53) ---
true stim:log(day) = +10.40 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.32 | 0.62 <- observed
5 | 0.78 | 0.87
8 | 0.93 | 0.97
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +5.20 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.10 | 0.20 <- observed
5 | 0.17 | 0.31
8 | 0.45 | 0.53
12 | 0.57 | 0.60
16 | 0.86 | 0.85
24 | 0.93 | 0.93
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,38)=3.165 p(resid)=0.08322 p(Satt)=0.08246 (df=42)
honest per-animal random slope (log-day): F(1,9.4)=2.59 p=0.1407
Cohen's f (interaction, partial eta^2=0.026) = 0.163 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+13.68, ratSD=0.00, resSD=14.91) ---
true stim:log(day) = +13.68 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.18 | 0.42 <- observed
5 | 0.43 | 0.59
8 | 0.68 | 0.75
12 | 0.81 | 0.88
16 | 0.96 | 0.97
24 | 0.98 | 0.99
true stim:log(day) = +6.84 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.14 <- observed
5 | 0.12 | 0.19
8 | 0.17 | 0.19
12 | 0.31 | 0.39
16 | 0.43 | 0.49
24 | 0.70 | 0.75
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,26)=0.850 p(resid)=0.365 p(Satt)=0.3657 (df=24)
honest per-animal random slope (log-day): F(1,6.0)=0.67 p=0.4431
Cohen's f (interaction, partial eta^2=0.014) = 0.118 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+18.32, ratSD=6.21, resSD=8.20) ---
true stim:log(day) = +18.32 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.06 | 0.17
5 | 0.17 | 0.26
8 | 0.27 | 0.36
12 | 0.41 | 0.42
16 | 0.57 | 0.63
24 | 0.75 | 0.78
true stim:log(day) = +9.16 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.11
5 | 0.10 | 0.12
8 | 0.07 | 0.13
12 | 0.16 | 0.15
16 | 0.20 | 0.20
24 | 0.31 | 0.31
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- matched_current_d0_3
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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.043 p(resid)=0.01524 p(Satt)=0.0139 (df=24)
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)
--- power simulation (log-day ground truth: stim:logday=+20.79, ratSD=0.00, resSD=9.99) ---
true stim:log(day) = +20.79 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.39 | 0.77 <- observed
5 | 0.84 | 0.93
8 | 0.97 | 0.99
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +10.40 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.11 | 0.23 <- observed
5 | 0.33 | 0.47
8 | 0.49 | 0.59
12 | 0.80 | 0.79
16 | 0.79 | 0.84
24 | 0.95 | 0.97
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,100)=4.633 p(resid)=0.03378 p(Satt)=0.03391 (df=95)
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)
--- power simulation (log-day ground truth: stim:logday=+9.59, ratSD=8.17, resSD=14.84) ---
true stim:log(day) = +9.59 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.23 | 0.44 <- observed
5 | 0.57 | 0.68
8 | 0.83 | 0.86
12 | 0.97 | 0.98
16 | 0.99 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +4.79 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.10 <- observed
5 | 0.11 | 0.17
8 | 0.30 | 0.37
12 | 0.34 | 0.38
16 | 0.62 | 0.64
24 | 0.75 | 0.78
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,120)=6.786 p(resid)=0.01035 p(Satt)=0.01039 (df=116)
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)
--- power simulation (log-day ground truth: stim:logday=+9.94, ratSD=8.36, resSD=14.40) ---
true stim:log(day) = +9.94 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.27 | 0.71 <- observed
5 | 0.68 | 0.82
8 | 0.97 | 0.97
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +4.97 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.25 <- observed
5 | 0.27 | 0.31
8 | 0.51 | 0.57
12 | 0.62 | 0.68
16 | 0.80 | 0.85
24 | 0.92 | 0.93
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,55)=9.364 p(resid)=0.003417 p(Satt)=0.003585 (df=49)
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)
--- power simulation (log-day ground truth: stim:logday=+19.33, ratSD=9.47, resSD=13.37) ---
true stim:log(day) = +19.33 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.38 | 0.73 <- observed
5 | 0.76 | 0.92
8 | 0.98 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +9.67 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.28 <- observed
5 | 0.31 | 0.40
8 | 0.46 | 0.53
12 | 0.67 | 0.68
16 | 0.83 | 0.84
24 | 0.95 | 0.96
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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)=0.662 p(resid)=0.4205 p(Satt)=0.4211 (df=36)
honest per-animal random slope (log-day): F(1,9.0)=0.46 p=0.5139
Cohen's f (interaction, partial eta^2=0.006) = 0.076 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=-16.93, ratSD=10.99, resSD=10.51) ---
true stim:log(day) = -16.93 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.12 <- observed
5 | 0.07 | 0.07
8 | 0.20 | 0.27
12 | 0.30 | 0.34
16 | 0.34 | 0.36
24 | 0.51 | 0.53
true stim:log(day) = -8.46 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.03 | 0.07 <- observed
5 | 0.07 | 0.07
8 | 0.08 | 0.11
12 | 0.07 | 0.08
16 | 0.11 | 0.12
24 | 0.13 | 0.14
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_a2_d6_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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.405 p(resid)=0.5267 p(Satt)=0.5268 (df=58)
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)
--- power simulation (log-day ground truth: stim:logday=+9.35, ratSD=9.78, resSD=11.90) ---
true stim:log(day) = +9.35 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.10 <- observed
5 | 0.19 | 0.19
8 | 0.17 | 0.19
12 | 0.25 | 0.25
16 | 0.30 | 0.32
24 | 0.42 | 0.42
true stim:log(day) = +4.68 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.09 <- observed
5 | 0.05 | 0.07
8 | 0.07 | 0.10
12 | 0.11 | 0.09
16 | 0.07 | 0.07
24 | 0.17 | 0.17
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_boxa_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,111)=3.705 p(resid)=0.05681 p(Satt)=0.05696 (df=105)
honest per-animal random slope (log-day): F(1,8.2)=2.03 p=0.1913
Cohen's f (interaction, partial eta^2=0.007) = 0.086 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+8.43, ratSD=8.13, resSD=14.95) ---
true stim:log(day) = +8.43 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.17 | 0.35 <- observed
5 | 0.44 | 0.56
8 | 0.68 | 0.69 <- observed
12 | 0.91 | 0.92
16 | 0.95 | 0.96
24 | 0.99 | 1.00
true stim:log(day) = +4.21 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.07 <- observed
5 | 0.08 | 0.14
8 | 0.23 | 0.29 <- observed
12 | 0.28 | 0.28
16 | 0.55 | 0.53
24 | 0.63 | 0.62
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_boxa_d0_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,134)=4.865 p(resid)=0.02911 p(Satt)=0.02918 (df=129)
honest per-animal random slope (log-day): F(1,7.9)=3.22 p=0.1109
Cohen's f (interaction, partial eta^2=0.007) = 0.083 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+8.27, ratSD=8.74, resSD=14.51) ---
true stim:log(day) = +8.27 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.21 | 0.55 <- observed
5 | 0.51 | 0.65
8 | 0.88 | 0.90 <- observed
12 | 0.97 | 1.00
16 | 0.98 | 0.98
24 | 1.00 | 1.00
true stim:log(day) = +4.13 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.08 | 0.17 <- observed
5 | 0.19 | 0.27
8 | 0.42 | 0.42 <- observed
12 | 0.47 | 0.51
16 | 0.64 | 0.68
24 | 0.76 | 0.78
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_boxa_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,61)=8.016 p(resid)=0.006273 p(Satt)=0.006505 (df=54)
honest per-animal random slope (log-day): F(1,15.5)=6.72 p=0.01995
Cohen's f (interaction, partial eta^2=0.048) = 0.223 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+18.30, ratSD=8.82, resSD=13.99) ---
true stim:log(day) = +18.30 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.33 | 0.66 <- observed
5 | 0.73 | 0.82
8 | 0.95 | 0.98 <- observed
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +9.15 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.11 | 0.24 <- observed
5 | 0.23 | 0.31
8 | 0.37 | 0.43 <- observed
12 | 0.56 | 0.63
16 | 0.74 | 0.78
24 | 0.93 | 0.94
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_boxa_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,46)=0.613 p(resid)=0.4376 p(Satt)=0.4382 (df=40)
honest per-animal random slope (log-day): F(1,10.0)=0.43 p=0.5247
Cohen's f (interaction, partial eta^2=0.005) = 0.069 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=-15.28, ratSD=11.24, resSD=10.10) ---
true stim:log(day) = -15.28 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.06 | 0.08 <- observed
5 | 0.07 | 0.07
8 | 0.18 | 0.21
12 | 0.28 | 0.30
16 | 0.30 | 0.35
24 | 0.45 | 0.49
true stim:log(day) = -7.64 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.03 | 0.09 <- observed
5 | 0.07 | 0.07
8 | 0.07 | 0.11
12 | 0.07 | 0.08
16 | 0.10 | 0.11
24 | 0.13 | 0.14
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- naive_boxa_d6_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,69)=0.239 p(resid)=0.6263 p(Satt)=0.6264 (df=65)
honest per-animal random slope (log-day): F(1,9.1)=0.07 p=0.8013
Cohen's f (interaction, partial eta^2=0.001) = 0.034 (small; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+6.70, ratSD=10.78, resSD=11.38) ---
true stim:log(day) = +6.70 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.07 <- observed
5 | 0.13 | 0.15
8 | 0.11 | 0.17
12 | 0.15 | 0.12
16 | 0.23 | 0.25
24 | 0.23 | 0.24
true stim:log(day) = +3.35 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.08 <- observed
5 | 0.04 | 0.07
8 | 0.06 | 0.10
12 | 0.05 | 0.05
16 | 0.06 | 0.07
24 | 0.12 | 0.11
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- prev_f_full
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=12, control n=12 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,227)=8.356 p(resid)=0.004218 p(Satt)=0.004253 (df=208)
honest per-animal random slope (log-day): F(1,24.3)=3.32 p=0.08079
Cohen's f (interaction, partial eta^2=0.012) = 0.112 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+2.39, ratSD=4.33, resSD=4.34) ---
true stim:log(day) = +2.39 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.15 | 0.28
5 | 0.33 | 0.46
8 | 0.61 | 0.65
12 | 0.82 | 0.83 <- observed
16 | 0.92 | 0.94
24 | 0.99 | 0.98
true stim:log(day) = +1.19 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.05 | 0.09
5 | 0.12 | 0.17
8 | 0.20 | 0.26
12 | 0.33 | 0.32 <- observed
16 | 0.46 | 0.49
24 | 0.46 | 0.49
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,68)=8.420 p(resid)=0.004998 p(Satt)=0.005011 (df=67)
honest per-animal random slope (log-day): F(1,22.8)=7.55 p=0.01151
Cohen's f (interaction, partial eta^2=0.016) = 0.126 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+11.36, ratSD=4.91, resSD=11.42) ---
true stim:log(day) = +11.36 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.43 | 0.79 <- observed
5 | 0.88 | 0.95
8 | 0.98 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +5.68 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.13 | 0.29 <- observed
5 | 0.30 | 0.42
8 | 0.65 | 0.62
12 | 0.78 | 0.82
16 | 0.93 | 0.93
24 | 0.99 | 0.99
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,80)=13.182 p(resid)=0.0004969 p(Satt)=0.0004951 (df=81)
honest per-animal random slope (log-day): F(1,18.9)=11.35 p=0.003237
Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+12.18, ratSD=5.40, resSD=10.82) ---
true stim:log(day) = +12.18 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.68 | 0.97 <- observed
5 | 1.00 | 1.00
8 | 1.00 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +6.09 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.25 | 0.47 <- observed
5 | 0.53 | 0.69
8 | 0.84 | 0.87
12 | 0.93 | 0.96
16 | 1.00 | 1.00
24 | 1.00 | 1.00
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,38)=5.667 p(resid)=0.02241 p(Satt)=0.0219 (df=42)
honest per-animal random slope (log-day): F(1,9.0)=4.64 p=0.05966
Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+16.20, ratSD=0.00, resSD=13.20) ---
true stim:log(day) = +16.20 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.33 | 0.65 <- observed
5 | 0.66 | 0.78
8 | 0.93 | 0.97
12 | 0.99 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +8.10 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.09 | 0.24 <- observed
5 | 0.21 | 0.29
8 | 0.29 | 0.38
12 | 0.50 | 0.60
16 | 0.69 | 0.69
24 | 0.92 | 0.92
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,26)=1.327 p(resid)=0.2599 p(Satt)=0.2608 (df=24)
honest per-animal random slope (log-day): F(1,6.0)=0.96 p=0.3641
Cohen's f (interaction, partial eta^2=0.021) = 0.145 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+23.04, ratSD=5.69, resSD=8.25) ---
true stim:log(day) = +23.04 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.28
5 | 0.28 | 0.38
8 | 0.44 | 0.50
12 | 0.61 | 0.61
16 | 0.76 | 0.78
24 | 0.97 | 0.97
true stim:log(day) = +11.52 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.13
5 | 0.11 | 0.14
8 | 0.10 | 0.15
12 | 0.19 | 0.23
16 | 0.25 | 0.25
24 | 0.45 | 0.52
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,38)=5.774 p(resid)=0.02125 p(Satt)=0.02133 (df=37)
honest per-animal random slope (log-day): F(1,0.0)=3.19 p=NaN
Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+28.42, ratSD=6.48, resSD=7.39) ---
true stim:log(day) = +28.42 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.57 | 0.86
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:log(day) = +14.21 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.17 | 0.30
5 | 0.37 | 0.55
8 | 0.62 | 0.68
12 | 0.81 | 0.85
16 | 0.90 | 0.92
24 | 1.00 | 1.00
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d0_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,57)=6.230 p(resid)=0.01548 p(Satt)=0.01528 (df=61)
honest per-animal random slope (log-day): F(1,15.0)=5.28 p=0.0363
Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+11.60, ratSD=0.00, resSD=12.93) ---
true stim:log(day) = +11.60 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.34 | 0.69 <- observed
5 | 0.82 | 0.92
8 | 0.97 | 0.97
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +5.80 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.10 | 0.25 <- observed
5 | 0.26 | 0.40
8 | 0.53 | 0.58
12 | 0.61 | 0.72
16 | 0.90 | 0.90
24 | 0.97 | 0.97
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d0_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,66)=8.712 p(resid)=0.004373 p(Satt)=0.004357 (df=67)
honest per-animal random slope (log-day): F(1,12.9)=6.96 p=0.02062
Cohen's f (interaction, partial eta^2=0.015) = 0.124 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+11.06, ratSD=5.68, resSD=11.36) ---
true stim:log(day) = +11.06 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.53 | 0.91 <- observed
5 | 0.99 | 1.00
8 | 1.00 | 1.00
12 | 1.00 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +5.53 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.18 | 0.38 <- observed
5 | 0.41 | 0.54
8 | 0.78 | 0.80
12 | 0.86 | 0.91
16 | 0.97 | 0.98
24 | 0.98 | 0.98
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d0_5
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
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,32)=4.849 p(resid)=0.03498 p(Satt)=0.03414 (df=36)
honest per-animal random slope (log-day): F(1,7.2)=3.65 p=0.09636
Cohen's f (interaction, partial eta^2=0.039) = 0.200 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+16.48, ratSD=0.00, resSD=13.58) ---
true stim:log(day) = +16.48 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.33 | 0.64 <- observed
5 | 0.66 | 0.78
8 | 0.93 | 0.97
12 | 0.99 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +8.24 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.09 | 0.24 <- observed
5 | 0.21 | 0.29
8 | 0.28 | 0.38
12 | 0.47 | 0.59
16 | 0.69 | 0.69
24 | 0.90 | 0.92
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d6_10
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=2 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,21)=0.581 p(resid)=0.4545 p(Satt)=0.4549 (df=20)
honest per-animal random slope (log-day): F(1,5.0)=0.44 p=0.5382
Cohen's f (interaction, partial eta^2=0.011) = 0.106 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+16.98, ratSD=6.08, resSD=8.72) ---
true stim:log(day) = +16.98 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.06 | 0.11 <- observed
5 | 0.17 | 0.23
8 | 0.22 | 0.31
12 | 0.33 | 0.33
16 | 0.47 | 0.50
24 | 0.63 | 0.66
true stim:log(day) = +8.49 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.07 | 0.11 <- observed
5 | 0.09 | 0.11
8 | 0.06 | 0.09
12 | 0.12 | 0.14
16 | 0.17 | 0.17
24 | 0.26 | 0.26
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end
@@ -0,0 +1,36 @@
==============================================================================
LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d6_13
==============================================================================
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
observed groups: stim n=3, control n=2 nrep=120, alpha=0.05
--- fitted on real data ---
stim x log(day) interaction: F(1,30)=3.007 p(resid)=0.09315 p(Satt)=0.0931 (df=30)
honest per-animal random slope (log-day): F(1,12.8)=2.47 p=0.1406
Cohen's f (interaction, partial eta^2=0.025) = 0.161 (small-medium; f: .10 small, .25 medium, .40 large)
--- power simulation (log-day ground truth: stim:logday=+22.97, ratSD=7.14, resSD=7.81) ---
true stim:log(day) = +22.97 (100% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.47 | 0.64 <- observed
5 | 0.79 | 0.90
8 | 0.96 | 0.99
12 | 0.99 | 1.00
16 | 1.00 | 1.00
24 | 1.00 | 1.00
true stim:log(day) = +11.49 (50% of observed)
N/group | per-animal power | LME power
------------------------------------------
3 | 0.12 | 0.22 <- observed
5 | 0.26 | 0.36
8 | 0.41 | 0.47
12 | 0.53 | 0.58
16 | 0.71 | 0.74
24 | 0.90 | 0.93
Read the per-animal column as the honest power; the LME column matches the
paper's power code (anova interaction p, observation-level DF) and is optimistic.
@@ -0,0 +1,148 @@
% Variation log-day analysis + Cohen's f + power simulation.
%
% The paper's power code models behavior against LOG training day, not raw day:
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
% "Day 1"), so log(day + 1) reproduces their transform exactly.
%
% This script (a) refits that log-day model on data.csv, (b) reports the
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
% scoring per-animal (cluster-honest) and LME power across N.
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
here = fileparts(mfilename('fullpath'));
if isempty(here); here = pwd; end
vname = regexprep(here, '.*[/\\]', '');
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
NS = [3 5 8 12 16 24];
EFFMULS = [1 0.5];
NREP = 120;
warnState = warning('off', 'all');
rng(1);
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
nStim = numel(unique(D.subject(D.stim == 1)));
nCtrl = numel(unique(D.subject(D.stim == 0)));
bar = repmat('=', 1, 78);
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
localFinish(s, here); warning(warnState); return
end
% ---- fitted model on the real data ----
full = fitlme(tbl0, FORMULA);
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
cohenf = sqrt(eta2part / (1 - eta2part));
% honest per-animal random-slope interaction
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
try
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
Ar = anova(mr, 'DFMethod', 'satterthwaite'); ri = strcmp(Ar.Term, 'day:stim');
rsF = Ar.FStat(ri); rsDf = Ar.DF2(ri); rsP = Ar.pValue(ri); rsOk = true;
catch
end
if cohenf < 0.10; mag = 'small';
elseif cohenf < 0.25; mag = 'small-medium';
elseif cohenf < 0.40; mag = 'medium';
else; mag = 'large'; end
s = [s sprintf('\n--- fitted on real data ---\n')];
s = [s sprintf('stim x log(day) interaction: F(1,%d)=%.3f p(resid)=%.4g p(Satt)=%.4g (df=%.0f)\n', ...
An.DF2(ii), An.FStat(ii), An.pValue(ii), Asatt.pValue(is), Asatt.DF2(is))];
if rsOk
s = [s sprintf(' honest per-animal random slope (log-day): F(1,%.1f)=%.2f p=%.4g\n', rsDf, rsF, rsP)];
else
s = [s sprintf(' honest per-animal random slope (log-day): did not converge\n')];
end
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
eta2part, cohenf, mag)];
% ---- power simulation on the log-day ground truth ----
cn = full.CoefficientNames; be = full.fixedEffects;
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
days = unique(tbl0.day); % the log(day) grid
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
bInt, sRat, sRes)];
for eMul = EFFMULS
bI = bInt * eMul;
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
for N = NS
sigPA = 0; sigL = 0;
for r = 1:NREP
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
sigPA = sigPA + (localPAp(tb) < 0.05);
try
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
catch
end
end
star = '';
if N == nStim || N == nCtrl; star = ' <- observed'; end
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
end
end
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
localFinish(s, here);
warning(warnState);
% ---------------------------------------------------------------- helpers
function localFinish(s, here)
fprintf('%s', s);
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
fprintf(fid, '%s', s); fclose(fid);
end
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
nd = numel(days); rows = 2 * N * nd;
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
k = 0;
for g = 0:1
for sIdx = 1:N
re = sRat * randn;
rid = sprintf('g%d_r%d', g, sIdx);
for d = 1:nd
k = k + 1;
rat(k) = rid; day(k) = days(d); stim(k) = g;
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
end
end
end
tbl = table(categorical(rat), day, stim, behavior, ...
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
end
function p = localPAp(tbl)
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
for i = 1:numel(rats)
r = tbl.rat == rats(i);
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
gr(i) = tbl.stim(find(r, 1));
end
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
end