analysis(matlab): rate + count outputs, variation batch 5
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
@@ -1,50 +1,52 @@
|
||||
% Variation log-day analysis + Cohen's f + power simulation.
|
||||
% Variation log-day analysis + Cohen's f + power simulation, for BOTH metrics:
|
||||
% metric = count : behavior = # successes -> logpower_result.txt
|
||||
% metric = rate : behavior = success / attempts -> logpower_result_rate.txt
|
||||
%
|
||||
% The paper's power code models behavior against LOG training day, not raw day:
|
||||
% The paper's power code models behavior against LOG training day:
|
||||
% behavior ~ stim + log(day) + stim:log(day) + (1|rat).
|
||||
% Their day is 1-indexed (1..10); our data.csv day is 0-indexed (day 0 = paper
|
||||
% "Day 1"), so log(day + 1) reproduces their transform exactly.
|
||||
%
|
||||
% This script (a) refits that log-day model on data.csv, (b) reports the
|
||||
% interaction (residual DF, Satterthwaite DF, and the honest per-animal
|
||||
% random-slope test) and Cohen's f -- the partial-eta^2 effect size of the
|
||||
% interaction, var(fitted_full) - var(fitted_no_interaction) over var(behavior)
|
||||
% -- and (c) runs the Monte-Carlo power simulation on the log-day model,
|
||||
% scoring per-animal (cluster-honest) and LME power across N.
|
||||
% Writes logpower_result.txt. Run: matlab -batch "logpowersim"
|
||||
% Their day is 1-indexed; our data.csv day is 0-indexed, so log(day + 1)
|
||||
% reproduces their transform (our day 0 = paper "Day 1"). For each metric this
|
||||
% refits that model, reports the interaction (residual / Satterthwaite / honest
|
||||
% per-animal random-slope DF) and Cohen's f (partial-eta^2 effect size), then
|
||||
% runs the Monte-Carlo power sim (per-animal cluster-honest + LME power).
|
||||
% Run: matlab -batch "logpowersim"
|
||||
% (Copy of analysis/matlab/variation_logpower.m; see make_variation_logpower.m.)
|
||||
|
||||
here = fileparts(mfilename('fullpath'));
|
||||
if isempty(here); here = pwd; end
|
||||
vname = regexprep(here, '.*[/\\]', '');
|
||||
|
||||
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' column = log(day+1)
|
||||
NS = [3 5 8 12 16 24];
|
||||
EFFMULS = [1 0.5];
|
||||
NREP = 120;
|
||||
|
||||
warnState = warning('off', 'all');
|
||||
rng(1);
|
||||
localLogPower(D, 'count', here, vname);
|
||||
localLogPower(D, 'rate', here, vname);
|
||||
|
||||
logday = log(D.day + 1); % 0-indexed day -> their log(1-indexed day)
|
||||
tbl0 = table(D.success, logday, double(D.stim), categorical(D.subject), ...
|
||||
% ---------------------------------------------------------------- per metric
|
||||
function localLogPower(D, metric, here, vname)
|
||||
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
||||
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)'; % 'day' = log(day+1)
|
||||
if strcmp(metric, 'rate')
|
||||
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
||||
else
|
||||
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
||||
end
|
||||
warnState = warning('off', 'all'); rng(1);
|
||||
logday = log(D.day + 1);
|
||||
tbl0 = table(beh, logday, double(D.stim), categorical(D.subject), ...
|
||||
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
||||
nStim = numel(unique(D.subject(D.stim == 1)));
|
||||
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
||||
|
||||
bar = repmat('=', 1, 78);
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s\n%s\n', bar, vname, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (success COUNT)\n')];
|
||||
s = sprintf('%s\nLOG-DAY MODEL + COHEN''S f + POWER -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
||||
s = [s sprintf('model: behavior ~ stim + log(day) + stim:log(day) + (1|rat) (behavior = %s)\n', mlabel)];
|
||||
s = [s sprintf('log(day) uses 1-indexed training day (our day 0 = paper "Day 1")\n')];
|
||||
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
||||
|
||||
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
||||
s = [s sprintf('\nInsufficient data for this analysis (need >=2 animals/group and >=2 days).\n')];
|
||||
localFinish(s, here); warning(warnState); return
|
||||
s = [s sprintf('\nInsufficient data for this analysis.\n')];
|
||||
localFinish(s, here, suffix); warning(warnState); return
|
||||
end
|
||||
|
||||
% ---- fitted model on the real data ----
|
||||
full = fitlme(tbl0, FORMULA);
|
||||
An = anova(full); Asatt = anova(full, 'DFMethod', 'satterthwaite');
|
||||
ii = strcmp(An.Term, 'day:stim'); is = strcmp(Asatt.Term, 'day:stim');
|
||||
@@ -52,7 +54,6 @@ reduced = fitlme(tbl0, 'behavior ~ stim + day + (1|rat)');
|
||||
eta2part = max((var(fitted(full)) - var(fitted(reduced))) / var(tbl0.behavior), 0);
|
||||
cohenf = sqrt(eta2part / (1 - eta2part));
|
||||
|
||||
% honest per-animal random-slope interaction
|
||||
rsP = NaN; rsDf = NaN; rsF = NaN; rsOk = false;
|
||||
try
|
||||
mr = fitlme(tbl0, 'behavior ~ stim + day + stim:day + (day|rat)');
|
||||
@@ -76,19 +77,18 @@ end
|
||||
s = [s sprintf('Cohen''s f (interaction, partial eta^2=%.3f) = %.3f (%s; f: .10 small, .25 medium, .40 large)\n', ...
|
||||
eta2part, cohenf, mag)];
|
||||
|
||||
% ---- power simulation on the log-day ground truth ----
|
||||
cn = full.CoefficientNames; be = full.fixedEffects;
|
||||
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||
bDay = be(strcmp(cn, 'day')); bInt = be(strcmp(cn, 'day:stim'));
|
||||
psi = covarianceParameters(full); sRat = sqrt(psi{1}); sRes = sqrt(full.MSE);
|
||||
days = unique(tbl0.day); % the log(day) grid
|
||||
days = unique(tbl0.day);
|
||||
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.2f, ratSD=%.2f, resSD=%.2f) ---\n', ...
|
||||
s = [s sprintf('\n--- power simulation (log-day ground truth: stim:logday=%+.4g, ratSD=%.4g, resSD=%.4g) ---\n', ...
|
||||
bInt, sRat, sRes)];
|
||||
for eMul = EFFMULS
|
||||
bI = bInt * eMul;
|
||||
s = [s sprintf('\n true stim:log(day) = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||
s = [s sprintf('\n true stim:log(day) = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
||||
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
||||
for N = NS
|
||||
sigPA = 0; sigL = 0;
|
||||
for r = 1:NREP
|
||||
@@ -102,19 +102,18 @@ for eMul = EFFMULS
|
||||
end
|
||||
star = '';
|
||||
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)];
|
||||
end
|
||||
end
|
||||
s = [s sprintf(['\nRead the per-animal column as the honest power; the LME column matches the\n' ...
|
||||
'paper''s power code (anova interaction p, observation-level DF) and is optimistic.\n'])];
|
||||
|
||||
localFinish(s, here);
|
||||
s = [s sprintf('\nRead per-animal as the honest power; LME matches the paper''s power code (optimistic).\n')];
|
||||
localFinish(s, here, suffix);
|
||||
warning(warnState);
|
||||
end
|
||||
|
||||
% ---------------------------------------------------------------- helpers
|
||||
function localFinish(s, here)
|
||||
function localFinish(s, here, suffix)
|
||||
fprintf('%s', s);
|
||||
fid = fopen(fullfile(here, 'logpower_result.txt'), 'w');
|
||||
fid = fopen(fullfile(here, ['logpower_result' suffix '.txt']), 'w');
|
||||
fprintf(fid, '%s', s); fclose(fid);
|
||||
end
|
||||
|
||||
|
||||
Reference in New Issue
Block a user