Experiments DB Dev
119 lines
4.9 KiB
Matlab
119 lines
4.9 KiB
Matlab
% Variation power simulation -- Monte-Carlo power for the paper's stim x day
|
|
% interaction, using THIS folder's data as the ground truth, for BOTH metrics:
|
|
% metric = count : behavior = # successes -> power_result.txt
|
|
% metric = rate : behavior = success / attempts -> power_result_rate.txt
|
|
%
|
|
% Ground truth: fitlme(behavior ~ stim + day + stim:day + (1|rat)) on data.csv
|
|
% (day within-window). Its fixed effects, per-rat intercept SD, and residual SD
|
|
% generate NREP synthetic datasets at each rats-per-group N and each true-effect
|
|
% multiplier (1 = observed, 0.5 = half). Each is scored at alpha=0.05 by:
|
|
% per-animal : Welch t on per-rat behavior~day slopes (cluster-honest power)
|
|
% LME : the fitlme stim:day p (observation-level DF -- optimistic)
|
|
% Writes power_result[_rate].txt. Run: matlab -batch "powersim"
|
|
% (Copy of analysis/matlab/variation_power.m; see make_variation_power.m.)
|
|
|
|
here = fileparts(mfilename('fullpath'));
|
|
if isempty(here); here = pwd; end
|
|
vname = regexprep(here, '.*[/\\]', '');
|
|
D = readtable(fullfile(here, 'data.csv'), 'TextType', 'string');
|
|
|
|
localPower(D, 'count', here, vname);
|
|
localPower(D, 'rate', here, vname);
|
|
|
|
% ---------------------------------------------------------------- per metric
|
|
function localPower(D, metric, here, vname)
|
|
NS = [3 5 8 12 16 24]; EFFMULS = [1 0.5]; NREP = 120;
|
|
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
|
if strcmp(metric, 'rate')
|
|
D = D(D.total > 0, :); beh = D.success ./ D.total; mlabel = 'success RATE'; suffix = '_rate';
|
|
else
|
|
beh = D.success; mlabel = '# successes (count)'; suffix = '';
|
|
end
|
|
warnState = warning('off', 'all'); rng(1);
|
|
day0 = min(D.day);
|
|
tbl0 = table(beh, D.day - day0, double(D.stim), categorical(D.subject), ...
|
|
'VariableNames', {'behavior', 'day', 'stim', 'rat'});
|
|
nStim = numel(unique(D.subject(D.stim == 1)));
|
|
nCtrl = numel(unique(D.subject(D.stim == 0)));
|
|
|
|
bar = repmat('=', 1, 78);
|
|
s = sprintf('%s\nPOWER SIMULATION -- %s [metric: %s]\n%s\n', bar, vname, mlabel, bar);
|
|
s = [s sprintf('model: %s (behavior = %s; day within-window)\n', FORMULA, mlabel)];
|
|
s = [s sprintf('observed groups: stim n=%d, control n=%d nrep=%d, alpha=0.05\n', nStim, nCtrl, NREP)];
|
|
|
|
if nStim < 2 || nCtrl < 2 || numel(unique(tbl0.day)) < 2
|
|
s = [s sprintf('\nInsufficient data for a power simulation.\n')];
|
|
localFinish(s, here, suffix); 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=%+.4g/day, rat SD=%.4g, residual SD=%.4g, days=%d\n', ...
|
|
bInt, sRat, sRes, numel(days))];
|
|
|
|
for eMul = EFFMULS
|
|
bI = bInt * eMul;
|
|
s = [s sprintf('\n true stim:day interaction = %+.4g (%.0f%% of observed)\n', bI, eMul * 100)];
|
|
s = [s sprintf(' %-8s | per-animal power | LME power\n %s\n', 'N/group', repmat('-', 1, 42))];
|
|
for N = NS
|
|
sigPA = 0; sigL = 0;
|
|
for r = 1:NREP
|
|
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)];
|
|
end
|
|
end
|
|
s = [s sprintf('\nRead the per-animal column as the honest power; LME is optimistic (obs-level DF).\n')];
|
|
localFinish(s, here, suffix);
|
|
warning(warnState);
|
|
end
|
|
|
|
% ---------------------------------------------------------------- helpers
|
|
function localFinish(s, here, suffix)
|
|
fprintf('%s', s);
|
|
fid = fopen(fullfile(here, ['power_result' suffix '.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
|