diff --git a/analysis/matlab/variations/right_only_d0_10/analyze.m b/analysis/matlab/variations/right_only_d0_10/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/right_only_d0_10/analyze.m +++ b/analysis/matlab/variations/right_only_d0_10/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/right_only_d0_10/learning_curve.png b/analysis/matlab/variations/right_only_d0_10/learning_curve.png index 7e11bfd..7354cbe 100644 Binary files a/analysis/matlab/variations/right_only_d0_10/learning_curve.png and b/analysis/matlab/variations/right_only_d0_10/learning_curve.png differ diff --git a/analysis/matlab/variations/right_only_d0_10/learning_curve_rate.png b/analysis/matlab/variations/right_only_d0_10/learning_curve_rate.png new file mode 100644 index 0000000..9c30c9e Binary files /dev/null and b/analysis/matlab/variations/right_only_d0_10/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/right_only_d0_10/logpower_result.txt b/analysis/matlab/variations/right_only_d0_10/logpower_result.txt index 2743f7f..8f55caf 100644 --- a/analysis/matlab/variations/right_only_d0_10/logpower_result.txt +++ b/analysis/matlab/variations/right_only_d0_10/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_10 +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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,7 +10,7 @@ stim x log(day) interaction: F(1,68)=8.420 p(resid)=0.004998 p(Satt)=0.005011 honest per-animal random slope (log-day): F(1,22.8)=7.55 p=0.01151 Cohen's f (interaction, partial eta^2=0.016) = 0.126 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+11.36, ratSD=4.91, resSD=11.42) --- +--- power simulation (log-day ground truth: stim:logday=+11.36, ratSD=4.911, resSD=11.42) --- true stim:log(day) = +11.36 (100% of observed) N/group | per-animal power | LME power @@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.016) = 0.126 (small-medium; f: .10 smal 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:log(day) = +5.68 (50% of observed) + true stim:log(day) = +5.682 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.13 | 0.29 <- observed @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.016) = 0.126 (small-medium; f: .10 smal 16 | 0.93 | 0.93 24 | 0.99 | 0.99 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_10/logpower_result_rate.txt b/analysis/matlab/variations/right_only_d0_10/logpower_result_rate.txt new file mode 100644 index 0000000..af17f14 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_10/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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)=13.480 p(resid)=0.0004756 p(Satt)=0.0004808 (df=67) + honest per-animal random slope (log-day): F(1,62.4)=13.42 p=0.0005158 +Cohen's f (interaction, partial eta^2=0.033) = 0.185 (small-medium; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.095, ratSD=0.04259, resSD=0.07524) --- + + true stim:log(day) = +0.095 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.63 | 0.95 <- observed + 5 | 0.99 | 1.00 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:log(day) = +0.0475 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.23 | 0.42 <- observed + 5 | 0.55 | 0.66 + 8 | 0.78 | 0.82 + 12 | 0.96 | 0.97 + 16 | 0.98 | 1.00 + 24 | 0.99 | 1.00 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_10/logpowersim.m b/analysis/matlab/variations/right_only_d0_10/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/right_only_d0_10/logpowersim.m +++ b/analysis/matlab/variations/right_only_d0_10/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_10/plotcurve.m b/analysis/matlab/variations/right_only_d0_10/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/right_only_d0_10/plotcurve.m +++ b/analysis/matlab/variations/right_only_d0_10/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/right_only_d0_10/power_result.txt b/analysis/matlab/variations/right_only_d0_10/power_result.txt index 01766dd..e6c5418 100644 --- a/analysis/matlab/variations/right_only_d0_10/power_result.txt +++ b/analysis/matlab/variations/right_only_d0_10/power_result.txt @@ -1,9 +1,9 @@ ============================================================================== -POWER SIMULATION -- right_only_d0_10 +POWER SIMULATION -- right_only_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.73/day, rat SD=5.35, residual SD=11.96, days=11 +ground truth: stim:day=+1.73/day, rat SD=5.346, residual SD=11.96, days=11 true stim:day interaction = +1.73 (100% of observed) N/group | per-animal power | LME power @@ -15,7 +15,7 @@ ground truth: stim:day=+1.73/day, rat SD=5.35, residual SD=11.96, days=11 16 | 0.98 | 0.98 24 | 1.00 | 1.00 - true stim:day interaction = +0.87 (50% of observed) + true stim:day interaction = +0.865 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.12 | 0.12 <- observed @@ -25,5 +25,4 @@ ground truth: stim:day=+1.73/day, rat SD=5.35, residual SD=11.96, days=11 16 | 0.62 | 0.62 24 | 0.74 | 0.71 -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). diff --git a/analysis/matlab/variations/right_only_d0_10/power_result_rate.txt b/analysis/matlab/variations/right_only_d0_10/power_result_rate.txt new file mode 100644 index 0000000..cf79116 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_10/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- right_only_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.01393/day, rat SD=0.04403, residual SD=0.08402, days=11 + + true stim:day interaction = +0.01393 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.28 | 0.53 <- observed + 5 | 0.65 | 0.81 + 8 | 0.89 | 0.90 + 12 | 0.99 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:day interaction = +0.006967 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.13 | 0.16 <- observed + 5 | 0.15 | 0.23 + 8 | 0.36 | 0.44 + 12 | 0.43 | 0.48 + 16 | 0.72 | 0.72 + 24 | 0.84 | 0.87 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/right_only_d0_10/powersim.m b/analysis/matlab/variations/right_only_d0_10/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/right_only_d0_10/powersim.m +++ b/analysis/matlab/variations/right_only_d0_10/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_10/result.txt b/analysis/matlab/variations/right_only_d0_10/result.txt index a1b04e2..49bda1e 100644 --- a/analysis/matlab/variations/right_only_d0_10/result.txt +++ b/analysis/matlab/variations/right_only_d0_10/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: right_only_d0_10 +VARIATION: right_only_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.83 F(1)= 3.336 p=0.07218 p=0.07218 (df=68) day (learning) t(68)= 10.55 F(1)=111.222 p=5.943e-16 p=4.533e-16 (df=70) stim (main, window start) t(68)= 0.37 F(1)= 0.134 p=0.7159 p=0.7187 (df=19) -interaction 95% CI: [-0.16, +3.62] +interaction 95% CI: [-0.1602, +3.62] HONEST LME (per-animal random slope, day|rat): interaction F(1,20.5)=2.56, p=0.1246 (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.07218, slope diff=+1.73) -Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008. +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d0_10/result_rate.txt b/analysis/matlab/variations/right_only_d0_10/result_rate.txt new file mode 100644 index 0000000..fe2622a --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_10/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: right_only_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, Right-Electrode +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 + -130.99 -117.33 71.495 -142.99 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF + {'(Intercept)'} 0.22977 0.037841 6.0719 68 + {'day' } 0.040867 0.0053407 7.652 68 + {'stim' } 0.00041563 0.049782 0.008349 68 + {'day:stim' } 0.013934 0.006676 2.0872 68 + + + pValue Lower Upper + 6.3632e-08 0.15426 0.30528 + 9.3879e-11 0.03021 0.051524 + 0.99336 -0.098922 0.099754 + 0.040617 0.00061265 0.027256 + +Random effects covariance parameters (95% CIs): +Group: rat (7 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.044032 + + + Lower Upper + 0.021261 0.09119 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.084023 0.070735 0.099809 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(68)= 2.09 F(1)= 4.357 p=0.04062 p=0.04064 (df=68) +day (learning) t(68)= 7.65 F(1)= 58.553 p=9.388e-11 p=8.647e-11 (df=69) +stim (main, window start) t(68)= 0.01 F(1)= 0.000 p=0.9934 p=0.9934 (df=17) +interaction 95% CI: [+0.0006126, +0.02726] +HONEST LME (per-animal random slope, day|rat): interaction F(1,67.1)=4.40, p=0.03968 + (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.04062, slope diff=+0.01393) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d0_13/analyze.m b/analysis/matlab/variations/right_only_d0_13/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/right_only_d0_13/analyze.m +++ b/analysis/matlab/variations/right_only_d0_13/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/right_only_d0_13/learning_curve.png b/analysis/matlab/variations/right_only_d0_13/learning_curve.png index fde70c7..f427c40 100644 Binary files a/analysis/matlab/variations/right_only_d0_13/learning_curve.png and b/analysis/matlab/variations/right_only_d0_13/learning_curve.png differ diff --git a/analysis/matlab/variations/right_only_d0_13/learning_curve_rate.png b/analysis/matlab/variations/right_only_d0_13/learning_curve_rate.png new file mode 100644 index 0000000..d54749e Binary files /dev/null and b/analysis/matlab/variations/right_only_d0_13/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/right_only_d0_13/logpower_result.txt b/analysis/matlab/variations/right_only_d0_13/logpower_result.txt index ff80186..c76cc4a 100644 --- a/analysis/matlab/variations/right_only_d0_13/logpower_result.txt +++ b/analysis/matlab/variations/right_only_d0_13/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_13 +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_13 [metric: # successes (count)] ============================================================================== -model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) +model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count)) log(day) uses 1-indexed training day (our day 0 = paper "Day 1") observed groups: stim n=4, control n=3 nrep=120, alpha=0.05 @@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,80)=13.182 p(resid)=0.0004969 p(Satt)=0.00049 honest per-animal random slope (log-day): F(1,18.9)=11.35 p=0.003237 Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+12.18, ratSD=5.40, resSD=10.82) --- +--- power simulation (log-day ground truth: stim:logday=+12.18, ratSD=5.4, resSD=10.82) --- true stim:log(day) = +12.18 (100% of observed) N/group | per-animal power | LME power @@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 smal 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:log(day) = +6.09 (50% of observed) + true stim:log(day) = +6.089 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.25 | 0.47 <- observed @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 smal 16 | 1.00 | 1.00 24 | 1.00 | 1.00 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_13/logpower_result_rate.txt b/analysis/matlab/variations/right_only_d0_13/logpower_result_rate.txt new file mode 100644 index 0000000..cdd2c9c --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_13/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_13 [metric: success RATE] +============================================================================== +model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE) +log(day) uses 1-indexed training day (our day 0 = paper "Day 1") +observed groups: stim n=4, control n=3 nrep=120, alpha=0.05 + +--- fitted on real data --- +stim x log(day) interaction: F(1,80)=20.061 p(resid)=2.465e-05 p(Satt)=2.466e-05 (df=80) + honest per-animal random slope (log-day): F(1,71.2)=19.98 p=2.882e-05 +Cohen's f (interaction, partial eta^2=0.035) = 0.191 (small-medium; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.1007, ratSD=0.04245, resSD=0.07232) --- + + true stim:log(day) = +0.1007 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.85 | 0.99 <- observed + 5 | 1.00 | 1.00 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:log(day) = +0.05035 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.40 | 0.66 <- observed + 5 | 0.74 | 0.86 + 8 | 0.97 | 0.99 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_13/logpowersim.m b/analysis/matlab/variations/right_only_d0_13/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/right_only_d0_13/logpowersim.m +++ b/analysis/matlab/variations/right_only_d0_13/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_13/plotcurve.m b/analysis/matlab/variations/right_only_d0_13/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/right_only_d0_13/plotcurve.m +++ b/analysis/matlab/variations/right_only_d0_13/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/right_only_d0_13/power_result.txt b/analysis/matlab/variations/right_only_d0_13/power_result.txt index ac5a438..29b02c8 100644 --- a/analysis/matlab/variations/right_only_d0_13/power_result.txt +++ b/analysis/matlab/variations/right_only_d0_13/power_result.txt @@ -1,11 +1,11 @@ ============================================================================== -POWER SIMULATION -- right_only_d0_13 +POWER SIMULATION -- right_only_d0_13 [metric: # successes (count)] ============================================================================== -model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) +model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window) observed groups: stim n=4, control n=3 nrep=120, alpha=0.05 -ground truth: stim:day=+1.70/day, rat SD=0.00, residual SD=14.12, days=14 +ground truth: stim:day=+1.7/day, rat SD=3.135e-15, residual SD=14.12, days=14 - true stim:day interaction = +1.70 (100% of observed) + true stim:day interaction = +1.7 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.28 | 0.67 <- observed @@ -15,7 +15,7 @@ ground truth: stim:day=+1.70/day, rat SD=0.00, residual SD=14.12, days=14 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:day interaction = +0.85 (50% of observed) + true stim:day interaction = +0.8501 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.15 | 0.24 <- observed @@ -25,5 +25,4 @@ ground truth: stim:day=+1.70/day, rat SD=0.00, residual SD=14.12, days=14 16 | 0.75 | 0.82 24 | 0.87 | 0.85 -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). diff --git a/analysis/matlab/variations/right_only_d0_13/power_result_rate.txt b/analysis/matlab/variations/right_only_d0_13/power_result_rate.txt new file mode 100644 index 0000000..ec546e3 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_13/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- right_only_d0_13 [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.01351/day, rat SD=0.03866, residual SD=0.09009, days=14 + + true stim:day interaction = +0.01351 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.37 | 0.82 <- observed + 5 | 0.85 | 0.93 + 8 | 0.99 | 0.99 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:day interaction = +0.006756 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.20 | 0.34 <- observed + 5 | 0.33 | 0.39 + 8 | 0.66 | 0.69 + 12 | 0.78 | 0.82 + 16 | 0.93 | 0.94 + 24 | 0.95 | 0.94 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/right_only_d0_13/powersim.m b/analysis/matlab/variations/right_only_d0_13/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/right_only_d0_13/powersim.m +++ b/analysis/matlab/variations/right_only_d0_13/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_13/result.txt b/analysis/matlab/variations/right_only_d0_13/result.txt index 752106e..a19fe47 100644 --- a/analysis/matlab/variations/right_only_d0_13/result.txt +++ b/analysis/matlab/variations/right_only_d0_13/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: right_only_d0_13 +VARIATION: right_only_d0_13 [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(80)= 2.01 F(1)= 4.031 p=0.04804 p=0.04788 (df=84) day (learning) t(80)= 9.45 F(1)= 89.224 p=1.171e-14 p=7.518e-15 (df=84) stim (main, window start) t(80)= 0.63 F(1)= 0.401 p=0.5285 p=0.5284 (df=84) -interaction 95% CI: [+0.02, +3.39] +interaction 95% CI: [+0.01501, +3.385] HONEST LME (per-animal random slope, day|rat): interaction F(1,20.2)=2.92, p=0.1029 (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.04804, slope diff=+1.70) -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.04804, slope diff=+1.7) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d0_13/result_rate.txt b/analysis/matlab/variations/right_only_d0_13/result_rate.txt new file mode 100644 index 0000000..e43cd19 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_13/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: right_only_d0_13 [metric: success RATE (success/attempts)] +============================================================================== +model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts)) +day = training day within window (0 = first analyzed day) +treatment (stim=1): Electrode-Box-B2, Right-Electrode +control (stim=0): Electrode-Box-A2 +N = 7 rats, 84 sessions raw day coverage: treat 0..13, control 0..13 +(equal day coverage over this window) + +============================================================================== +FULL MODEL SUMMARY -- fitlme +============================================================================== + +Linear mixed-effects model fit by ML + +Model information: + Number of observations 84 + 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 + -145.93 -131.35 78.967 -157.93 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF pValue + {'(Intercept)'} 0.25632 0.035875 7.1447 80 3.7358e-10 + {'day' } 0.032358 0.0046124 7.0154 80 6.6273e-10 + {'stim' } 0.0062182 0.046891 0.13261 80 0.89483 + {'day:stim' } 0.013512 0.0056067 2.41 80 0.018249 + + + Lower Upper + 0.18492 0.32771 + 0.023179 0.041537 + -0.087097 0.099534 + 0.0023542 0.02467 + +Random effects covariance parameters (95% CIs): +Group: rat (7 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.038665 + + + Lower Upper + 0.017826 0.083864 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.090086 0.076927 0.1055 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(80)= 2.41 F(1)= 5.808 p=0.01825 p=0.01818 (df=82) +day (learning) t(80)= 7.02 F(1)= 49.216 p=6.627e-10 p=5.568e-10 (df=83) +stim (main, window start) t(80)= 0.13 F(1)= 0.018 p=0.8948 p=0.8959 (df=19) +interaction 95% CI: [+0.002354, +0.02467] +HONEST LME (per-animal random slope, day|rat): interaction F(1,81.6)=5.93, p=0.01709 + (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.01825, slope diff=+0.01351) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d0_5/analyze.m b/analysis/matlab/variations/right_only_d0_5/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/right_only_d0_5/analyze.m +++ b/analysis/matlab/variations/right_only_d0_5/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/right_only_d0_5/learning_curve.png b/analysis/matlab/variations/right_only_d0_5/learning_curve.png index 01c8d4f..54d94cf 100644 Binary files a/analysis/matlab/variations/right_only_d0_5/learning_curve.png and b/analysis/matlab/variations/right_only_d0_5/learning_curve.png differ diff --git a/analysis/matlab/variations/right_only_d0_5/learning_curve_rate.png b/analysis/matlab/variations/right_only_d0_5/learning_curve_rate.png new file mode 100644 index 0000000..a1ba618 Binary files /dev/null and b/analysis/matlab/variations/right_only_d0_5/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/right_only_d0_5/logpower_result.txt b/analysis/matlab/variations/right_only_d0_5/logpower_result.txt index 963643d..eff0965 100644 --- a/analysis/matlab/variations/right_only_d0_5/logpower_result.txt +++ b/analysis/matlab/variations/right_only_d0_5/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d0_5 +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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,9 +10,9 @@ stim x log(day) interaction: F(1,38)=5.667 p(resid)=0.02241 p(Satt)=0.0219 (df honest per-animal random slope (log-day): F(1,9.0)=4.64 p=0.05966 Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+16.20, ratSD=0.00, resSD=13.20) --- +--- power simulation (log-day ground truth: stim:logday=+16.2, ratSD=0, resSD=13.2) --- - true stim:log(day) = +16.20 (100% of observed) + true stim:log(day) = +16.2 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.33 | 0.65 <- observed @@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 smal 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:log(day) = +8.10 (50% of observed) + true stim:log(day) = +8.101 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.09 | 0.24 <- observed @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 smal 16 | 0.69 | 0.69 24 | 0.92 | 0.92 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_5/logpower_result_rate.txt b/analysis/matlab/variations/right_only_d0_5/logpower_result_rate.txt new file mode 100644 index 0000000..d113d38 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_5/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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)=15.000 p(resid)=0.0004109 p(Satt)=0.0004502 (df=35) + honest per-animal random slope (log-day): F(1,10.8)=11.68 p=0.005855 +Cohen's f (interaction, partial eta^2=0.111) = 0.353 (medium; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.1583, ratSD=0.04781, resSD=0.07931) --- + + true stim:log(day) = +0.1583 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.62 | 0.93 <- observed + 5 | 0.95 | 0.98 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:log(day) = +0.07917 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.16 | 0.42 <- observed + 5 | 0.52 | 0.62 + 8 | 0.81 | 0.82 + 12 | 0.93 | 0.92 + 16 | 0.98 | 0.98 + 24 | 1.00 | 1.00 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d0_5/logpowersim.m b/analysis/matlab/variations/right_only_d0_5/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/right_only_d0_5/logpowersim.m +++ b/analysis/matlab/variations/right_only_d0_5/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_5/plotcurve.m b/analysis/matlab/variations/right_only_d0_5/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/right_only_d0_5/plotcurve.m +++ b/analysis/matlab/variations/right_only_d0_5/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/right_only_d0_5/power_result.txt b/analysis/matlab/variations/right_only_d0_5/power_result.txt index ce0bc1d..8baf17f 100644 --- a/analysis/matlab/variations/right_only_d0_5/power_result.txt +++ b/analysis/matlab/variations/right_only_d0_5/power_result.txt @@ -1,11 +1,11 @@ ============================================================================== -POWER SIMULATION -- right_only_d0_5 +POWER SIMULATION -- right_only_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=+5.31/day, rat SD=5.25, residual SD=9.96, days=6 +ground truth: stim:day=+5.307/day, rat SD=5.248, residual SD=9.964, days=6 - true stim:day interaction = +5.31 (100% of observed) + true stim:day interaction = +5.307 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.38 | 0.80 <- observed @@ -15,7 +15,7 @@ ground truth: stim:day=+5.31/day, rat SD=5.25, residual SD=9.96, days=6 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:day interaction = +2.65 (50% of observed) + true stim:day interaction = +2.654 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.11 | 0.28 <- observed @@ -25,5 +25,4 @@ ground truth: stim:day=+5.31/day, rat SD=5.25, residual SD=9.96, days=6 16 | 0.83 | 0.85 24 | 0.95 | 0.97 -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). diff --git a/analysis/matlab/variations/right_only_d0_5/power_result_rate.txt b/analysis/matlab/variations/right_only_d0_5/power_result_rate.txt new file mode 100644 index 0000000..f221424 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_5/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- right_only_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.04981/day, rat SD=0.04929, residual SD=0.07366, days=6 + + true stim:day interaction = +0.04981 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.60 | 0.91 <- observed + 5 | 0.97 | 0.98 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:day interaction = +0.0249 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.15 | 0.37 <- observed + 5 | 0.49 | 0.55 + 8 | 0.73 | 0.81 + 12 | 0.88 | 0.89 + 16 | 0.97 | 0.97 + 24 | 1.00 | 1.00 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/right_only_d0_5/powersim.m b/analysis/matlab/variations/right_only_d0_5/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/right_only_d0_5/powersim.m +++ b/analysis/matlab/variations/right_only_d0_5/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d0_5/result.txt b/analysis/matlab/variations/right_only_d0_5/result.txt index a2f22a5..630bf5f 100644 --- a/analysis/matlab/variations/right_only_d0_5/result.txt +++ b/analysis/matlab/variations/right_only_d0_5/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: right_only_d0_5 +VARIATION: right_only_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.92 F(1)= 8.511 p=0.005897 p=0.006128 (df=35) day (learning) t(38)= 7.17 F(1)= 51.382 p=1.458e-08 p=2.321e-08 (df=35) stim (main, window start) t(38)= -0.73 F(1)= 0.534 p=0.4696 p=0.4735 (df=20) -interaction 95% CI: [+1.62, +8.99] +interaction 95% CI: [+1.625, +8.99] HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=3.94, p=0.08763 (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.005897, slope diff=+5.31) -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.005897, slope diff=+5.307) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d0_5/result_rate.txt b/analysis/matlab/variations/right_only_d0_5/result_rate.txt new file mode 100644 index 0000000..ac77b29 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_5/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: right_only_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, Right-Electrode +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 + -78.771 -68.345 45.385 -90.771 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF pValue + {'(Intercept)'} 0.2186 0.041919 5.2147 38 6.7762e-06 + {'day' } 0.043366 0.010166 4.2656 38 0.00012743 + {'stim' } -0.074496 0.055454 -1.3434 38 0.18712 + {'day:stim' } 0.049806 0.013449 3.7033 38 0.00067357 + + + Lower Upper + 0.13374 0.30346 + 0.022785 0.063946 + -0.18676 0.037766 + 0.02258 0.077032 + +Random effects covariance parameters (95% CIs): +Group: rat (7 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.049289 + + + Lower Upper + 0.023894 0.10168 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.073663 0.058279 0.093108 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(38)= 3.70 F(1)= 13.715 p=0.0006736 p=0.0007295 (df=35) +day (learning) t(38)= 4.27 F(1)= 18.195 p=0.0001274 p=0.0001438 (df=35) +stim (main, window start) t(38)= -1.34 F(1)= 1.805 p=0.1871 p=0.1975 (df=16) +interaction 95% CI: [+0.02258, +0.07703] +HONEST LME (per-animal random slope, day|rat): interaction F(1,7.0)=9.42, p=0.0181 + (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.0006736, slope diff=+0.04981) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d6_10/analyze.m b/analysis/matlab/variations/right_only_d6_10/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/right_only_d6_10/analyze.m +++ b/analysis/matlab/variations/right_only_d6_10/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/right_only_d6_10/learning_curve.png b/analysis/matlab/variations/right_only_d6_10/learning_curve.png index beec2c1..8cedcdd 100644 Binary files a/analysis/matlab/variations/right_only_d6_10/learning_curve.png and b/analysis/matlab/variations/right_only_d6_10/learning_curve.png differ diff --git a/analysis/matlab/variations/right_only_d6_10/learning_curve_rate.png b/analysis/matlab/variations/right_only_d6_10/learning_curve_rate.png new file mode 100644 index 0000000..0b6182e Binary files /dev/null and b/analysis/matlab/variations/right_only_d6_10/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/right_only_d6_10/logpower_result.txt b/analysis/matlab/variations/right_only_d6_10/logpower_result.txt index 1a2bb3c..226dc36 100644 --- a/analysis/matlab/variations/right_only_d6_10/logpower_result.txt +++ b/analysis/matlab/variations/right_only_d6_10/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_10 +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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)=1.327 p(resid)=0.2599 p(Satt)=0.2608 (df= honest per-animal random slope (log-day): F(1,6.0)=0.96 p=0.3641 Cohen's f (interaction, partial eta^2=0.021) = 0.145 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+23.04, ratSD=5.69, resSD=8.25) --- +--- power simulation (log-day ground truth: stim:logday=+23.04, ratSD=5.689, resSD=8.25) --- true stim:log(day) = +23.04 (100% of observed) N/group | per-animal power | LME power @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.021) = 0.145 (small-medium; f: .10 smal 16 | 0.25 | 0.25 24 | 0.45 | 0.52 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d6_10/logpower_result_rate.txt b/analysis/matlab/variations/right_only_d6_10/logpower_result_rate.txt new file mode 100644 index 0000000..d1368af --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_10/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_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.485 p(resid)=0.4924 p(Satt)=0.4929 (df=24) + honest per-animal random slope (log-day): F(1,6.0)=0.24 p=0.644 +Cohen's f (interaction, partial eta^2=0.009) = 0.094 (small; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.09523, ratSD=0.03826, resSD=0.05641) --- + + true stim:log(day) = +0.09523 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.04 | 0.10 + 5 | 0.15 | 0.19 + 8 | 0.17 | 0.23 + 12 | 0.20 | 0.27 + 16 | 0.40 | 0.38 + 24 | 0.54 | 0.54 + + true stim:log(day) = +0.04761 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.06 | 0.11 + 5 | 0.09 | 0.09 + 8 | 0.06 | 0.07 + 12 | 0.12 | 0.12 + 16 | 0.14 | 0.14 + 24 | 0.22 | 0.22 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d6_10/logpowersim.m b/analysis/matlab/variations/right_only_d6_10/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/right_only_d6_10/logpowersim.m +++ b/analysis/matlab/variations/right_only_d6_10/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d6_10/plotcurve.m b/analysis/matlab/variations/right_only_d6_10/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/right_only_d6_10/plotcurve.m +++ b/analysis/matlab/variations/right_only_d6_10/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/right_only_d6_10/power_result.txt b/analysis/matlab/variations/right_only_d6_10/power_result.txt index 5374fde..fe348a6 100644 --- a/analysis/matlab/variations/right_only_d6_10/power_result.txt +++ b/analysis/matlab/variations/right_only_d6_10/power_result.txt @@ -1,11 +1,11 @@ ============================================================================== -POWER SIMULATION -- right_only_d6_10 +POWER SIMULATION -- right_only_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=+2.42/day, rat SD=5.71, residual SD=8.18, days=5 +ground truth: stim:day=+2.425/day, rat SD=5.709, residual SD=8.18, days=5 - true stim:day interaction = +2.42 (100% of observed) + true stim:day interaction = +2.425 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.06 | 0.25 @@ -15,7 +15,7 @@ ground truth: stim:day=+2.42/day, rat SD=5.71, residual SD=8.18, days=5 16 | 0.71 | 0.74 24 | 0.94 | 0.94 - true stim:day interaction = +1.21 (50% of observed) + true stim:day interaction = +1.212 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.07 | 0.12 @@ -25,5 +25,4 @@ ground truth: stim:day=+2.42/day, rat SD=5.71, residual SD=8.18, days=5 16 | 0.23 | 0.23 24 | 0.40 | 0.47 -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). diff --git a/analysis/matlab/variations/right_only_d6_10/power_result_rate.txt b/analysis/matlab/variations/right_only_d6_10/power_result_rate.txt new file mode 100644 index 0000000..bdd92d5 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_10/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- right_only_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.0102/day, rat SD=0.03828, residual SD=0.05632, days=5 + + true stim:day interaction = +0.0102 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.03 | 0.09 + 5 | 0.13 | 0.18 + 8 | 0.16 | 0.24 + 12 | 0.19 | 0.24 + 16 | 0.36 | 0.36 + 24 | 0.47 | 0.47 + + true stim:day interaction = +0.0051 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.07 | 0.11 + 5 | 0.07 | 0.10 + 8 | 0.06 | 0.08 + 12 | 0.10 | 0.12 + 16 | 0.12 | 0.13 + 24 | 0.19 | 0.19 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/right_only_d6_10/powersim.m b/analysis/matlab/variations/right_only_d6_10/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/right_only_d6_10/powersim.m +++ b/analysis/matlab/variations/right_only_d6_10/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d6_10/result.txt b/analysis/matlab/variations/right_only_d6_10/result.txt index e952e59..6c7527b 100644 --- a/analysis/matlab/variations/right_only_d6_10/result.txt +++ b/analysis/matlab/variations/right_only_d6_10/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: right_only_d6_10 +VARIATION: right_only_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) @@ -64,5 +64,5 @@ interaction 95% CI: [-2.18, +7.03] HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.90, p=0.3784 (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.289, slope diff=+2.42) -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.289, slope diff=+2.425) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d6_10/result_rate.txt b/analysis/matlab/variations/right_only_d6_10/result_rate.txt new file mode 100644 index 0000000..6f97191 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_10/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: right_only_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, Right-Electrode +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 + -68.28 -59.872 40.14 -80.28 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF pValue + {'(Intercept)'} 0.52619 0.041043 12.82 26 9.5461e-13 + {'day' } 0.0099391 0.012595 0.78916 26 0.43716 + {'stim' } 0.090314 0.050267 1.7967 26 0.084015 + {'day:stim' } 0.0102 0.015425 0.66127 26 0.51426 + + + Lower Upper + 0.44182 0.61055 + -0.015949 0.035828 + -0.013012 0.19364 + -0.021507 0.041907 + +Random effects covariance parameters (95% CIs): +Group: rat (6 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.038283 + + + Lower Upper + 0.016862 0.086916 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.056325 0.042446 0.074741 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(26)= 0.66 F(1)= 0.437 p=0.5143 p=0.5147 (df=24) +day (learning) t(26)= 0.79 F(1)= 0.623 p=0.4372 p=0.4377 (df=24) +stim (main, window start) t(26)= 1.80 F(1)= 3.228 p=0.08402 p=0.09376 (df=14) +interaction 95% CI: [-0.02151, +0.04191] +HONEST LME (per-animal random slope, day|rat): interaction F(1,6.0)=0.23, p=0.6495 + (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.5143, slope diff=+0.0102) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d6_13/analyze.m b/analysis/matlab/variations/right_only_d6_13/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/right_only_d6_13/analyze.m +++ b/analysis/matlab/variations/right_only_d6_13/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/right_only_d6_13/learning_curve.png b/analysis/matlab/variations/right_only_d6_13/learning_curve.png index d731fff..7b1fc7f 100644 Binary files a/analysis/matlab/variations/right_only_d6_13/learning_curve.png and b/analysis/matlab/variations/right_only_d6_13/learning_curve.png differ diff --git a/analysis/matlab/variations/right_only_d6_13/learning_curve_rate.png b/analysis/matlab/variations/right_only_d6_13/learning_curve_rate.png new file mode 100644 index 0000000..455bff5 Binary files /dev/null and b/analysis/matlab/variations/right_only_d6_13/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/right_only_d6_13/logpower_result.txt b/analysis/matlab/variations/right_only_d6_13/logpower_result.txt index 1f4088d..dea8a5c 100644 --- a/analysis/matlab/variations/right_only_d6_13/logpower_result.txt +++ b/analysis/matlab/variations/right_only_d6_13/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_13 +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_13 [metric: # successes (count)] ============================================================================== -model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT) +model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = # successes (count)) log(day) uses 1-indexed training day (our day 0 = paper "Day 1") observed groups: stim n=4, control n=2 nrep=120, alpha=0.05 @@ -10,7 +10,7 @@ stim x log(day) interaction: F(1,38)=5.774 p(resid)=0.02125 p(Satt)=0.02133 (d honest per-animal random slope (log-day): F(1,0.0)=3.19 p=NaN Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+28.42, ratSD=6.48, resSD=7.39) --- +--- power simulation (log-day ground truth: stim:logday=+28.42, ratSD=6.482, resSD=7.389) --- true stim:log(day) = +28.42 (100% of observed) N/group | per-animal power | LME power @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.036) = 0.194 (small-medium; f: .10 smal 16 | 0.90 | 0.92 24 | 1.00 | 1.00 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d6_13/logpower_result_rate.txt b/analysis/matlab/variations/right_only_d6_13/logpower_result_rate.txt new file mode 100644 index 0000000..82fa8b3 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_13/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- right_only_d6_13 [metric: success RATE] +============================================================================== +model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = success RATE) +log(day) uses 1-indexed training day (our day 0 = paper "Day 1") +observed groups: stim n=4, control n=2 nrep=120, alpha=0.05 + +--- fitted on real data --- +stim x log(day) interaction: F(1,38)=4.975 p(resid)=0.0317 p(Satt)=0.03178 (df=38) + honest per-animal random slope (log-day): F(1,4.9)=0.25 p=0.6396 +Cohen's f (interaction, partial eta^2=0.039) = 0.201 (small-medium; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.1853, ratSD=0.04282, resSD=0.05196) --- + + true stim:log(day) = +0.1853 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.53 | 0.80 + 5 | 0.88 | 0.96 + 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.09266 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.16 | 0.27 + 5 | 0.31 | 0.42 + 8 | 0.52 | 0.62 + 12 | 0.74 | 0.79 + 16 | 0.86 | 0.89 + 24 | 0.98 | 0.99 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/right_only_d6_13/logpowersim.m b/analysis/matlab/variations/right_only_d6_13/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/right_only_d6_13/logpowersim.m +++ b/analysis/matlab/variations/right_only_d6_13/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d6_13/plotcurve.m b/analysis/matlab/variations/right_only_d6_13/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/right_only_d6_13/plotcurve.m +++ b/analysis/matlab/variations/right_only_d6_13/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/right_only_d6_13/power_result.txt b/analysis/matlab/variations/right_only_d6_13/power_result.txt index 2e054dd..34247c9 100644 --- a/analysis/matlab/variations/right_only_d6_13/power_result.txt +++ b/analysis/matlab/variations/right_only_d6_13/power_result.txt @@ -1,11 +1,11 @@ ============================================================================== -POWER SIMULATION -- right_only_d6_13 +POWER SIMULATION -- right_only_d6_13 [metric: # successes (count)] ============================================================================== -model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window) +model: behavior ~ stim + day + stim:day + (1|rat) (behavior = # successes (count); day within-window) observed groups: stim n=4, control n=2 nrep=120, alpha=0.05 -ground truth: stim:day=+2.76/day, rat SD=6.47, residual SD=7.37, days=8 +ground truth: stim:day=+2.758/day, rat SD=6.469, residual SD=7.367, days=8 - true stim:day interaction = +2.76 (100% of observed) + true stim:day interaction = +2.758 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.56 | 0.83 @@ -15,7 +15,7 @@ ground truth: stim:day=+2.76/day, rat SD=6.47, residual SD=7.37, days=8 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:day interaction = +1.38 (50% of observed) + true stim:day interaction = +1.379 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.15 | 0.28 @@ -25,5 +25,4 @@ ground truth: stim:day=+2.76/day, rat SD=6.47, residual SD=7.37, days=8 16 | 0.93 | 0.93 24 | 0.99 | 0.99 -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). diff --git a/analysis/matlab/variations/right_only_d6_13/power_result_rate.txt b/analysis/matlab/variations/right_only_d6_13/power_result_rate.txt new file mode 100644 index 0000000..6a86d4c --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_13/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- right_only_d6_13 [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.01878/day, rat SD=0.04274, residual SD=0.05164, days=8 + + true stim:day interaction = +0.01878 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.54 | 0.82 + 5 | 0.92 | 0.97 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:day interaction = +0.009388 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.15 | 0.28 + 5 | 0.33 | 0.51 + 8 | 0.58 | 0.65 + 12 | 0.77 | 0.81 + 16 | 0.88 | 0.90 + 24 | 0.98 | 0.99 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/right_only_d6_13/powersim.m b/analysis/matlab/variations/right_only_d6_13/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/right_only_d6_13/powersim.m +++ b/analysis/matlab/variations/right_only_d6_13/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/right_only_d6_13/result.txt b/analysis/matlab/variations/right_only_d6_13/result.txt index e0a7616..d611f97 100644 --- a/analysis/matlab/variations/right_only_d6_13/result.txt +++ b/analysis/matlab/variations/right_only_d6_13/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: right_only_d6_13 +VARIATION: right_only_d6_13 [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.32 F(1)= 5.379 p=0.02585 p=0.02591 (df=38) day (learning) t(38)= 1.70 F(1)= 2.885 p=0.09761 p=0.09766 (df=38) stim (main, window start) t(38)= 1.64 F(1)= 2.685 p=0.1096 p=0.1304 (df=11) -interaction 95% CI: [+0.35, +5.17] +interaction 95% CI: [+0.3507, +5.165] HONEST LME (per-animal random slope, day|rat): interaction F(1,21.2)=4.96, p=0.03697 (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.02585, slope diff=+2.76) -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.02585, slope diff=+2.758) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/right_only_d6_13/result_rate.txt b/analysis/matlab/variations/right_only_d6_13/result_rate.txt new file mode 100644 index 0000000..5233cd6 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_13/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: right_only_d6_13 [metric: success RATE (success/attempts)] +============================================================================== +model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts)) +day = training day within window (0 = first analyzed day) +treatment (stim=1): Electrode-Box-B2, Right-Electrode +control (stim=0): Electrode-Box-A2 +N = 6 rats, 42 sessions raw day coverage: treat 6..13, control 6..13 +(equal day coverage over this window) + +============================================================================== +FULL MODEL SUMMARY -- fitlme +============================================================================== + +Linear mixed-effects model fit by ML + +Model information: + Number of observations 42 + 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 + -107.29 -96.867 59.646 -119.29 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF pValue + {'(Intercept)'} 0.53543 0.038925 13.755 38 2.4576e-16 + {'day' } 0.0034157 0.0070712 0.48304 38 0.63184 + {'stim' } 0.077975 0.047524 1.6407 38 0.1091 + {'day:stim' } 0.018776 0.0083266 2.2549 38 0.029981 + + + Lower Upper + 0.45663 0.61423 + -0.010899 0.017731 + -0.018233 0.17418 + 0.0019195 0.035632 + +Random effects covariance parameters (95% CIs): +Group: rat (6 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.042741 + + + Lower Upper + 0.021774 0.083895 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.051638 0.041022 0.065002 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(38)= 2.25 F(1)= 5.085 p=0.02998 p=0.03003 (df=38) +day (learning) t(38)= 0.48 F(1)= 0.233 p=0.6318 p=0.6318 (df=38) +stim (main, window start) t(38)= 1.64 F(1)= 2.692 p=0.1091 p=0.1287 (df=11) +interaction 95% CI: [+0.00192, +0.03563] +HONEST LME (per-animal random slope, day|rat): interaction F(1,4.0)=0.24, p=0.6492 + (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.02998, slope diff=+0.01878) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/unmerge_d0_10/analyze.m b/analysis/matlab/variations/unmerge_d0_10/analyze.m index 82da5ed..0d40402 100644 --- a/analysis/matlab/variations/unmerge_d0_10/analyze.m +++ b/analysis/matlab/variations/unmerge_d0_10/analyze.m @@ -1,30 +1,58 @@ -% Variation analysis -- the paper's linear mixed model on the successful-reach -% COUNT, fit on this folder's curated data subset. +% Variation analysis -- the paper's linear mixed model on this folder's data, +% for BOTH metrics: +% metric = count : behavior = # successes -> result.txt +% metric = rate : behavior = success / attempts -> result_rate.txt +% (rate uses only sessions with attempts > 0) % % model: behavior ~ stim + day + stim:day + (1|rat) -% behavior = successful reaches (count per session) -% stim = 1 for the treatment group(s), 0 for the control group(s) -% day = training day within this window (0 = first analyzed day) -% rat = subject (random intercept) -% -% Self-contained: reads data.csv beside this script and writes result.txt. -% Run headless from this folder with: matlab -batch "analyze" -% (This is a copy of analysis/matlab/variation_analyze.m; see make_variations.m.) +% stim = 1 treatment / 0 control; day = training day within window (0 = +% first analyzed day); rat = subject (random intercept). +% For the interaction we report residual DF, Satterthwaite DF, and the honest +% per-animal random-slope test. Self-contained: reads data.csv beside this +% script. Run headless with: matlab -batch "analyze" +% (Copy of analysis/matlab/variation_analyze.m; see make_variations.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end -vname = regexprep(here, '.*[/\\]', ''); % folder name = variation id +vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -tbl = table(D.success, D.day - min(D.day), double(D.stim), categorical(D.subject), ... +Rc = localAnalyze(D, 'count', here, vname); +Rr = localAnalyze(D, 'rate', here, vname); + +% Machine-readable handoff for SUMMARY.csv (count drives it; rate appended). +VARRESULT = struct('name', vname, 'nRats', Rc.nRats, 'nObs', Rc.nObs, ... + 'interP', Rc.interP, 'interEst', Rc.interEst, ... + 'interPsatt', Rc.interPsatt, 'interPrs', Rc.interPrs, ... + 'stimP', Rc.stimP, 'dayP', Rc.dayP, 'covEqual', Rc.covEqual, ... + 'interPrate', Rr.interP, 'interEstRate', Rr.interEst, 'interPrsRate', Rr.interPrs); + +% ------------------------------------------------------------------ helper +function R = localAnalyze(D, metric, here, vname) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); + beh = D.success ./ D.total; + mlabel = 'success RATE (success/attempts)'; suffix = '_rate'; +else + beh = D.success; + mlabel = 'success COUNT'; suffix = ''; +end +R = struct('interP', NaN, 'interEst', NaN, 'interPsatt', NaN, 'interPrs', NaN, ... + 'stimP', NaN, 'dayP', NaN, 'nRats', numel(unique(D.subject)), ... + 'nObs', height(D), 'covEqual', false); +if numel(unique(D.stim)) < 2 || numel(unique(D.day)) < 2 + localWrite(sprintf('VARIATION: %s [metric: %s]\nInsufficient data for this metric.\n', ... + vname, mlabel), here, suffix); + return +end + +tbl = table(beh, D.day - min(D.day), double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); m = fitlme(tbl, 'behavior ~ stim + day + stim:day + (1|rat)'); C = m.Coefficients; A = anova(m); ci = coefCI(m); -As = anova(m, 'DFMethod', 'satterthwaite'); % Satterthwaite denominator DF +As = anova(m, 'DFMethod', 'satterthwaite'); -% Honest test: refit with a per-animal random SLOPE so the interaction DF -% collapses toward the animal count (guarded -- may not converge in short windows). rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; wst = warning('off', 'all'); try @@ -43,6 +71,7 @@ row = @(nm, t) sprintf('%-26s t(%d)=%6.2f F(%d)=%7.3f p=%.4g p=%.4g (df=%.0f C.DF(gi(t)), C.tStat(gi(t)), A.DF1(ga(t)), A.FStat(ga(t)), C.pValue(gi(t)), ... As.pValue(gs(t)), As.DF2(gs(t))); +ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); maxT = max(D.day(D.stim == 1)); minT = min(D.day(D.stim == 1)); maxC = max(D.day(D.stim == 0)); minC = min(D.day(D.stim == 0)); if abs(maxT - maxC) > 2 @@ -50,20 +79,14 @@ if abs(maxT - maxC) > 2 else cov = '(equal day coverage over this window)'; end - -ii = gi('day:stim'); pI = C.pValue(ii); eI = C.Estimate(ii); -if pI >= 0.05 - verdict = 'n.s. -- slopes parallel (no differential learning rate)'; -elseif eI > 0 - verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; -else - verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; -end +if pI >= 0.05; verdict = 'n.s. -- slopes parallel (no differential learning rate)'; +elseif eI > 0; verdict = 'SIGNIFICANT positive -- treatment improves FASTER (benefit accumulates)'; +else; verdict = 'SIGNIFICANT negative -- treatment improves SLOWER (groups converge)'; end bar = repmat('=', 1, 78); raw = regexprep(evalc('disp(m)'), '', ''); -s = sprintf('%s\nVARIATION: %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success COUNT)\n')]; +s = sprintf('%s\nVARIATION: %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + day + stim:day + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('day = training day within window (0 = first analyzed day)\n')]; s = [s sprintf('treatment (stim=1): %s\n', strjoin(cellstr(unique(D.group(D.stim == 1))), ', '))]; s = [s sprintf('control (stim=0): %s\n', strjoin(cellstr(unique(D.group(D.stim == 0))), ', '))]; @@ -74,7 +97,7 @@ s = [s sprintf('%-26s %-18s %-12s %s\n%s\n', 'effect', 't(df) / F(df1)', 'p (res s = [s row('stim x day (interaction)', 'day:stim')]; s = [s row('day (learning)', 'day')]; s = [s row('stim (main, window start)', 'stim')]; -s = [s sprintf('interaction 95%% CI: [%+.2f, %+.2f]\n', ci(ii, 1), ci(ii, 2))]; +s = [s sprintf('interaction 95%% CI: [%+.4g, %+.4g]\n', ci(ii, 1), ci(ii, 2))]; if rsOk s = [s sprintf('HONEST LME (per-animal random slope, day|rat): interaction F(1,%.1f)=%.2f, p=%.4g\n', rsDf, rsF, rsP)]; else @@ -82,17 +105,18 @@ else end s = [s sprintf([' (Satterthwaite DF ~= residual on this random-intercept model; the random-slope\n' ... ' model above is the honest learning-rate test -- DF collapses toward the animal count.)\n'])]; -s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.2f)\n', verdict, pI, eI)]; -s = [s sprintf('Paper (N=24): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +s = [s sprintf('INTERPRETATION: stim x day interaction %s (p=%.4g, slope diff=%+.4g)\n', verdict, pI, eI)]; +s = [s sprintf('Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.\n')]; +localWrite(s, here, suffix); + +R = struct('interP', pI, 'interEst', eI, 'interPsatt', As.pValue(gs('day:stim')), ... + 'interPrs', rsP, 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... + 'nRats', numel(unique(D.subject)), 'nObs', height(D), 'covEqual', abs(maxT - maxC) <= 2); +end + +function localWrite(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'result.txt'), 'w'); -fprintf(fid, '%s', s); -fclose(fid); - -% Machine-readable handoff for the summary table (see make_variations.m). -VARRESULT = struct('name', vname, 'nRats', numel(unique(D.subject)), ... - 'nObs', height(D), 'interP', pI, 'interEst', eI, ... - 'interPsatt', As.pValue(gs('day:stim')), 'interPrs', rsP, ... - 'stimP', C.pValue(gi('stim')), 'dayP', C.pValue(gi('day')), ... - 'covEqual', abs(maxT - maxC) <= 2); +fid = fopen(fullfile(here, ['result' suffix '.txt']), 'w'); +fprintf(fid, '%s', s); fclose(fid); +end diff --git a/analysis/matlab/variations/unmerge_d0_10/learning_curve.png b/analysis/matlab/variations/unmerge_d0_10/learning_curve.png index 0d415cc..5089669 100644 Binary files a/analysis/matlab/variations/unmerge_d0_10/learning_curve.png and b/analysis/matlab/variations/unmerge_d0_10/learning_curve.png differ diff --git a/analysis/matlab/variations/unmerge_d0_10/learning_curve_rate.png b/analysis/matlab/variations/unmerge_d0_10/learning_curve_rate.png new file mode 100644 index 0000000..0b30787 Binary files /dev/null and b/analysis/matlab/variations/unmerge_d0_10/learning_curve_rate.png differ diff --git a/analysis/matlab/variations/unmerge_d0_10/logpower_result.txt b/analysis/matlab/variations/unmerge_d0_10/logpower_result.txt index 85ffb18..330d817 100644 --- a/analysis/matlab/variations/unmerge_d0_10/logpower_result.txt +++ b/analysis/matlab/variations/unmerge_d0_10/logpower_result.txt @@ -1,7 +1,7 @@ ============================================================================== -LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_d0_10 +LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_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=3 nrep=120, alpha=0.05 @@ -10,9 +10,9 @@ stim x log(day) interaction: F(1,57)=6.230 p(resid)=0.01548 p(Satt)=0.01528 (d honest per-animal random slope (log-day): F(1,15.0)=5.28 p=0.0363 Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 small, .25 medium, .40 large) ---- power simulation (log-day ground truth: stim:logday=+11.60, ratSD=0.00, resSD=12.93) --- +--- power simulation (log-day ground truth: stim:logday=+11.6, ratSD=0, resSD=12.93) --- - true stim:log(day) = +11.60 (100% of observed) + true stim:log(day) = +11.6 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.34 | 0.69 <- observed @@ -22,7 +22,7 @@ Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 smal 16 | 1.00 | 1.00 24 | 1.00 | 1.00 - true stim:log(day) = +5.80 (50% of observed) + true stim:log(day) = +5.801 (50% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.10 | 0.25 <- observed @@ -32,5 +32,4 @@ Cohen's f (interaction, partial eta^2=0.017) = 0.131 (small-medium; f: .10 smal 16 | 0.90 | 0.90 24 | 0.97 | 0.97 -Read the per-animal column as the honest power; the LME column matches the -paper's power code (anova interaction p, observation-level DF) and is optimistic. +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/unmerge_d0_10/logpower_result_rate.txt b/analysis/matlab/variations/unmerge_d0_10/logpower_result_rate.txt new file mode 100644 index 0000000..c75ff5b --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_10/logpower_result_rate.txt @@ -0,0 +1,35 @@ +============================================================================== +LOG-DAY MODEL + COHEN'S f + POWER -- unmerge_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=3 nrep=120, alpha=0.05 + +--- fitted on real data --- +stim x log(day) interaction: F(1,57)=9.174 p(resid)=0.003684 p(Satt)=0.003695 (df=57) + honest per-animal random slope (log-day): F(1,35.1)=8.84 p=0.005302 +Cohen's f (interaction, partial eta^2=0.030) = 0.177 (small-medium; f: .10 small, .25 medium, .40 large) + +--- power simulation (log-day ground truth: stim:logday=+0.08869, ratSD=0.04059, resSD=0.08029) --- + + true stim:log(day) = +0.08869 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.54 | 0.85 <- observed + 5 | 0.95 | 0.98 + 8 | 1.00 | 1.00 + 12 | 1.00 | 1.00 + 16 | 1.00 | 1.00 + 24 | 1.00 | 1.00 + + true stim:log(day) = +0.04434 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.17 | 0.33 <- observed + 5 | 0.40 | 0.52 + 8 | 0.69 | 0.68 + 12 | 0.89 | 0.90 + 16 | 0.95 | 0.96 + 24 | 0.99 | 0.99 + +Read per-animal as the honest power; LME matches the paper's power code (optimistic). diff --git a/analysis/matlab/variations/unmerge_d0_10/logpowersim.m b/analysis/matlab/variations/unmerge_d0_10/logpowersim.m index e0686ee..5e38720 100644 --- a/analysis/matlab/variations/unmerge_d0_10/logpowersim.m +++ b/analysis/matlab/variations/unmerge_d0_10/logpowersim.m @@ -1,50 +1,52 @@ -% Variation log-day analysis + Cohen's f + power simulation. +% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics: +% metric = count : behavior = # successes -> logpower_result.txt +% metric = rate : behavior = success / attempts -> logpower_result_rate.txt % -% The paper's power code models behavior against LOG training day, not raw day: +% The paper's power code models behavior against LOG training day: % behavior ~ stim + log(day) + stim:log(day) + (1|rat). -% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper -% "Day 1"), so log(day + 1) reproduces their transform exactly. -% -% This script (a) refits that log-day model on data.csv, (b) reports the -% interaction (residual DF, Satterthwaite DF, and the honest per-animal -% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the -% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior) -% -- and (c) runs the Monte-Carlo power simulation on the log-day model, -% scoring per-animal (cluster-honest) and LME power across N. -% Writes logpower_result.txt. Run: matlab -batch "logpowersim" +% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1) +% reproduces their transform (our day 0 = paper "Day 1"). For each metric this +% refits that model, reports the interaction (residual / Satterthwaite / honest +% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then +% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power). +% Run: matlab -batch "logpowersim" % (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1) -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; -warnState = warning('off', 'all'); -rng(1); +localLogPower(D, 'count', here, vname); +localLogPower(D, 'rate', here, vname); -logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day) -tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ... +% ---------------------------------------------------------------- per metric +function localLogPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; +FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1) +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); +logday = log(D.day + 1); +tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')]; +s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)]; s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for this analysis.\n')]; + localFinish(s, here, suffix); warning(warnState); return end -% ---- fitted model on the real data ---- full = fitlme(tbl0, FORMULA); An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite'); ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim'); @@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)'); eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0); cohenf = sqrt(eta2part / (1 - eta2part)); -% honest per-animal random-slope interaction rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false; try mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)'); @@ -76,19 +77,18 @@ end s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ... eta2part, cohenf, mag)]; -% ---- power simulation on the log-day ground truth ---- cn = full.CoefficientNames; be = full.fixedEffects; b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim')); bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE); -days = unique(tbl0.day); % the log(day) grid +days = unique(tbl0.day); -s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ... +s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ... bInt, sRat, sRes)]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -102,19 +102,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end -s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ... - 'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'logpower_result.txt'), 'w'); +fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/unmerge_d0_10/plotcurve.m b/analysis/matlab/variations/unmerge_d0_10/plotcurve.m index bfc1386..193a44b 100644 --- a/analysis/matlab/variations/unmerge_d0_10/plotcurve.m +++ b/analysis/matlab/variations/unmerge_d0_10/plotcurve.m @@ -1,11 +1,12 @@ -% Variation learning-curve plot, in the style of the paper: +% Variation learning-curve plots, in the style of the paper: % "Lines indicate mean (and SEM) across animals in the anodal (red) and % control (blue) groups." -% Plots mean +/- SEM successful reaches per training day for the treatment / -% anodal group (stim = 1, red) and the control group (stim = 0, blue), reading -% this folder's data.csv and saving learning_curve.png. The per-group N is read -% from the data (each variation pools different groups), so the legend shows the -% actual counts. Training day is 1-indexed (our day 0 = the paper's "Day 1"). +% Produces TWO figures from this folder's data.csv: +% learning_curve.png # successes (count) per training day +% learning_curve_rate.png success rate (success/attempts) per training day +% anodal / treatment = stim 1 (red); control = stim 0 (blue). Per-group N is +% read from the data. Training day is 1-indexed (our day 0 = paper "Day 1") and +% the x-axis tick labels are drawn vertically. % Run: matlab -batch "plotcurve" % (Copy of analysis/matlab/variation_plot.m; see make_variation_plot.m.) @@ -14,15 +15,20 @@ if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); -days = unique(D.day); % 0-indexed -xd = days + 1; % plot as 1-indexed training day (paper axis) - +days = unique(D.day); +xd = days + 1; % plot as 1-indexed training day (paper axis) red = [0.85 0.10 0.10]; blue = [0.10 0.30 0.85]; -[Ma, Sa, na] = localCurve(D, 1, days); % anodal / treatment (stim = 1) -[Mc, Sc, nc] = localCurve(D, 0, days); % control (stim = 0) +localPlot(D, days, xd, 'count', '# successes', ... + fullfile(here, 'learning_curve.png'), vname, red, blue); +localPlot(D, days, xd, 'rate', 'success rate', ... + fullfile(here, 'learning_curve_rate.png'), vname, red, blue); +% ------------------------------------------------------------------ helpers +function localPlot(D, days, xd, metric, ylab, outFile, vname, red, blue) +[Ma, Sa, na] = localCurve(D, 1, days, metric); % anodal / treatment +[Mc, Sc, nc] = localCurve(D, 0, days, metric); % control fig = figure('Visible', 'off', 'Color', 'w', 'Position', [100 100 560 460]); hold on e1 = errorbar(xd, Ma, Sa, '-o', 'Color', red, 'MarkerFaceColor', red, 'LineWidth', 2); @@ -30,26 +36,30 @@ e2 = errorbar(xd, Mc, Sc, '-o', 'Color', blue, 'MarkerFaceColor', blue, 'LineWid hold off legend([e1 e2], {sprintf('anodal, N = %d', na), sprintf('control, N = %d', nc)}, ... 'Location', 'northwest', 'Box', 'off'); -xlabel('training day'); -ylabel('# successes'); +xlabel('training day'); ylabel(ylab); title(vname, 'Interpreter', 'none'); set(gca, 'XTick', xd, 'FontName', 'Arial', 'FontSize', 13, 'LineWidth', 1.5, 'Box', 'off'); - -outFile = fullfile(here, 'learning_curve.png'); +xtickangle(90); % vertical x-axis tick labels exportgraphics(fig, outFile, 'Resolution', 150); close(fig); -fprintf('%s: wrote learning_curve.png (anodal N=%d, control N=%d)\n', vname, na, nc); +fprintf('%s: wrote %s (anodal N=%d, control N=%d)\n', vname, outFile, na, nc); +end -% ------------------------------------------------------------------ helper -function [M, S, n] = localCurve(D, stimVal, days) -%LOCALCURVE Per-day mean and SEM of successes across the animals in a group. +function [M, S, n] = localCurve(D, stimVal, days, metric) +%LOCALCURVE Per-day mean and SEM across the animals in a group, for a metric. subs = unique(D.subject(D.stim == stimVal)); n = numel(subs); X = nan(numel(days), n); for j = 1:n for i = 1:numel(days) r = D.subject == subs(j) & D.day == days(i); - if any(r); X(i, j) = mean(D.success(r)); end + if ~any(r); continue; end + if strcmp(metric, 'rate') + tot = sum(D.total(r)); + if tot > 0; X(i, j) = sum(D.success(r)) / tot; end + else + X(i, j) = mean(D.success(r)); + end end end M = mean(X, 2, 'omitnan'); diff --git a/analysis/matlab/variations/unmerge_d0_10/power_result.txt b/analysis/matlab/variations/unmerge_d0_10/power_result.txt index dfc1f55..91a393c 100644 --- a/analysis/matlab/variations/unmerge_d0_10/power_result.txt +++ b/analysis/matlab/variations/unmerge_d0_10/power_result.txt @@ -1,11 +1,11 @@ ============================================================================== -POWER SIMULATION -- unmerge_d0_10 +POWER SIMULATION -- unmerge_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=3 nrep=120, alpha=0.05 -ground truth: stim:day=+1.56/day, rat SD=5.35, residual SD=12.51, days=11 +ground truth: stim:day=+1.558/day, rat SD=5.354, residual SD=12.51, days=11 - true stim:day interaction = +1.56 (100% of observed) + true stim:day interaction = +1.558 (100% of observed) N/group | per-animal power | LME power ------------------------------------------ 3 | 0.15 | 0.31 <- observed @@ -15,7 +15,7 @@ ground truth: stim:day=+1.56/day, rat SD=5.35, residual SD=12.51, days=11 16 | 0.93 | 0.95 24 | 0.99 | 1.00 - true stim:day interaction = +0.78 (50% of observed) + true stim:day interaction = +0.7788 (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.56/day, rat SD=5.35, residual SD=12.51, days=11 16 | 0.47 | 0.47 24 | 0.63 | 0.60 -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). diff --git a/analysis/matlab/variations/unmerge_d0_10/power_result_rate.txt b/analysis/matlab/variations/unmerge_d0_10/power_result_rate.txt new file mode 100644 index 0000000..ce7f9f3 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_10/power_result_rate.txt @@ -0,0 +1,28 @@ +============================================================================== +POWER SIMULATION -- unmerge_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=3 nrep=120, alpha=0.05 +ground truth: stim:day=+0.01218/day, rat SD=0.04284, residual SD=0.08704, days=11 + + true stim:day interaction = +0.01218 (100% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.20 | 0.36 <- observed + 5 | 0.51 | 0.64 + 8 | 0.77 | 0.79 + 12 | 0.96 | 0.97 + 16 | 0.97 | 0.97 + 24 | 1.00 | 1.00 + + true stim:day interaction = +0.006089 (50% of observed) + N/group | per-animal power | LME power + ------------------------------------------ + 3 | 0.12 | 0.11 <- observed + 5 | 0.11 | 0.18 + 8 | 0.28 | 0.31 + 12 | 0.36 | 0.41 + 16 | 0.58 | 0.62 + 24 | 0.71 | 0.72 + +Read the per-animal column as the honest power; LME is optimistic (obs-level DF). diff --git a/analysis/matlab/variations/unmerge_d0_10/powersim.m b/analysis/matlab/variations/unmerge_d0_10/powersim.m index 554d032..ecd9c6c 100644 --- a/analysis/matlab/variations/unmerge_d0_10/powersim.m +++ b/analysis/matlab/variations/unmerge_d0_10/powersim.m @@ -1,45 +1,49 @@ % Variation power simulation -- Monte-Carlo power for the paper's stim x day -% interaction, using THIS folder's data as the ground truth. +% interaction, using THIS folder's data as the ground truth, for BOTH metrics: +% metric = count : behavior = # successes -> power_result.txt +% metric = rate : behavior = success / attempts -> power_result_rate.txt % % Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv -% (success COUNT; day within-window). Its fixed effects, per-rat intercept SD, -% and residual SD generate NREP synthetic datasets at each rats-per-group N and -% each true-effect multiplier (1 = observed slope, 0.5 = half). Each dataset is -% scored at alpha = 0.05 two ways: -% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest -- the -% honest power, matching the random-slope / per-animal inference) +% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD +% generate NREP synthetic datasets at each rats-per-group N and each true-effect +% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by: +% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power) % LME : the fitlme stim:day p (observation-level DF -- optimistic) -% Writes power_result.txt beside this script. Run: matlab -batch "powersim" +% Writes power_result[_rate].txt. Run: matlab -batch "powersim" % (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.) here = fileparts(mfilename('fullpath')); if isempty(here); here = pwd; end vname = regexprep(here, '.*[/\\]', ''); - D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string'); + +localPower(D, 'count', here, vname); +localPower(D, 'rate', here, vname); + +% ---------------------------------------------------------------- per metric +function localPower(D, metric, here, vname) +NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120; FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; -NS = [3 5 8 12 16 24]; -EFFMULS = [1 0.5]; -NREP = 120; - -warnState = warning('off', 'all'); -rng(1); - +if strcmp(metric, 'rate') + D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate'; +else + beh = D.success; mlabel = '# successes (count)'; suffix = ''; +end +warnState = warning('off', 'all'); rng(1); day0 = min(D.day); -tbl0 = table(D.success, D.day - day0, double(D.stim), categorical(D.subject), ... +tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ... 'VariableNames', {'behavior', 'day', 'stim', 'rat'}); - nStim = numel(unique(D.subject(D.stim == 1))); nCtrl = numel(unique(D.subject(D.stim == 0))); bar = repmat('=', 1, 78); -s = sprintf('%s\nPOWER SIMULATION -- %s\n%s\n', bar, vname, bar); -s = [s sprintf('model: %s (success COUNT; day within-window)\n', FORMULA)]; +s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar); +s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)]; s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)]; if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2 - s = [s sprintf('\nInsufficient data for a power simulation (need >=2 animals/group and >=2 days).\n')]; - localFinish(s, here); warning(warnState); return + s = [s sprintf('\nInsufficient data for a power simulation.\n')]; + localFinish(s, here, suffix); warning(warnState); return end lme = fitlme(tbl0, FORMULA); @@ -49,14 +53,13 @@ bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim')); psi = covarianceParameters(lme); sRat = sqrt(psi{1}); sRes = sqrt(lme.MSE); days = (0:max(tbl0.day))'; -s = [s sprintf('ground truth: stim:day=%+.2f/day, rat SD=%.2f, residual SD=%.2f, days=%d\n', ... +s = [s sprintf('ground truth: stim:day=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ... bInt, sRat, sRes, numel(days))]; for eMul = EFFMULS bI = bInt * eMul; - s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok - s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok - s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)]; + s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; for N = NS sigPA = 0; sigL = 0; for r = 1:NREP @@ -70,20 +73,18 @@ for eMul = EFFMULS end star = ''; if N == nStim || N == nCtrl; star = ' <- observed'; end - s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok + s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; end end - -s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ... - 'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])]; - -localFinish(s, here); +s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')]; +localFinish(s, here, suffix); warning(warnState); +end % ---------------------------------------------------------------- helpers -function localFinish(s, here) +function localFinish(s, here, suffix) fprintf('%s', s); -fid = fopen(fullfile(here, 'power_result.txt'), 'w'); +fid = fopen(fullfile(here, ['power_result' suffix '.txt']), 'w'); fprintf(fid, '%s', s); fclose(fid); end diff --git a/analysis/matlab/variations/unmerge_d0_10/result.txt b/analysis/matlab/variations/unmerge_d0_10/result.txt index d9b1736..6c8b6c9 100644 --- a/analysis/matlab/variations/unmerge_d0_10/result.txt +++ b/analysis/matlab/variations/unmerge_d0_10/result.txt @@ -1,5 +1,5 @@ ============================================================================== -VARIATION: unmerge_d0_10 +VARIATION: unmerge_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(57)= 1.49 F(1)= 2.208 p=0.1428 p=0.1428 (df=57) day (learning) t(57)= 10.08 F(1)=101.556 p=2.832e-14 p=1.931e-14 (df=59) stim (main, window start) t(57)= 0.68 F(1)= 0.466 p=0.4974 p=0.5037 (df=17) -interaction 95% CI: [-0.54, +3.66] +interaction 95% CI: [-0.5416, +3.657] HONEST LME (per-animal random slope, day|rat): interaction F(1,13.0)=1.46, p=0.2479 (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.1428, slope diff=+1.56) -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.1428, slope diff=+1.558) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008. diff --git a/analysis/matlab/variations/unmerge_d0_10/result_rate.txt b/analysis/matlab/variations/unmerge_d0_10/result_rate.txt new file mode 100644 index 0000000..95601b0 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_10/result_rate.txt @@ -0,0 +1,68 @@ +============================================================================== +VARIATION: unmerge_d0_10 [metric: success RATE (success/attempts)] +============================================================================== +model: behavior ~ stim + day + stim:day + (1|rat) (behavior = success RATE (success/attempts)) +day = training day within window (0 = first analyzed day) +treatment (stim=1): Electrode-Box-B2 +control (stim=0): Electrode-Box-A2 +N = 6 rats, 61 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 61 + 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 + -105.35 -92.689 58.677 -117.35 + +Fixed effects coefficients (95% CIs): + Name Estimate SE tStat DF pValue + {'(Intercept)'} 0.22998 0.038137 6.0304 57 1.2826e-07 + {'day' } 0.040747 0.0055227 7.378 57 7.4583e-10 + {'stim' } 0.02493 0.053569 0.46538 57 0.64343 + {'day:stim' } 0.012178 0.0073113 1.6656 57 0.10127 + + + Lower Upper + 0.15361 0.30635 + 0.029687 0.051806 + -0.082339 0.1322 + -0.0024627 0.026819 + +Random effects covariance parameters (95% CIs): +Group: rat (6 Levels) + Name1 Name2 Type Estimate + {'(Intercept)'} {'(Intercept)'} {'std'} 0.04284 + + + Lower Upper + 0.01883 0.097462 + +Group: Error + Name Estimate Lower Upper + {'Res Std'} 0.087035 0.072173 0.10496 + + +effect t(df) / F(df1) p (resid) Satterthwaite: p (df) +---------------------------------------------------------------------------- +stim x day (interaction) t(57)= 1.67 F(1)= 2.774 p=0.1013 p=0.1013 (df=57) +day (learning) t(57)= 7.38 F(1)= 54.435 p=7.458e-10 p=6.427e-10 (df=59) +stim (main, window start) t(57)= 0.47 F(1)= 0.217 p=0.6434 p=0.6483 (df=15) +interaction 95% CI: [-0.002463, +0.02682] +HONEST LME (per-animal random slope, day|rat): interaction F(1,40.5)=2.61, p=0.1136 + (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.1013, slope diff=+0.01218) +Paper (N=24, count): interaction t(227)=2.68, F(1)=7.12, p=0.008.