analysis(matlab): rate + count outputs, variation batch 1
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 41 KiB After Width: | Height: | Size: 41 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 42 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_10
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,68)=3.821 p(resid)=0.05474 p(Satt)=0.05451 (d
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+8.69, ratSD=2.979e-15, resSD=13.42) ---
|
||||
|
||||
true stim:log(day) = +8.69 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.009) = 0.097 (small; f: .10 small, .25
|
||||
16 | 0.99 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +4.34 (50% of observed)
|
||||
true stim:log(day) = +4.345 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.12 <- observed
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.009) = 0.097 (small; f: .10 small, .25
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,68)=5.247 p(resid)=0.02509 p(Satt)=0.02516 (df=66)
|
||||
honest per-animal random slope (log-day): F(1,25.2)=4.88 p=0.0364
|
||||
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=+0.06695, ratSD=0.03723, resSD=0.0875) ---
|
||||
|
||||
true stim:log(day) = +0.06695 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.25 | 0.57 <- observed
|
||||
5 | 0.71 | 0.84
|
||||
8 | 0.92 | 0.93
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +0.03347 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.09 | 0.16 <- observed
|
||||
5 | 0.14 | 0.27
|
||||
8 | 0.41 | 0.46
|
||||
12 | 0.52 | 0.55
|
||||
16 | 0.78 | 0.77
|
||||
24 | 0.88 | 0.87
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d0_10
|
||||
POWER SIMULATION -- boxa_a2_d0_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+1.20/day, rat SD=0.00, residual SD=13.47, days=11
|
||||
ground truth: stim:day=+1.198/day, rat SD=0, residual SD=13.47, days=11
|
||||
|
||||
true stim:day interaction = +1.20 (100% of observed)
|
||||
true stim:day interaction = +1.198 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.17 <- observed
|
||||
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.20/day, rat SD=0.00, residual SD=13.47, days=11
|
||||
16 | 0.72 | 0.72
|
||||
24 | 0.88 | 0.88
|
||||
|
||||
true stim:day interaction = +0.60 (50% of observed)
|
||||
true stim:day interaction = +0.5992 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.08 <- observed
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.20/day, rat SD=0.00, residual SD=13.47, days=11
|
||||
16 | 0.29 | 0.30
|
||||
24 | 0.31 | 0.33
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d0_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+0.007205/day, rat SD=0.04049, residual SD=0.08883, days=11
|
||||
|
||||
true stim:day interaction = +0.007205 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.16 <- observed
|
||||
5 | 0.18 | 0.19
|
||||
8 | 0.28 | 0.33
|
||||
12 | 0.53 | 0.51
|
||||
16 | 0.62 | 0.62
|
||||
24 | 0.84 | 0.85
|
||||
|
||||
true stim:day interaction = +0.003602 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.05 | 0.07 <- observed
|
||||
5 | 0.04 | 0.06
|
||||
8 | 0.09 | 0.10
|
||||
12 | 0.15 | 0.16
|
||||
16 | 0.25 | 0.26
|
||||
24 | 0.23 | 0.26
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d0_10
|
||||
VARIATION: boxa_a2_d0_10 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(68)= 1.18 F(1)= 1.397 p=0.2414 p=0.2412 (df=72)
|
||||
day (learning) t(68)= 12.04 F(1)=144.862 p=1.64e-18 p=6.575e-19 (df=72)
|
||||
stim (main, window start) t(68)= 1.16 F(1)= 1.336 p=0.2517 p=0.2515 (df=72)
|
||||
interaction 95% CI: [-0.83, +3.22]
|
||||
interaction 95% CI: [-0.8251, +3.222]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,14.0)=0.68, p=0.4245
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.2414, slope diff=+1.20)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.2414, slope diff=+1.198)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d0_10 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A, Electrode-Box-A2
|
||||
N = 7 rats, 72 sessions raw day coverage: treat 0..10, control 0..10
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 72
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 7
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-124.36 -110.7 68.18 -136.36
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF pValue
|
||||
{'(Intercept)'} 0.21277 0.032537 6.5394 68 9.4885e-09
|
||||
{'day' } 0.04572 0.0046771 9.7753 68 1.3713e-14
|
||||
{'stim' } 0.042139 0.049417 0.85273 68 0.3968
|
||||
{'day:stim' } 0.0072048 0.0067664 1.0648 68 0.29074
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.14785 0.2777
|
||||
0.036387 0.055053
|
||||
-0.056471 0.14075
|
||||
-0.0062973 0.020707
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (7 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.040493
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.018248 0.089857
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.088826 0.074767 0.10553
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(68)= 1.06 F(1)= 1.134 p=0.2907 p=0.2908 (df=67)
|
||||
day (learning) t(68)= 9.78 F(1)= 95.556 p=1.371e-14 p=1.275e-14 (df=68)
|
||||
stim (main, window start) t(68)= 0.85 F(1)= 0.727 p=0.3968 p=0.4044 (df=19)
|
||||
interaction 95% CI: [-0.006297, +0.02071]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,29.1)=0.97, p=0.3332
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.2907, slope diff=+0.007205)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 36 KiB After Width: | Height: | Size: 35 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 37 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_5
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_5 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,38)=4.320 p(resid)=0.04447 p(Satt)=0.04382 (d
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+15.16, ratSD=0, resSD=14.15) ---
|
||||
|
||||
true stim:log(day) = +15.16 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.032) = 0.181 (small-medium; f: .10 smal
|
||||
16 | 1.00 | 0.99
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +7.58 (50% of observed)
|
||||
true stim:log(day) = +7.581 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.08 | 0.23 <- observed
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.032) = 0.181 (small-medium; f: .10 smal
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d0_5 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,38)=8.750 p(resid)=0.005301 p(Satt)=0.005517 (df=35)
|
||||
honest per-animal random slope (log-day): F(1,12.5)=7.35 p=0.01836
|
||||
Cohen's f (interaction, partial eta^2=0.099) = 0.331 (medium; f: .10 small, .25 medium, .40 large)
|
||||
|
||||
--- power simulation (log-day ground truth: stim:logday=+0.143, ratSD=0.0361, resSD=0.09377) ---
|
||||
|
||||
true stim:log(day) = +0.143 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.41 | 0.80 <- observed
|
||||
5 | 0.80 | 0.93
|
||||
8 | 0.99 | 1.00
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +0.07149 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.12 | 0.29 <- observed
|
||||
5 | 0.33 | 0.42
|
||||
8 | 0.50 | 0.58
|
||||
12 | 0.72 | 0.72
|
||||
16 | 0.86 | 0.87
|
||||
24 | 0.95 | 0.97
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d0_5
|
||||
POWER SIMULATION -- boxa_a2_d0_5 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+4.68/day, rat SD=0.00, residual SD=12.01, days=6
|
||||
ground truth: stim:day=+4.676/day, rat SD=0, residual SD=12.01, days=6
|
||||
|
||||
true stim:day interaction = +4.68 (100% of observed)
|
||||
true stim:day interaction = +4.676 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.23 | 0.54 <- observed
|
||||
@@ -15,7 +15,7 @@ ground truth: stim:day=+4.68/day, rat SD=0.00, residual SD=12.01, days=6
|
||||
16 | 0.99 | 0.99
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +2.34 (50% of observed)
|
||||
true stim:day interaction = +2.338 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.08 | 0.20 <- observed
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+4.68/day, rat SD=0.00, residual SD=12.01, days=6
|
||||
16 | 0.56 | 0.60
|
||||
24 | 0.82 | 0.84
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d0_5 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+0.04308/day, rat SD=0.03963, residual SD=0.08479, days=6
|
||||
|
||||
true stim:day interaction = +0.04308 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.36 | 0.73 <- observed
|
||||
5 | 0.78 | 0.91
|
||||
8 | 0.98 | 0.98
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +0.02154 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.11 | 0.25 <- observed
|
||||
5 | 0.26 | 0.38
|
||||
8 | 0.44 | 0.46
|
||||
12 | 0.68 | 0.72
|
||||
16 | 0.81 | 0.81
|
||||
24 | 0.95 | 0.94
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d0_5
|
||||
VARIATION: boxa_a2_d0_5 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(38)= 2.13 F(1)= 4.551 p=0.03941 p=0.03878 (df=42)
|
||||
day (learning) t(38)= 7.35 F(1)= 53.982 p=8.377e-09 p=4.651e-09 (df=42)
|
||||
stim (main, window start) t(38)= -0.08 F(1)= 0.007 p=0.9342 p=0.9342 (df=42)
|
||||
interaction 95% CI: [+0.24, +9.11]
|
||||
interaction 95% CI: [+0.2389, +9.113]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=2.87, p=0.1341
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.03941, slope diff=+4.68)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.03941, slope diff=+4.676)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d0_5 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A, Electrode-Box-A2
|
||||
N = 7 rats, 42 sessions raw day coverage: treat 0..5, control 0..5
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 42
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 7
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-70.223 -59.797 41.112 -82.223
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF pValue
|
||||
{'(Intercept)'} 0.19895 0.036526 5.4467 38 3.2657e-06
|
||||
{'day' } 0.049304 0.010134 4.8651 38 2.0216e-05
|
||||
{'stim' } -0.034105 0.055794 -0.61127 38 0.54466
|
||||
{'day:stim' } 0.043083 0.01548 2.7831 38 0.008341
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.125 0.27289
|
||||
0.028788 0.06982
|
||||
-0.14705 0.078844
|
||||
0.011745 0.074421
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (7 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.039633
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.015473 0.10151
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.084789 0.067081 0.10717
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(38)= 2.78 F(1)= 7.746 p=0.008341 p=0.008621 (df=35)
|
||||
day (learning) t(38)= 4.87 F(1)= 23.669 p=2.022e-05 p=2.41e-05 (df=35)
|
||||
stim (main, window start) t(38)= -0.61 F(1)= 0.374 p=0.5447 p=0.5472 (df=22)
|
||||
interaction 95% CI: [+0.01175, +0.07442]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,9.1)=5.64, p=0.0412
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.008341, slope diff=+0.04308)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 35 KiB After Width: | Height: | Size: 36 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 38 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d6_10
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d6_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,26)=0.252 p(resid)=0.6199 p(Satt)=0.6202 (df=
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+9.522, ratSD=6.321, resSD=8.298) ---
|
||||
|
||||
true stim:log(day) = +9.52 (100% of observed)
|
||||
true stim:log(day) = +9.522 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.02 | 0.06 <- observed
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.004) = 0.065 (small; f: .10 small, .25
|
||||
16 | 0.14 | 0.16
|
||||
24 | 0.32 | 0.34
|
||||
|
||||
true stim:log(day) = +4.76 (50% of observed)
|
||||
true stim:log(day) = +4.761 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.06 | 0.10 <- observed
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.004) = 0.065 (small; f: .10 small, .25
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_a2_d6_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,26)=0.064 p(resid)=0.8023 p(Satt)=0.8025 (df=24)
|
||||
honest per-animal random slope (log-day): F(1,6.0)=0.04 p=0.8568
|
||||
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=-0.03226, ratSD=0.04467, resSD=0.05581) ---
|
||||
|
||||
true stim:log(day) = -0.03226 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.05 | 0.05 <- observed
|
||||
5 | 0.03 | 0.03
|
||||
8 | 0.07 | 0.08
|
||||
12 | 0.12 | 0.12
|
||||
16 | 0.15 | 0.12
|
||||
24 | 0.08 | 0.09
|
||||
|
||||
true stim:log(day) = -0.01613 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.03 | 0.07 <- observed
|
||||
5 | 0.07 | 0.04
|
||||
8 | 0.04 | 0.03
|
||||
12 | 0.06 | 0.07
|
||||
16 | 0.06 | 0.05
|
||||
24 | 0.07 | 0.07
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d6_10
|
||||
POWER SIMULATION -- boxa_a2_d6_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+1.00/day, rat SD=6.34, residual SD=8.21, days=5
|
||||
ground truth: stim:day=+1/day, rat SD=6.344, residual SD=8.208, days=5
|
||||
|
||||
true stim:day interaction = +1.00 (100% of observed)
|
||||
true stim:day interaction = +1 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.03 | 0.05 <- observed
|
||||
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.00/day, rat SD=6.34, residual SD=8.21, days=5
|
||||
16 | 0.12 | 0.14
|
||||
24 | 0.27 | 0.32
|
||||
|
||||
true stim:day interaction = +0.50 (50% of observed)
|
||||
true stim:day interaction = +0.5 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.05 | 0.09 <- observed
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.00/day, rat SD=6.34, residual SD=8.21, days=5
|
||||
16 | 0.08 | 0.11
|
||||
24 | 0.12 | 0.12
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_a2_d6_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=-0.003749/day, rat SD=0.04471, residual SD=0.05566, days=5
|
||||
|
||||
true stim:day interaction = -0.003749 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.03 | 0.04 <- observed
|
||||
5 | 0.03 | 0.04
|
||||
8 | 0.07 | 0.08
|
||||
12 | 0.11 | 0.13
|
||||
16 | 0.15 | 0.12
|
||||
24 | 0.09 | 0.09
|
||||
|
||||
true stim:day interaction = -0.001874 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.03 | 0.07 <- observed
|
||||
5 | 0.07 | 0.04
|
||||
8 | 0.04 | 0.03
|
||||
12 | 0.05 | 0.07
|
||||
16 | 0.06 | 0.05
|
||||
24 | 0.07 | 0.06
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d6_10
|
||||
VARIATION: boxa_a2_d6_10 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(26)= 0.47 F(1)= 0.223 p=0.6409 p=0.6413 (df=24)
|
||||
day (learning) t(26)= 2.11 F(1)= 4.466 p=0.04433 p=0.04517 (df=24)
|
||||
stim (main, window start) t(26)= 1.73 F(1)= 2.983 p=0.09599 p=0.1083 (df=13)
|
||||
interaction 95% CI: [-3.36, +5.36]
|
||||
interaction 95% CI: [-3.356, +5.356]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.18, p=0.6872
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.6409, slope diff=+1.00)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.6409, slope diff=+1)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_a2_d6_10 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A, Electrode-Box-A2
|
||||
N = 6 rats, 30 sessions raw day coverage: treat 6..10, control 6..10
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 30
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 6
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-67.521 -59.114 39.761 -79.521
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF
|
||||
{'(Intercept)'} 0.54193 0.035861 15.112 26
|
||||
{'day' } 0.015949 0.010163 1.5693 26
|
||||
{'stim' } 0.10168 0.050715 2.0049 26
|
||||
{'day:stim' } -0.0037489 0.014373 -0.26083 26
|
||||
|
||||
|
||||
pValue Lower Upper
|
||||
2.1656e-14 0.46821 0.61564
|
||||
0.12868 -0.0049418 0.036839
|
||||
0.05549 -0.0025676 0.20592
|
||||
0.79628 -0.033292 0.025795
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (6 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.044709
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.021196 0.094306
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.055665 0.041949 0.073866
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(26)= -0.26 F(1)= 0.068 p=0.7963 p=0.7964 (df=24)
|
||||
day (learning) t(26)= 1.57 F(1)= 2.463 p=0.1287 p=0.1297 (df=24)
|
||||
stim (main, window start) t(26)= 2.00 F(1)= 4.020 p=0.05549 p=0.06743 (df=12)
|
||||
interaction 95% CI: [-0.03329, +0.02579]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.04, p=0.8467
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.7963, slope diff=-0.003749)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 42 KiB After Width: | Height: | Size: 42 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 43 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_10
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,68)=5.859 p(resid)=0.01817 p(Satt)=0.0182 (df
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+10.4, ratSD=5.517, resSD=12.53) ---
|
||||
|
||||
true stim:log(day) = +10.40 (100% of observed)
|
||||
true stim:log(day) = +10.4 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.32 | 0.62 <- observed
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.013) = 0.115 (small-medium; f: .10 smal
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +5.20 (50% of observed)
|
||||
true stim:log(day) = +5.199 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.10 | 0.20 <- observed
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.013) = 0.115 (small-medium; f: .10 smal
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,68)=8.517 p(resid)=0.004766 p(Satt)=0.004783 (df=67)
|
||||
honest per-animal random slope (log-day): F(1,40.6)=8.18 p=0.006648
|
||||
Cohen's f (interaction, partial eta^2=0.026) = 0.164 (small-medium; f: .10 small, .25 medium, .40 large)
|
||||
|
||||
--- power simulation (log-day ground truth: stim:logday=+0.08588, ratSD=0.04132, resSD=0.08572) ---
|
||||
|
||||
true stim:log(day) = +0.08588 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.43 | 0.79 <- observed
|
||||
5 | 0.89 | 0.95
|
||||
8 | 0.98 | 0.99
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +0.04294 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.13 | 0.28 <- observed
|
||||
5 | 0.30 | 0.43
|
||||
8 | 0.65 | 0.62
|
||||
12 | 0.79 | 0.82
|
||||
16 | 0.93 | 0.93
|
||||
24 | 0.99 | 0.99
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d0_10
|
||||
POWER SIMULATION -- boxa_b2_d0_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+1.60/day, rat SD=6.05, residual SD=12.42, days=11
|
||||
ground truth: stim:day=+1.604/day, rat SD=6.047, residual SD=12.42, days=11
|
||||
|
||||
true stim:day interaction = +1.60 (100% of observed)
|
||||
true stim:day interaction = +1.604 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.17 | 0.33 <- observed
|
||||
@@ -15,7 +15,7 @@ ground truth: stim:day=+1.60/day, rat SD=6.05, residual SD=12.42, days=11
|
||||
16 | 0.95 | 0.97
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +0.80 (50% of observed)
|
||||
true stim:day interaction = +0.802 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.11 | 0.11 <- observed
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.60/day, rat SD=6.05, residual SD=12.42, days=11
|
||||
16 | 0.51 | 0.53
|
||||
24 | 0.67 | 0.66
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d0_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+0.01317/day, rat SD=0.04399, residual SD=0.08748, days=11
|
||||
|
||||
true stim:day interaction = +0.01317 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.23 | 0.45 <- observed
|
||||
5 | 0.56 | 0.71
|
||||
8 | 0.82 | 0.86
|
||||
12 | 0.99 | 0.99
|
||||
16 | 0.98 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +0.006586 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.12 | 0.13 <- observed
|
||||
5 | 0.12 | 0.19
|
||||
8 | 0.33 | 0.34
|
||||
12 | 0.38 | 0.42
|
||||
16 | 0.64 | 0.64
|
||||
24 | 0.77 | 0.77
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d0_10
|
||||
VARIATION: boxa_b2_d0_10 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(68)= 1.63 F(1)= 2.647 p=0.1084 p=0.1084 (df=68)
|
||||
day (learning) t(68)= 10.15 F(1)=103.100 p=2.916e-15 p=2.368e-15 (df=69)
|
||||
stim (main, window start) t(68)= 0.26 F(1)= 0.068 p=0.7946 p=0.7968 (df=18)
|
||||
interaction 95% CI: [-0.36, +3.57]
|
||||
interaction 95% CI: [-0.3634, +3.571]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,22.2)=2.05, p=0.1666
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.1084, slope diff=+1.60)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.1084, slope diff=+1.604)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d0_10 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A2
|
||||
N = 7 rats, 72 sessions raw day coverage: treat 0..10, control 0..10
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 72
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 7
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-125.6 -111.94 68.799 -137.6
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF pValue
|
||||
{'(Intercept)'} 0.22991 0.038687 5.9428 68 1.07e-07
|
||||
{'day' } 0.040788 0.0055544 7.3434 68 3.402e-10
|
||||
{'stim' } 0.0034269 0.050885 0.067345 68 0.9465
|
||||
{'day:stim' } 0.013172 0.0069458 1.8964 68 0.062154
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.15271 0.30711
|
||||
0.029705 0.051872
|
||||
-0.098113 0.10497
|
||||
-0.00068802 0.027032
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (7 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.043995
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.020897 0.092622
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.087482 0.073647 0.10392
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(68)= 1.90 F(1)= 3.596 p=0.06215 p=0.06217 (df=68)
|
||||
day (learning) t(68)= 7.34 F(1)= 53.926 p=3.402e-10 p=3.124e-10 (df=69)
|
||||
stim (main, window start) t(68)= 0.07 F(1)= 0.005 p=0.9465 p=0.9471 (df=18)
|
||||
interaction 95% CI: [-0.000688, +0.02703]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,62.1)=3.56, p=0.06388
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.06215, slope diff=+0.01317)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 36 KiB After Width: | Height: | Size: 36 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 38 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_5
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_5 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,38)=3.165 p(resid)=0.08322 p(Satt)=0.08246 (d
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+13.68, ratSD=0, resSD=14.91) ---
|
||||
|
||||
true stim:log(day) = +13.68 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.026) = 0.163 (small-medium; f: .10 smal
|
||||
16 | 0.96 | 0.97
|
||||
24 | 0.98 | 0.99
|
||||
|
||||
true stim:log(day) = +6.84 (50% of observed)
|
||||
true stim:log(day) = +6.838 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.14 <- observed
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.026) = 0.163 (small-medium; f: .10 smal
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d0_5 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,38)=6.782 p(resid)=0.01307 p(Satt)=0.01342 (df=35)
|
||||
honest per-animal random slope (log-day): F(1,9.8)=4.90 p=0.05187
|
||||
Cohen's f (interaction, partial eta^2=0.073) = 0.280 (medium; f: .10 small, .25 medium, .40 large)
|
||||
|
||||
--- power simulation (log-day ground truth: stim:logday=+0.1288, ratSD=0.045, resSD=0.09595) ---
|
||||
|
||||
true stim:log(day) = +0.1288 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.35 | 0.69 <- observed
|
||||
5 | 0.75 | 0.86
|
||||
8 | 0.96 | 0.99
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:log(day) = +0.0644 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.12 | 0.23 <- observed
|
||||
5 | 0.24 | 0.30
|
||||
8 | 0.38 | 0.42
|
||||
12 | 0.61 | 0.64
|
||||
16 | 0.77 | 0.78
|
||||
24 | 0.93 | 0.95
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d0_5
|
||||
POWER SIMULATION -- boxa_b2_d0_5 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+4.71/day, rat SD=5.87, residual SD=11.30, days=6
|
||||
ground truth: stim:day=+4.707/day, rat SD=5.872, residual SD=11.3, days=6
|
||||
|
||||
true stim:day interaction = +4.71 (100% of observed)
|
||||
true stim:day interaction = +4.707 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.27 | 0.59 <- observed
|
||||
@@ -15,7 +15,7 @@ ground truth: stim:day=+4.71/day, rat SD=5.87, residual SD=11.30, days=6
|
||||
16 | 0.99 | 0.99
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +2.35 (50% of observed)
|
||||
true stim:day interaction = +2.354 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.10 | 0.20 <- observed
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+4.71/day, rat SD=5.87, residual SD=11.30, days=6
|
||||
16 | 0.62 | 0.62
|
||||
24 | 0.85 | 0.87
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d0_5 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+0.0427/day, rat SD=0.04855, residual SD=0.08495, days=6
|
||||
|
||||
true stim:day interaction = +0.0427 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.35 | 0.74 <- observed
|
||||
5 | 0.78 | 0.91
|
||||
8 | 0.98 | 0.98
|
||||
12 | 1.00 | 1.00
|
||||
16 | 1.00 | 1.00
|
||||
24 | 1.00 | 1.00
|
||||
|
||||
true stim:day interaction = +0.02135 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.11 | 0.25 <- observed
|
||||
5 | 0.25 | 0.36
|
||||
8 | 0.44 | 0.47
|
||||
12 | 0.68 | 0.73
|
||||
16 | 0.80 | 0.82
|
||||
24 | 0.95 | 0.94
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d0_5
|
||||
VARIATION: boxa_b2_d0_5 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(38)= 2.28 F(1)= 5.206 p=0.02821 p=0.02871 (df=35)
|
||||
day (learning) t(38)= 6.32 F(1)= 39.947 p=2.071e-07 p=2.928e-07 (df=35)
|
||||
stim (main, window start) t(38)= -0.61 F(1)= 0.378 p=0.5425 p=0.5456 (df=20)
|
||||
interaction 95% CI: [+0.53, +8.88]
|
||||
interaction 95% CI: [+0.5306, +8.884]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=2.92, p=0.131
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.02821, slope diff=+4.71)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.02821, slope diff=+4.707)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d0_5 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A2
|
||||
N = 7 rats, 42 sessions raw day coverage: treat 0..5, control 0..5
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 42
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 7
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-68.33 -57.904 40.165 -80.33
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF pValue
|
||||
{'(Intercept)'} 0.2186 0.045229 4.8332 38 2.2327e-05
|
||||
{'day' } 0.043366 0.011725 3.6987 38 0.00068267
|
||||
{'stim' } -0.059971 0.059832 -1.0023 38 0.32253
|
||||
{'day:stim' } 0.042705 0.01551 2.7533 38 0.0089978
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.12704 0.31016
|
||||
0.01963 0.067101
|
||||
-0.1811 0.061153
|
||||
0.011306 0.074103
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (7 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.048546
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.021809 0.10806
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.084953 0.067211 0.10738
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(38)= 2.75 F(1)= 7.581 p=0.008998 p=0.00929 (df=35)
|
||||
day (learning) t(38)= 3.70 F(1)= 13.680 p=0.0006827 p=0.0007391 (df=35)
|
||||
stim (main, window start) t(38)= -1.00 F(1)= 1.005 p=0.3225 p=0.3289 (df=19)
|
||||
interaction 95% CI: [+0.01131, +0.0741]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,8.8)=5.05, p=0.05179
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates) (p=0.008998, slope diff=+0.0427)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
@@ -1,30 +1,58 @@
|
||||
% Variation analysis -- the paper's linear mixed model on the successful-reach
|
||||
% COUNT, fit on this folder's curated data subset.
|
||||
% Variation analysis -- the paper's linear mixed model on this folder's data,
|
||||
% for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> result.txt
|
||||
% metric = rate : behavior = success / attempts -> result_rate.txt
|
||||
% (rate uses only sessions with attempts > 0)
|
||||
%
|
||||
% model: behavior ~ stim + day + stim:day + (1|rat)
|
||||
% behavior = successful reaches (count per session)
|
||||
% stim = 1 for the treatment group(s), 0 for the control group(s)
|
||||
% day = training day within this window (0 = first analyzed day)
|
||||
% rat = subject (random intercept)
|
||||
%
|
||||
% Self-contained: reads data.csv beside this script and writes result.txt.
|
||||
% Run headless from this folder with: matlab -batch "analyze"
|
||||
% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
% stim = 1 treatment / 0 control; day = training day within window (0 =
|
||||
% first analyzed day); rat = subject (random intercept).
|
||||
% For the interaction we report residual DF, Satterthwaite DF, and the honest
|
||||
% per-animal random-slope test. Self-contained: reads data.csv beside this
|
||||
% script. Run headless with: matlab -batch "analyze"
|
||||
% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
Rc = localAnalyze(D, 'count', here, vname);
|
||||
Rr = localAnalyze(D, 'rate', here, vname);
|
||||
|
||||
% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended).
|
||||
VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ...
|
||||
'interP', Rc.interP, 'interEst', Rc.interEst, ...
|
||||
'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ...
|
||||
'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ...
|
||||
'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs);
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function R = localAnalyze(D, metric, here, vname)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :);
|
||||
beh = D.success ./ D.total;
|
||||
mlabel = 'success RATE (success/attempts)'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success;
|
||||
mlabel = 'success COUNT'; suffix = '';
|
||||
end
|
||||
R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ...
|
||||
'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'covEqual', false);
|
||||
if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2
|
||||
localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ...
|
||||
vname, mlabel), here, suffix);
|
||||
return
|
||||
end
|
||||
|
||||
tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)');
|
||||
C = m.Coefficients; A = anova(m); ci = coefCI(m);
|
||||
As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF
|
||||
As = anova(m, 'DFMethod', 'satterthwaite');
|
||||
|
||||
% Honest test: refit with a per-animal random SLOPE so the interaction DF
|
||||
% collapses toward the animal count (guarded -- may not converge in short windows).
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
wst = warning('off', 'all');
|
||||
try
|
||||
@@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f
|
||||
C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ...
|
||||
As.pValue(gs(t)), As.DF2(gs(t)));
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1));
|
||||
maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0));
|
||||
if abs(maxT - maxC) > 2
|
||||
@@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2
|
||||
else
|
||||
cov = '(equal day coverage over this window)';
|
||||
end
|
||||
|
||||
ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii);
|
||||
if pI >= 0.05
|
||||
verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0
|
||||
verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else
|
||||
verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)';
|
||||
end
|
||||
if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)';
|
||||
elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)';
|
||||
else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
raw = regexprep(evalc('disp(m)'), '</?strong>', '');
|
||||
s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')];
|
||||
s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('day = training day within window (0 = first analyzed day)\n')];
|
||||
s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))];
|
||||
s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))];
|
||||
@@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res
|
||||
s = [s row('stim x day (interaction)', 'day:stim')];
|
||||
s = [s row('day (learning)', 'day')];
|
||||
s = [s row('stim (main, window start)', 'stim')];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))];
|
||||
s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))];
|
||||
if rsOk
|
||||
s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)];
|
||||
else
|
||||
@@ -82,17 +105,18 @@ else
|
||||
end
|
||||
s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ...
|
||||
' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)];
|
||||
s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')];
|
||||
|
||||
localWrite(s, here, suffix);
|
||||
|
||||
R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ...
|
||||
'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2);
|
||||
end
|
||||
|
||||
function localWrite(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'result.txt'), 'w');
|
||||
fprintf(fid, '%s', s);
|
||||
fclose(fid);
|
||||
|
||||
% Machine-readable handoff for the summary table (see make_variations.m).
|
||||
VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ...
|
||||
'nObs', height(D), 'interP', pI, 'interEst', eI, ...
|
||||
'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ...
|
||||
'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ...
|
||||
'covEqual', abs(maxT - maxC) <= 2);
|
||||
fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 35 KiB After Width: | Height: | Size: 35 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 37 KiB |
@@ -1,7 +1,7 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d6_10
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d6_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count))
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
|
||||
|
||||
@@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,26)=0.850 p(resid)=0.365 p(Satt)=0.3657 (df=2
|
||||
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) ---
|
||||
--- power simulation (log-day ground truth: stim:logday=+18.32, ratSD=6.212, resSD=8.198) ---
|
||||
|
||||
true stim:log(day) = +18.32 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
@@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.014) = 0.118 (small-medium; f: .10 smal
|
||||
16 | 0.57 | 0.63
|
||||
24 | 0.75 | 0.78
|
||||
|
||||
true stim:log(day) = +9.16 (50% of observed)
|
||||
true stim:log(day) = +9.162 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.07 | 0.11
|
||||
@@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.014) = 0.118 (small-medium; f: .10 smal
|
||||
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.
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
|
||||
@@ -0,0 +1,35 @@
|
||||
==============================================================================
|
||||
LOG-DAY MODEL + COHEN'S f + POWER -- boxa_b2_d6_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE)
|
||||
log(day) uses 1-indexed training day (our day 0 = paper "Day 1")
|
||||
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
|
||||
|
||||
--- fitted on real data ---
|
||||
stim x log(day) interaction: F(1,26)=0.195 p(resid)=0.6626 p(Satt)=0.6629 (df=24)
|
||||
honest per-animal random slope (log-day): F(1,6.0)=0.11 p=0.7527
|
||||
Cohen's f (interaction, partial eta^2=0.004) = 0.059 (small; f: .10 small, .25 medium, .40 large)
|
||||
|
||||
--- power simulation (log-day ground truth: stim:logday=+0.05955, ratSD=0.03763, resSD=0.05566) ---
|
||||
|
||||
true stim:log(day) = +0.05955 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.02 | 0.04
|
||||
5 | 0.06 | 0.12
|
||||
8 | 0.11 | 0.09
|
||||
12 | 0.11 | 0.14
|
||||
16 | 0.12 | 0.15
|
||||
24 | 0.27 | 0.31
|
||||
|
||||
true stim:log(day) = +0.02978 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.06 | 0.09
|
||||
5 | 0.07 | 0.07
|
||||
8 | 0.04 | 0.07
|
||||
12 | 0.09 | 0.08
|
||||
16 | 0.09 | 0.10
|
||||
24 | 0.12 | 0.15
|
||||
|
||||
Read per-animal as the honest power; LME matches the paper's power code (optimistic).
|
||||
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,11 +1,12 @@
|
||||
% Variation learning-curve plot, in the style of the paper:
|
||||
% Variation learning-curve plots, in the style of the paper:
|
||||
% "Lines indicate mean (and SEM) across animals in the anodal (red) and
|
||||
% control (blue) groups."
|
||||
% Plots mean +/- SEM successful reaches per training day for the treatment /
|
||||
% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading
|
||||
% this folder's data.csv and saving learning_curve.png. The per-group N is read
|
||||
% from the data (each variation pools different groups), so the legend shows the
|
||||
% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1").
|
||||
% Produces TWO figures from this folder's data.csv:
|
||||
% learning_curve.png # successes (count) per training day
|
||||
% learning_curve_rate.png success rate (success/attempts) per training day
|
||||
% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is
|
||||
% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and
|
||||
% the x-axis tick labels are drawn vertically.
|
||||
% Run: matlab -batch "plotcurve"
|
||||
% (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.)
|
||||
|
||||
@@ -14,15 +15,20 @@ if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
days = unique(D.day); % 0-indexed
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
|
||||
days = unique(D.day);
|
||||
xd = days + 1; % plot as 1-indexed training day (paper axis)
|
||||
red = [0.85 0.10 0.10];
|
||||
blue = [0.10 0.30 0.85];
|
||||
|
||||
[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1)
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0)
|
||||
localPlot(D, days, xd, 'count', '# successes', ...
|
||||
fullfile(here, 'learning_curve.png'), vname, red, blue);
|
||||
localPlot(D, days, xd, 'rate', 'success rate', ...
|
||||
fullfile(here, 'learning_curve_rate.png'), vname, red, blue);
|
||||
|
||||
% ------------------------------------------------------------------ helpers
|
||||
function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue)
|
||||
[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment
|
||||
[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control
|
||||
fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]);
|
||||
hold on
|
||||
e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2);
|
||||
@@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid
|
||||
hold off
|
||||
legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ...
|
||||
'Location', 'northwest', 'Box', 'off');
|
||||
xlabel('training day');
|
||||
ylabel('# successes');
|
||||
xlabel('training day'); ylabel(ylab);
|
||||
title(vname, 'Interpreter', 'none');
|
||||
set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off');
|
||||
|
||||
outFile = fullfile(here, 'learning_curve.png');
|
||||
xtickangle(90); % vertical x-axis tick labels
|
||||
exportgraphics(fig, outFile, 'Resolution', 150);
|
||||
close(fig);
|
||||
fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc);
|
||||
fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc);
|
||||
end
|
||||
|
||||
% ------------------------------------------------------------------ helper
|
||||
function [M, S, n] = localCurve(D, stimVal, days)
|
||||
%LOCALCURVE Per-day mean and SEM of successes across the animals in a group.
|
||||
function [M, S, n] = localCurve(D, stimVal, days, metric)
|
||||
%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric.
|
||||
subs = unique(D.subject(D.stim == stimVal));
|
||||
n = numel(subs);
|
||||
X = nan(numel(days), n);
|
||||
for j = 1:n
|
||||
for i = 1:numel(days)
|
||||
r = D.subject == subs(j) & D.day == days(i);
|
||||
if any(r); X(i, j) = mean(D.success(r)); end
|
||||
if ~any(r); continue; end
|
||||
if strcmp(metric, 'rate')
|
||||
tot = sum(D.total(r));
|
||||
if tot > 0; X(i, j) = sum(D.success(r)) / tot; end
|
||||
else
|
||||
X(i, j) = mean(D.success(r));
|
||||
end
|
||||
end
|
||||
end
|
||||
M = mean(X, 2, 'omitnan');
|
||||
|
||||
@@ -1,11 +1,11 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d6_10
|
||||
POWER SIMULATION -- boxa_b2_d6_10 [metric: # successes (count)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window)
|
||||
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+1.90/day, rat SD=6.23, residual SD=8.12, days=5
|
||||
ground truth: stim:day=+1.9/day, rat SD=6.231, residual SD=8.123, days=5
|
||||
|
||||
true stim:day interaction = +1.90 (100% of observed)
|
||||
true stim:day interaction = +1.9 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.04 | 0.17
|
||||
@@ -25,5 +25,4 @@ ground truth: stim:day=+1.90/day, rat SD=6.23, residual SD=8.12, days=5
|
||||
16 | 0.18 | 0.19
|
||||
24 | 0.25 | 0.28
|
||||
|
||||
Read the per-animal column as the honest power. At the observed N this study
|
||||
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
==============================================================================
|
||||
POWER SIMULATION -- boxa_b2_d6_10 [metric: success RATE]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE; day within-window)
|
||||
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
|
||||
ground truth: stim:day=+0.006202/day, rat SD=0.03766, residual SD=0.05555, days=5
|
||||
|
||||
true stim:day interaction = +0.006202 (100% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.03 | 0.04
|
||||
5 | 0.08 | 0.13
|
||||
8 | 0.08 | 0.08
|
||||
12 | 0.10 | 0.13
|
||||
16 | 0.12 | 0.11
|
||||
24 | 0.24 | 0.24
|
||||
|
||||
true stim:day interaction = +0.003101 (50% of observed)
|
||||
N/group | per-animal power | LME power
|
||||
------------------------------------------
|
||||
3 | 0.05 | 0.09
|
||||
5 | 0.07 | 0.07
|
||||
8 | 0.05 | 0.07
|
||||
12 | 0.07 | 0.08
|
||||
16 | 0.08 | 0.10
|
||||
24 | 0.10 | 0.12
|
||||
|
||||
Read the per-animal column as the honest power; LME is optimistic (obs-level DF).
|
||||
@@ -1,45 +1,49 @@
|
||||
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||
% interaction, using THIS folder's data as the ground truth.
|
||||
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> power_result.txt
|
||||
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
||||
%
|
||||
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
||||
% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD,
|
||||
% and residual SD generate NREP synthetic datasets at each rats-per-group N and
|
||||
% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is
|
||||
% scored at alpha = 0.05 two ways:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the
|
||||
% honest power, matching the random-slope / per-animal inference)
|
||||
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
||||
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
||||
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
||||
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
||||
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||
% Writes power_result.txt beside this script. Run: matlab -batch "powersim"
|
||||
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
||||
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
|
||||
localPower(D, 'count', here, vname);
|
||||
localPower(D, 'rate', here, vname);
|
||||
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
day0 = min(D.day);
|
||||
tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)];
|
||||
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
lme = fitlme(tbl0, FORMULA);
|
||||
@@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
||||
days = (0:max(tbl0.day))';
|
||||
|
||||
s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ...
|
||||
s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
||||
bInt, sRat, sRes, numel(days))];
|
||||
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -70,20 +73,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d6_10
|
||||
VARIATION: boxa_b2_d6_10 [metric: success COUNT]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)
|
||||
day = training day within window (0 = first analyzed day)
|
||||
@@ -60,9 +60,9 @@ effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
stim x day (interaction) t(26)= 0.85 F(1)= 0.729 p=0.4009 p=0.4015 (df=24)
|
||||
day (learning) t(26)= 1.32 F(1)= 1.746 p=0.1979 p=0.1989 (df=24)
|
||||
stim (main, window start) t(26)= 1.56 F(1)= 2.448 p=0.1297 p=0.142 (df=13)
|
||||
interaction 95% CI: [-2.67, +6.47]
|
||||
interaction 95% CI: [-2.673, +6.473]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.61, p=0.4631
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.4009, slope diff=+1.90)
|
||||
Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.4009, slope diff=+1.9)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
|
||||
@@ -0,0 +1,68 @@
|
||||
==============================================================================
|
||||
VARIATION: boxa_b2_d6_10 [metric: success RATE (success/attempts)]
|
||||
==============================================================================
|
||||
model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts))
|
||||
day = training day within window (0 = first analyzed day)
|
||||
treatment (stim=1): Electrode-Box-A, Electrode-Box-B2
|
||||
control (stim=0): Electrode-Box-A2
|
||||
N = 6 rats, 30 sessions raw day coverage: treat 6..10, control 6..10
|
||||
(equal day coverage over this window)
|
||||
|
||||
==============================================================================
|
||||
FULL MODEL SUMMARY -- fitlme
|
||||
==============================================================================
|
||||
|
||||
Linear mixed-effects model fit by ML
|
||||
|
||||
Model information:
|
||||
Number of observations 30
|
||||
Fixed effects coefficients 4
|
||||
Random effects coefficients 6
|
||||
Covariance parameters 2
|
||||
|
||||
Formula:
|
||||
behavior ~ 1 + day*stim + (1 | rat)
|
||||
|
||||
Model fit statistics:
|
||||
AIC BIC LogLikelihood Deviance
|
||||
-69.131 -60.724 40.566 -81.131
|
||||
|
||||
Fixed effects coefficients (95% CIs):
|
||||
Name Estimate SE tStat DF pValue
|
||||
{'(Intercept)'} 0.52619 0.040432 13.014 26 6.8005e-13
|
||||
{'day' } 0.0099391 0.012422 0.80014 26 0.43088
|
||||
{'stim' } 0.099864 0.049519 2.0167 26 0.054168
|
||||
{'day:stim' } 0.0062025 0.015213 0.4077 26 0.68683
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.44308 0.6093
|
||||
-0.015594 0.035472
|
||||
-0.0019243 0.20165
|
||||
-0.025069 0.037474
|
||||
|
||||
Random effects covariance parameters (95% CIs):
|
||||
Group: rat (6 Levels)
|
||||
Name1 Name2 Type Estimate
|
||||
{'(Intercept)'} {'(Intercept)'} {'std'} 0.037656
|
||||
|
||||
|
||||
Lower Upper
|
||||
0.016562 0.085613
|
||||
|
||||
Group: Error
|
||||
Name Estimate Lower Upper
|
||||
{'Res Std'} 0.055552 0.041864 0.073715
|
||||
|
||||
|
||||
effect t(df) / F(df1) p (resid) Satterthwaite: p (df)
|
||||
----------------------------------------------------------------------------
|
||||
stim x day (interaction) t(26)= 0.41 F(1)= 0.166 p=0.6868 p=0.6871 (df=24)
|
||||
day (learning) t(26)= 0.80 F(1)= 0.640 p=0.4309 p=0.4315 (df=24)
|
||||
stim (main, window start) t(26)= 2.02 F(1)= 4.067 p=0.05417 p=0.06308 (df=14)
|
||||
interaction 95% CI: [-0.02507, +0.03747]
|
||||
HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.10, p=0.7623
|
||||
(Satterthwaite DF ~= residual on this random-intercept model; the random-slope
|
||||
model above is the honest learning-rate test -- DF collapses toward the animal count.)
|
||||
INTERPRETATION: stim x day interaction n.s. -- slopes parallel (no differential learning rate) (p=0.6868, slope diff=+0.006202)
|
||||
Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.
|
||||
Reference in New Issue
Block a user