diff --git a/analysis/matlab/make_variation_power.m b/analysis/matlab/make_variation_power.m new file mode 100644 index 0000000..7f430ca --- /dev/null +++ b/analysis/matlab/make_variation_power.m @@ -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// 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 diff --git a/analysis/matlab/variation_power.m b/analysis/matlab/variation_power.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variation_power.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_a2_d0_10/power_result.txt b/analysis/matlab/variations/boxa_a2_d0_10/power_result.txt new file mode 100644 index 0000000..9543b93 --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_a2_d0_10/powersim.m b/analysis/matlab/variations/boxa_a2_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_a2_d0_5/power_result.txt b/analysis/matlab/variations/boxa_a2_d0_5/power_result.txt new file mode 100644 index 0000000..2b2b2ae --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_a2_d0_5/powersim.m b/analysis/matlab/variations/boxa_a2_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_a2_d6_10/power_result.txt b/analysis/matlab/variations/boxa_a2_d6_10/power_result.txt new file mode 100644 index 0000000..a531d82 --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_a2_d6_10/powersim.m b/analysis/matlab/variations/boxa_a2_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_a2_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_b2_d0_10/power_result.txt b/analysis/matlab/variations/boxa_b2_d0_10/power_result.txt new file mode 100644 index 0000000..0d19c3b --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_b2_d0_10/powersim.m b/analysis/matlab/variations/boxa_b2_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_b2_d0_5/power_result.txt b/analysis/matlab/variations/boxa_b2_d0_5/power_result.txt new file mode 100644 index 0000000..acb4cba --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_b2_d0_5/powersim.m b/analysis/matlab/variations/boxa_b2_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/boxa_b2_d6_10/power_result.txt b/analysis/matlab/variations/boxa_b2_d6_10/power_result.txt new file mode 100644 index 0000000..8d66ba9 --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/boxa_b2_d6_10/powersim.m b/analysis/matlab/variations/boxa_b2_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/boxa_b2_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/matched_current_d0_3/power_result.txt b/analysis/matlab/variations/matched_current_d0_3/power_result.txt new file mode 100644 index 0000000..1ec16ec --- /dev/null +++ b/analysis/matlab/variations/matched_current_d0_3/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/matched_current_d0_3/powersim.m b/analysis/matlab/variations/matched_current_d0_3/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/matched_current_d0_3/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_a2_d0_10/power_result.txt b/analysis/matlab/variations/naive_a2_d0_10/power_result.txt new file mode 100644 index 0000000..853c454 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_a2_d0_10/powersim.m b/analysis/matlab/variations/naive_a2_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_a2_d0_13/power_result.txt b/analysis/matlab/variations/naive_a2_d0_13/power_result.txt new file mode 100644 index 0000000..4de41e2 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_a2_d0_13/powersim.m b/analysis/matlab/variations/naive_a2_d0_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_a2_d0_5/power_result.txt b/analysis/matlab/variations/naive_a2_d0_5/power_result.txt new file mode 100644 index 0000000..c57b265 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_a2_d0_5/powersim.m b/analysis/matlab/variations/naive_a2_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_a2_d6_10/power_result.txt b/analysis/matlab/variations/naive_a2_d6_10/power_result.txt new file mode 100644 index 0000000..ae734b8 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_a2_d6_10/powersim.m b/analysis/matlab/variations/naive_a2_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_a2_d6_13/power_result.txt b/analysis/matlab/variations/naive_a2_d6_13/power_result.txt new file mode 100644 index 0000000..7b2e583 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d6_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_a2_d6_13/powersim.m b/analysis/matlab/variations/naive_a2_d6_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_a2_d6_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_boxa_d0_10/power_result.txt b/analysis/matlab/variations/naive_boxa_d0_10/power_result.txt new file mode 100644 index 0000000..64b5503 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_boxa_d0_10/powersim.m b/analysis/matlab/variations/naive_boxa_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_boxa_d0_13/power_result.txt b/analysis/matlab/variations/naive_boxa_d0_13/power_result.txt new file mode 100644 index 0000000..37e92dd --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_boxa_d0_13/powersim.m b/analysis/matlab/variations/naive_boxa_d0_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_boxa_d0_5/power_result.txt b/analysis/matlab/variations/naive_boxa_d0_5/power_result.txt new file mode 100644 index 0000000..6223e39 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_boxa_d0_5/powersim.m b/analysis/matlab/variations/naive_boxa_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_boxa_d6_10/power_result.txt b/analysis/matlab/variations/naive_boxa_d6_10/power_result.txt new file mode 100644 index 0000000..0da438b --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_boxa_d6_10/powersim.m b/analysis/matlab/variations/naive_boxa_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/naive_boxa_d6_13/power_result.txt b/analysis/matlab/variations/naive_boxa_d6_13/power_result.txt new file mode 100644 index 0000000..f391e92 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d6_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/naive_boxa_d6_13/powersim.m b/analysis/matlab/variations/naive_boxa_d6_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/naive_boxa_d6_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/prev_f_full/power_result.txt b/analysis/matlab/variations/prev_f_full/power_result.txt new file mode 100644 index 0000000..9368d88 --- /dev/null +++ b/analysis/matlab/variations/prev_f_full/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/prev_f_full/powersim.m b/analysis/matlab/variations/prev_f_full/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/prev_f_full/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/right_only_d0_10/power_result.txt b/analysis/matlab/variations/right_only_d0_10/power_result.txt new file mode 100644 index 0000000..01766dd --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/right_only_d0_10/powersim.m b/analysis/matlab/variations/right_only_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/right_only_d0_13/power_result.txt b/analysis/matlab/variations/right_only_d0_13/power_result.txt new file mode 100644 index 0000000..ac5a438 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/right_only_d0_13/powersim.m b/analysis/matlab/variations/right_only_d0_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/right_only_d0_5/power_result.txt b/analysis/matlab/variations/right_only_d0_5/power_result.txt new file mode 100644 index 0000000..ce0bc1d --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/right_only_d0_5/powersim.m b/analysis/matlab/variations/right_only_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/right_only_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/right_only_d6_10/power_result.txt b/analysis/matlab/variations/right_only_d6_10/power_result.txt new file mode 100644 index 0000000..5374fde --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/right_only_d6_10/powersim.m b/analysis/matlab/variations/right_only_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/right_only_d6_13/power_result.txt b/analysis/matlab/variations/right_only_d6_13/power_result.txt new file mode 100644 index 0000000..2e054dd --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/right_only_d6_13/powersim.m b/analysis/matlab/variations/right_only_d6_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/right_only_d6_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/unmerge_d0_10/power_result.txt b/analysis/matlab/variations/unmerge_d0_10/power_result.txt new file mode 100644 index 0000000..dfc1f55 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/unmerge_d0_10/powersim.m b/analysis/matlab/variations/unmerge_d0_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/unmerge_d0_13/power_result.txt b/analysis/matlab/variations/unmerge_d0_13/power_result.txt new file mode 100644 index 0000000..ee890a5 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/unmerge_d0_13/powersim.m b/analysis/matlab/variations/unmerge_d0_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/unmerge_d0_5/power_result.txt b/analysis/matlab/variations/unmerge_d0_5/power_result.txt new file mode 100644 index 0000000..ba0a995 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_5/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/unmerge_d0_5/powersim.m b/analysis/matlab/variations/unmerge_d0_5/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d0_5/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/unmerge_d6_10/power_result.txt b/analysis/matlab/variations/unmerge_d6_10/power_result.txt new file mode 100644 index 0000000..f3d18a4 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d6_10/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/unmerge_d6_10/powersim.m b/analysis/matlab/variations/unmerge_d6_10/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d6_10/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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 diff --git a/analysis/matlab/variations/unmerge_d6_13/power_result.txt b/analysis/matlab/variations/unmerge_d6_13/power_result.txt new file mode 100644 index 0000000..ed3c599 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d6_13/power_result.txt @@ -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. diff --git a/analysis/matlab/variations/unmerge_d6_13/powersim.m b/analysis/matlab/variations/unmerge_d6_13/powersim.m new file mode 100644 index 0000000..554d032 --- /dev/null +++ b/analysis/matlab/variations/unmerge_d6_13/powersim.m @@ -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 + s = [s sprintf(' %-8s | per-animal power | LME power\n', 'N/group')]; %#ok + s = [s sprintf(' %s\n', repmat('-', 1, 42))]; %#ok + 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 + 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