analysis(matlab): add per-variation power simulation (powersim.m + power_result.txt)
Add variation_power.m (self-contained Monte-Carlo power for the paper's stim x day interaction, using each folder's own data as ground truth; scores per-animal cluster-honest power + LME power across N=[3..24] and effect multipliers 1/0.5) and make_variation_power.m, which drops powersim.m into every variations/<name>/ folder and runs it, writing power_result.txt beside the existing data.csv/analyze.m/result.txt. Named powersim (not power) to avoid shadowing the MATLAB builtin. All 28 folders processed. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Experiments DB Dev
@@ -0,0 +1,29 @@
|
|||||||
|
function make_variation_power()
|
||||||
|
%MAKE_VARIATION_POWER Add a power-simulation script + result to every variation.
|
||||||
|
% MAKE_VARIATION_POWER() copies variation_power.m into each
|
||||||
|
% variations/<name>/ folder that contains a data.csv (as powersim.m) and runs
|
||||||
|
% it, producing power_result.txt beside the existing data.csv / analyze.m /
|
||||||
|
% result.txt. Re-runnable; also picks up any folders added later.
|
||||||
|
|
||||||
|
thisDir = fileparts(mfilename('fullpath'));
|
||||||
|
template = fullfile(thisDir, 'variation_power.m');
|
||||||
|
root = fullfile(thisDir, 'variations');
|
||||||
|
|
||||||
|
d = dir(root);
|
||||||
|
n = 0;
|
||||||
|
for i = 1:numel(d)
|
||||||
|
if ~d(i).isdir || ismember(d(i).name, {'.', '..'}); continue; end
|
||||||
|
folder = fullfile(root, d(i).name);
|
||||||
|
if ~exist(fullfile(folder, 'data.csv'), 'file'); continue; end
|
||||||
|
copyfile(template, fullfile(folder, 'powersim.m'));
|
||||||
|
fprintf('\n### %s ###\n', d(i).name);
|
||||||
|
localRun(fullfile(folder, 'powersim.m'));
|
||||||
|
n = n + 1;
|
||||||
|
end
|
||||||
|
fprintf('\nAdded powersim.m + power_result.txt to %d variation folders.\n', n);
|
||||||
|
|
||||||
|
end
|
||||||
|
|
||||||
|
function localRun(scriptPath)
|
||||||
|
run(scriptPath);
|
||||||
|
end
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_a2_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.20/day, rat SD=0.00, residual SD=13.47, days=11
|
||||||
|
|
||||||
|
true stim:day interaction = +1.20 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.17 <- observed
|
||||||
|
5 | 0.20 | 0.24
|
||||||
|
8 | 0.35 | 0.43
|
||||||
|
12 | 0.62 | 0.67
|
||||||
|
16 | 0.72 | 0.72
|
||||||
|
24 | 0.88 | 0.88
|
||||||
|
|
||||||
|
true stim:day interaction = +0.60 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.08 <- observed
|
||||||
|
5 | 0.05 | 0.07
|
||||||
|
8 | 0.10 | 0.11
|
||||||
|
12 | 0.17 | 0.18
|
||||||
|
16 | 0.29 | 0.30
|
||||||
|
24 | 0.31 | 0.33
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_a2_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=4 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+4.68/day, rat SD=0.00, residual SD=12.01, days=6
|
||||||
|
|
||||||
|
true stim:day interaction = +4.68 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.23 | 0.54 <- observed
|
||||||
|
5 | 0.57 | 0.70
|
||||||
|
8 | 0.85 | 0.88
|
||||||
|
12 | 0.98 | 0.99
|
||||||
|
16 | 0.99 | 0.99
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +2.34 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.08 | 0.20 <- observed
|
||||||
|
5 | 0.17 | 0.23
|
||||||
|
8 | 0.26 | 0.32
|
||||||
|
12 | 0.44 | 0.52
|
||||||
|
16 | 0.56 | 0.60
|
||||||
|
24 | 0.82 | 0.84
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_a2_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.00/day, rat SD=6.34, residual SD=8.21, days=5
|
||||||
|
|
||||||
|
true stim:day interaction = +1.00 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.03 | 0.05 <- observed
|
||||||
|
5 | 0.09 | 0.13
|
||||||
|
8 | 0.10 | 0.11
|
||||||
|
12 | 0.12 | 0.14
|
||||||
|
16 | 0.12 | 0.14
|
||||||
|
24 | 0.27 | 0.32
|
||||||
|
|
||||||
|
true stim:day interaction = +0.50 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.05 | 0.09 <- observed
|
||||||
|
5 | 0.07 | 0.07
|
||||||
|
8 | 0.05 | 0.07
|
||||||
|
12 | 0.09 | 0.09
|
||||||
|
16 | 0.08 | 0.11
|
||||||
|
24 | 0.12 | 0.12
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_b2_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.60/day, rat SD=6.05, residual SD=12.42, days=11
|
||||||
|
|
||||||
|
true stim:day interaction = +1.60 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.17 | 0.33 <- observed
|
||||||
|
5 | 0.46 | 0.58
|
||||||
|
8 | 0.67 | 0.70
|
||||||
|
12 | 0.90 | 0.93
|
||||||
|
16 | 0.95 | 0.97
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.80 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.11 <- observed
|
||||||
|
5 | 0.09 | 0.14
|
||||||
|
8 | 0.23 | 0.25
|
||||||
|
12 | 0.33 | 0.35
|
||||||
|
16 | 0.51 | 0.53
|
||||||
|
24 | 0.67 | 0.66
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_b2_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=4, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+4.71/day, rat SD=5.87, residual SD=11.30, days=6
|
||||||
|
|
||||||
|
true stim:day interaction = +4.71 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.27 | 0.59 <- observed
|
||||||
|
5 | 0.63 | 0.72
|
||||||
|
8 | 0.91 | 0.93
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 0.99 | 0.99
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +2.35 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.10 | 0.20 <- observed
|
||||||
|
5 | 0.20 | 0.24
|
||||||
|
8 | 0.31 | 0.33
|
||||||
|
12 | 0.52 | 0.54
|
||||||
|
16 | 0.62 | 0.62
|
||||||
|
24 | 0.85 | 0.87
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- boxa_b2_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=4, control n=2 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.90/day, rat SD=6.23, residual SD=8.12, days=5
|
||||||
|
|
||||||
|
true stim:day interaction = +1.90 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.04 | 0.17
|
||||||
|
5 | 0.17 | 0.24
|
||||||
|
8 | 0.24 | 0.33
|
||||||
|
12 | 0.38 | 0.36
|
||||||
|
16 | 0.51 | 0.55
|
||||||
|
24 | 0.68 | 0.71
|
||||||
|
|
||||||
|
true stim:day interaction = +0.95 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.11
|
||||||
|
5 | 0.09 | 0.12
|
||||||
|
8 | 0.06 | 0.12
|
||||||
|
12 | 0.14 | 0.14
|
||||||
|
16 | 0.18 | 0.19
|
||||||
|
24 | 0.25 | 0.28
|
||||||
|
|
||||||
|
Read the per-animal column as the honest power. At the observed N this study
|
||||||
|
is typically underpowered; per-animal power reaches ~0.8 only at larger N.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- matched_current_d0_3
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+10.20/day, rat SD=2.92, residual SD=8.54, days=4
|
||||||
|
|
||||||
|
true stim:day interaction = +10.20 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.55 | 0.91 <- observed
|
||||||
|
5 | 0.93 | 0.98
|
||||||
|
8 | 1.00 | 1.00
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +5.10 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.33 <- observed
|
||||||
|
5 | 0.42 | 0.57
|
||||||
|
8 | 0.67 | 0.76
|
||||||
|
12 | 0.89 | 0.90
|
||||||
|
16 | 0.97 | 0.97
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_a2_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.43/day, rat SD=8.75, residual SD=13.78, days=11
|
||||||
|
|
||||||
|
true stim:day interaction = +1.43 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.09 | 0.23 <- observed
|
||||||
|
5 | 0.33 | 0.38
|
||||||
|
8 | 0.48 | 0.53
|
||||||
|
12 | 0.74 | 0.78
|
||||||
|
16 | 0.84 | 0.84
|
||||||
|
24 | 0.95 | 0.94
|
||||||
|
|
||||||
|
true stim:day interaction = +0.72 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.09 <- observed
|
||||||
|
5 | 0.07 | 0.12
|
||||||
|
8 | 0.12 | 0.16
|
||||||
|
12 | 0.21 | 0.25
|
||||||
|
16 | 0.38 | 0.36
|
||||||
|
24 | 0.47 | 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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_a2_d0_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.61/day, rat SD=8.46, residual SD=15.26, days=14
|
||||||
|
|
||||||
|
true stim:day interaction = +1.61 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.22 | 0.53 <- observed
|
||||||
|
5 | 0.53 | 0.68
|
||||||
|
8 | 0.83 | 0.88
|
||||||
|
12 | 0.97 | 1.00
|
||||||
|
16 | 0.98 | 0.98
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.81 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.21 <- observed
|
||||||
|
5 | 0.21 | 0.21
|
||||||
|
8 | 0.42 | 0.42
|
||||||
|
12 | 0.50 | 0.53
|
||||||
|
16 | 0.64 | 0.72
|
||||||
|
24 | 0.78 | 0.75
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_a2_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+6.96/day, rat SD=9.64, residual SD=12.11, days=6
|
||||||
|
|
||||||
|
true stim:day interaction = +6.96 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.43 | 0.85 <- observed
|
||||||
|
5 | 0.87 | 0.96
|
||||||
|
8 | 1.00 | 1.00
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +3.48 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.31 <- observed
|
||||||
|
5 | 0.33 | 0.42
|
||||||
|
8 | 0.57 | 0.67
|
||||||
|
12 | 0.74 | 0.81
|
||||||
|
16 | 0.90 | 0.89
|
||||||
|
24 | 0.97 | 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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_a2_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=-1.80/day, rat SD=10.99, residual SD=10.55, days=5
|
||||||
|
|
||||||
|
true stim:day interaction = -1.80 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.06 | 0.09 <- observed
|
||||||
|
5 | 0.06 | 0.07
|
||||||
|
8 | 0.17 | 0.24
|
||||||
|
12 | 0.30 | 0.30
|
||||||
|
16 | 0.33 | 0.34
|
||||||
|
24 | 0.45 | 0.49
|
||||||
|
|
||||||
|
true stim:day interaction = -0.90 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.03 | 0.07 <- observed
|
||||||
|
5 | 0.07 | 0.07
|
||||||
|
8 | 0.05 | 0.12
|
||||||
|
12 | 0.07 | 0.08
|
||||||
|
16 | 0.11 | 0.11
|
||||||
|
24 | 0.13 | 0.14
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_a2_d6_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=6 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.16/day, rat SD=9.76, residual SD=12.03, days=8
|
||||||
|
|
||||||
|
true stim:day interaction = +1.16 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.08 | 0.15 <- observed
|
||||||
|
5 | 0.19 | 0.23
|
||||||
|
8 | 0.19 | 0.23
|
||||||
|
12 | 0.31 | 0.34
|
||||||
|
16 | 0.42 | 0.47
|
||||||
|
24 | 0.58 | 0.59
|
||||||
|
|
||||||
|
true stim:day interaction = +0.58 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.06 | 0.09 <- observed
|
||||||
|
5 | 0.06 | 0.08
|
||||||
|
8 | 0.09 | 0.11
|
||||||
|
12 | 0.12 | 0.10
|
||||||
|
16 | 0.12 | 0.12
|
||||||
|
24 | 0.21 | 0.23
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_boxa_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.18/day, rat SD=8.67, residual SD=13.71, days=11
|
||||||
|
|
||||||
|
true stim:day interaction = +1.18 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.17 <- observed
|
||||||
|
5 | 0.20 | 0.22
|
||||||
|
8 | 0.32 | 0.38 <- observed
|
||||||
|
12 | 0.58 | 0.62
|
||||||
|
16 | 0.68 | 0.70
|
||||||
|
24 | 0.88 | 0.88
|
||||||
|
|
||||||
|
true stim:day interaction = +0.59 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.05 | 0.07 <- observed
|
||||||
|
5 | 0.04 | 0.07
|
||||||
|
8 | 0.09 | 0.12 <- observed
|
||||||
|
12 | 0.16 | 0.17
|
||||||
|
16 | 0.28 | 0.27
|
||||||
|
24 | 0.28 | 0.29
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_boxa_d0_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.31/day, rat SD=8.76, residual SD=15.18, days=14
|
||||||
|
|
||||||
|
true stim:day interaction = +1.31 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.32 <- observed
|
||||||
|
5 | 0.36 | 0.47
|
||||||
|
8 | 0.63 | 0.68 <- observed
|
||||||
|
12 | 0.94 | 0.94
|
||||||
|
16 | 0.95 | 0.96
|
||||||
|
24 | 0.97 | 0.98
|
||||||
|
|
||||||
|
true stim:day interaction = +0.65 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.17 <- observed
|
||||||
|
5 | 0.17 | 0.18
|
||||||
|
8 | 0.33 | 0.31 <- observed
|
||||||
|
12 | 0.37 | 0.38
|
||||||
|
16 | 0.42 | 0.51
|
||||||
|
24 | 0.62 | 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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_boxa_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=8 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+6.40/day, rat SD=9.10, residual SD=12.48, days=6
|
||||||
|
|
||||||
|
true stim:day interaction = +6.40 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.36 | 0.77 <- observed
|
||||||
|
5 | 0.79 | 0.92
|
||||||
|
8 | 0.98 | 0.98 <- observed
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +3.20 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.24 <- observed
|
||||||
|
5 | 0.26 | 0.39
|
||||||
|
8 | 0.44 | 0.51 <- observed
|
||||||
|
12 | 0.69 | 0.74
|
||||||
|
16 | 0.82 | 0.84
|
||||||
|
24 | 0.95 | 0.94
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_boxa_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=-1.62/day, rat SD=11.24, residual SD=10.14, days=5
|
||||||
|
|
||||||
|
true stim:day interaction = -1.62 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.05 | 0.08 <- observed
|
||||||
|
5 | 0.06 | 0.07
|
||||||
|
8 | 0.17 | 0.21
|
||||||
|
12 | 0.26 | 0.27
|
||||||
|
16 | 0.28 | 0.30
|
||||||
|
24 | 0.41 | 0.45
|
||||||
|
|
||||||
|
true stim:day interaction = -0.81 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.03 | 0.08 <- observed
|
||||||
|
5 | 0.07 | 0.07
|
||||||
|
8 | 0.05 | 0.10
|
||||||
|
12 | 0.07 | 0.07
|
||||||
|
16 | 0.11 | 0.09
|
||||||
|
24 | 0.12 | 0.13
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- naive_boxa_d6_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=7 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+0.87/day, rat SD=10.75, residual SD=11.53, days=8
|
||||||
|
|
||||||
|
true stim:day interaction = +0.87 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.11 <- observed
|
||||||
|
5 | 0.17 | 0.21
|
||||||
|
8 | 0.17 | 0.19
|
||||||
|
12 | 0.23 | 0.23
|
||||||
|
16 | 0.28 | 0.33
|
||||||
|
24 | 0.38 | 0.41
|
||||||
|
|
||||||
|
true stim:day interaction = +0.44 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.06 | 0.10 <- observed
|
||||||
|
5 | 0.05 | 0.08
|
||||||
|
8 | 0.07 | 0.11
|
||||||
|
12 | 0.06 | 0.08
|
||||||
|
16 | 0.07 | 0.07
|
||||||
|
24 | 0.16 | 0.15
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- prev_f_full
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=12, control n=12 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+0.56/day, rat SD=4.32, residual SD=4.49, days=10
|
||||||
|
|
||||||
|
true stim:day interaction = +0.56 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.13 | 0.27
|
||||||
|
5 | 0.30 | 0.41
|
||||||
|
8 | 0.53 | 0.59
|
||||||
|
12 | 0.76 | 0.79 <- observed
|
||||||
|
16 | 0.91 | 0.93
|
||||||
|
24 | 0.98 | 0.98
|
||||||
|
|
||||||
|
true stim:day interaction = +0.28 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.04 | 0.07
|
||||||
|
5 | 0.12 | 0.14
|
||||||
|
8 | 0.17 | 0.23
|
||||||
|
12 | 0.33 | 0.35 <- observed
|
||||||
|
16 | 0.33 | 0.37
|
||||||
|
24 | 0.45 | 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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- right_only_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +1.73 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.22 | 0.39 <- observed
|
||||||
|
5 | 0.52 | 0.65
|
||||||
|
8 | 0.80 | 0.82
|
||||||
|
12 | 0.97 | 0.98
|
||||||
|
16 | 0.98 | 0.98
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.87 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.12 <- observed
|
||||||
|
5 | 0.11 | 0.18
|
||||||
|
8 | 0.31 | 0.32
|
||||||
|
12 | 0.38 | 0.41
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- right_only_d0_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +1.70 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.28 | 0.67 <- observed
|
||||||
|
5 | 0.67 | 0.82
|
||||||
|
8 | 0.92 | 0.94
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.85 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.15 | 0.24 <- observed
|
||||||
|
5 | 0.23 | 0.28
|
||||||
|
8 | 0.48 | 0.53
|
||||||
|
12 | 0.56 | 0.62
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- right_only_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +5.31 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.38 | 0.80 <- observed
|
||||||
|
5 | 0.82 | 0.94
|
||||||
|
8 | 0.99 | 1.00
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +2.65 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.28 <- observed
|
||||||
|
5 | 0.28 | 0.40
|
||||||
|
8 | 0.49 | 0.57
|
||||||
|
12 | 0.72 | 0.77
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- right_only_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +2.42 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.06 | 0.25
|
||||||
|
5 | 0.25 | 0.33
|
||||||
|
8 | 0.38 | 0.44
|
||||||
|
12 | 0.53 | 0.57
|
||||||
|
16 | 0.71 | 0.74
|
||||||
|
24 | 0.94 | 0.94
|
||||||
|
|
||||||
|
true stim:day interaction = +1.21 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.12
|
||||||
|
5 | 0.11 | 0.14
|
||||||
|
8 | 0.09 | 0.14
|
||||||
|
12 | 0.18 | 0.20
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- right_only_d6_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +2.76 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.56 | 0.83
|
||||||
|
5 | 0.93 | 0.98
|
||||||
|
8 | 1.00 | 1.00
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +1.38 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.15 | 0.28
|
||||||
|
5 | 0.36 | 0.53
|
||||||
|
8 | 0.61 | 0.69
|
||||||
|
12 | 0.80 | 0.85
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- unmerge_d0_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success 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
|
||||||
|
|
||||||
|
true stim:day interaction = +1.56 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.15 | 0.31 <- observed
|
||||||
|
5 | 0.44 | 0.56
|
||||||
|
8 | 0.59 | 0.66
|
||||||
|
12 | 0.88 | 0.90
|
||||||
|
16 | 0.93 | 0.95
|
||||||
|
24 | 0.99 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.78 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.11 <- observed
|
||||||
|
5 | 0.09 | 0.13
|
||||||
|
8 | 0.22 | 0.23
|
||||||
|
12 | 0.32 | 0.32
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- unmerge_d0_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.55/day, rat SD=0.00, residual SD=14.68, days=14
|
||||||
|
|
||||||
|
true stim:day interaction = +1.55 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.22 | 0.56 <- observed
|
||||||
|
5 | 0.53 | 0.69
|
||||||
|
8 | 0.82 | 0.88
|
||||||
|
12 | 0.97 | 1.00
|
||||||
|
16 | 0.98 | 0.98
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +0.77 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.22 <- observed
|
||||||
|
5 | 0.21 | 0.23
|
||||||
|
8 | 0.42 | 0.43
|
||||||
|
12 | 0.50 | 0.52
|
||||||
|
16 | 0.64 | 0.71
|
||||||
|
24 | 0.78 | 0.74
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- unmerge_d0_5
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=3 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+5.36/day, rat SD=4.85, residual SD=10.70, days=6
|
||||||
|
|
||||||
|
true stim:day interaction = +5.36 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.35 | 0.73 <- observed
|
||||||
|
5 | 0.78 | 0.90
|
||||||
|
8 | 0.97 | 0.97
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +2.68 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.11 | 0.24 <- observed
|
||||||
|
5 | 0.25 | 0.37
|
||||||
|
8 | 0.44 | 0.46
|
||||||
|
12 | 0.68 | 0.72
|
||||||
|
16 | 0.79 | 0.81
|
||||||
|
24 | 0.95 | 0.94
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- unmerge_d6_10
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=2 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+1.77/day, rat SD=6.11, residual SD=8.63, days=5
|
||||||
|
|
||||||
|
true stim:day interaction = +1.77 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.05 | 0.10 <- observed
|
||||||
|
5 | 0.16 | 0.20
|
||||||
|
8 | 0.17 | 0.28
|
||||||
|
12 | 0.26 | 0.30
|
||||||
|
16 | 0.44 | 0.42
|
||||||
|
24 | 0.57 | 0.58
|
||||||
|
|
||||||
|
true stim:day interaction = +0.88 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.07 | 0.11 <- observed
|
||||||
|
5 | 0.07 | 0.11
|
||||||
|
8 | 0.06 | 0.08
|
||||||
|
12 | 0.12 | 0.13
|
||||||
|
16 | 0.13 | 0.17
|
||||||
|
24 | 0.22 | 0.23
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
@@ -0,0 +1,29 @@
|
|||||||
|
==============================================================================
|
||||||
|
POWER SIMULATION -- unmerge_d6_13
|
||||||
|
==============================================================================
|
||||||
|
model: behavior ~ stim + day + stim:day + (1|rat) (success COUNT; day within-window)
|
||||||
|
observed groups: stim n=3, control n=2 nrep=120, alpha=0.05
|
||||||
|
ground truth: stim:day=+2.27/day, rat SD=7.12, residual SD=7.73, days=8
|
||||||
|
|
||||||
|
true stim:day interaction = +2.27 (100% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.46 | 0.62 <- observed
|
||||||
|
5 | 0.77 | 0.90
|
||||||
|
8 | 0.95 | 0.99
|
||||||
|
12 | 1.00 | 1.00
|
||||||
|
16 | 1.00 | 1.00
|
||||||
|
24 | 1.00 | 1.00
|
||||||
|
|
||||||
|
true stim:day interaction = +1.13 (50% of observed)
|
||||||
|
N/group | per-animal power | LME power
|
||||||
|
------------------------------------------
|
||||||
|
3 | 0.12 | 0.20 <- observed
|
||||||
|
5 | 0.27 | 0.38
|
||||||
|
8 | 0.41 | 0.47
|
||||||
|
12 | 0.54 | 0.62
|
||||||
|
16 | 0.72 | 0.74
|
||||||
|
24 | 0.91 | 0.93
|
||||||
|
|
||||||
|
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.
|
||||||
@@ -0,0 +1,117 @@
|
|||||||
|
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
||||||
|
% interaction, using THIS folder's data as the ground truth.
|
||||||
|
%
|
||||||
|
% 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)
|
||||||
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
||||||
|
% Writes power_result.txt beside this script. 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');
|
||||||
|
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);
|
||||||
|
|
||||||
|
day0 = min(D.day);
|
||||||
|
tbl0 = table(D.success, 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 = [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
|
||||||
|
end
|
||||||
|
|
||||||
|
lme = fitlme(tbl0, FORMULA);
|
||||||
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
||||||
|
b0 = be(strcmp(cn, '(Intercept)')); bStim = be(strcmp(cn, 'stim'));
|
||||||
|
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', ...
|
||||||
|
bInt, sRat, sRes, numel(days))];
|
||||||
|
|
||||||
|
for eMul = EFFMULS
|
||||||
|
bI = bInt * eMul;
|
||||||
|
s = [s sprintf('\n true stim:day interaction = %+.2f (%.0f%% of observed)\n', bI, eMul * 100)]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok<AGROW>
|
||||||
|
s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok<AGROW>
|
||||||
|
for N = NS
|
||||||
|
sigPA = 0; sigL = 0;
|
||||||
|
for r = 1:NREP
|
||||||
|
tb = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes);
|
||||||
|
sigPA = sigPA + (localPAp(tb) < 0.05);
|
||||||
|
try
|
||||||
|
m = fitlme(tb, FORMULA); Cm = m.Coefficients;
|
||||||
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name, 'day:stim')) < 0.05);
|
||||||
|
catch
|
||||||
|
end
|
||||||
|
end
|
||||||
|
star = '';
|
||||||
|
if N == nStim || N == nCtrl; star = ' <- observed'; end
|
||||||
|
s = [s sprintf(' %-8d | %5.2f | %5.2f%s\n', N, sigPA / NREP, sigL / NREP, star)]; %#ok<AGROW>
|
||||||
|
end
|
||||||
|
end
|
||||||
|
|
||||||
|
s = [s sprintf(['\nRead the per-animal column as the honest power. At the observed N this study\n' ...
|
||||||
|
'is typically underpowered; per-animal power reaches ~0.8 only at larger N.\n'])];
|
||||||
|
|
||||||
|
localFinish(s, here);
|
||||||
|
warning(warnState);
|
||||||
|
|
||||||
|
% ---------------------------------------------------------------- helpers
|
||||||
|
function localFinish(s, here)
|
||||||
|
fprintf('%s', s);
|
||||||
|
fid = fopen(fullfile(here, 'power_result.txt'), 'w');
|
||||||
|
fprintf(fid, '%s', s); fclose(fid);
|
||||||
|
end
|
||||||
|
|
||||||
|
function tbl = localSim(N, days, b0, bStim, bDay, bI, sRat, sRes)
|
||||||
|
nd = numel(days); rows = 2 * N * nd;
|
||||||
|
rat = strings(rows, 1); day = zeros(rows, 1); stim = zeros(rows, 1); behavior = zeros(rows, 1);
|
||||||
|
k = 0;
|
||||||
|
for g = 0:1
|
||||||
|
for sIdx = 1:N
|
||||||
|
re = sRat * randn;
|
||||||
|
rid = sprintf('g%d_r%d', g, sIdx);
|
||||||
|
for d = 1:nd
|
||||||
|
k = k + 1;
|
||||||
|
rat(k) = rid; day(k) = days(d); stim(k) = g;
|
||||||
|
behavior(k) = b0 + bStim * g + bDay * days(d) + bI * days(d) * g + re + sRes * randn;
|
||||||
|
end
|
||||||
|
end
|
||||||
|
end
|
||||||
|
tbl = table(categorical(rat), day, stim, behavior, ...
|
||||||
|
'VariableNames', {'rat', 'day', 'stim', 'behavior'});
|
||||||
|
end
|
||||||
|
|
||||||
|
function p = localPAp(tbl)
|
||||||
|
rats = unique(tbl.rat); sl = zeros(numel(rats), 1); gr = zeros(numel(rats), 1);
|
||||||
|
for i = 1:numel(rats)
|
||||||
|
r = tbl.rat == rats(i);
|
||||||
|
c = polyfit(tbl.day(r), tbl.behavior(r), 1); sl(i) = c(1);
|
||||||
|
gr(i) = tbl.stim(find(r, 1));
|
||||||
|
end
|
||||||
|
[~, p] = ttest2(sl(gr == 1), sl(gr == 0), 'Vartype', 'unequal');
|
||||||
|
end
|
||||||
Reference in New Issue
Block a user