Experiments DB Dev
e730aa0020
Adds tdcs_phase_lme (same LME refit within learning phases 0-5/6-10/6-13, with per-phase slopes + interaction 95% CI) as phase_* switch cases, and tdcs_power_sim (Monte-Carlo power for the early-phase interaction, scored by the cluster-honest per-animal test and the LME test) + run_power. Findings: the Box-B2 faster-acquisition signal is in the early phase; at n=3-5/group honest power is 0.3-0.7 even at the observed effect (need ~8/group if effect is as observed, ~20/group if half). Adds tPhasePower tests. Suite 35/35. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
98 lines
4.1 KiB
Matlab
98 lines
4.1 KiB
Matlab
function PW = tdcs_power_sim(mergeKey, phase, Ns, effMuls, nrep, cfg)
|
|
%TDCS_POWER_SIM Monte-Carlo power for the phased days x tDCS interaction.
|
|
% PW = TDCS_POWER_SIM(MERGEKEY, PHASE, NS, EFFMULS, NREP, CFG) estimates the
|
|
% power to detect the day x tDCS interaction of the phased LME
|
|
% (success ~ dayp*tDCS + (1|subject)). The fitted model for MERGEKEY over the
|
|
% PHASE = [lo hi] window is used as ground truth (its fixed effects, subject
|
|
% random-intercept SD, and residual SD); NREP datasets are simulated at each
|
|
% subjects-per-group in NS, for each true-effect multiplier in EFFMULS (e.g.
|
|
% [1 0.5] = observed and half the observed interaction). Each dataset is
|
|
% scored two ways at alpha = 0.05:
|
|
% - per-animal (cluster-honest): per-subject slope, Welch t (Box-B2 vs A2)
|
|
% - LME: the fitlme dayp:tDCS interaction p (observation-level DF)
|
|
% PW is a table (effMul, N, powerPerAnimal, powerLME); it is also printed.
|
|
% A fixed RNG seed makes the estimate reproducible.
|
|
%
|
|
% Defaults: PHASE=[0 5], NS=[3 5 8 12 16 24 30], EFFMULS=[1 0.5], NREP=200.
|
|
|
|
if nargin < 1 || isempty(mergeKey); mergeKey = 'mergeA2'; end
|
|
if nargin < 2 || isempty(phase); phase = [0 5]; end
|
|
if nargin < 3 || isempty(Ns); Ns = [3 5 8 12 16 24 30]; end
|
|
if nargin < 4 || isempty(effMuls); effMuls = [1 0.5]; end
|
|
if nargin < 5 || isempty(nrep); nrep = 200; end
|
|
if nargin < 6 || isempty(cfg); cfg = tdcs_config(); end
|
|
|
|
warnState = warning('off', 'all'); cleanupW = onCleanup(@() warning(warnState)); %#ok<NASGU>
|
|
rng(1); % reproducible
|
|
|
|
% Ground truth = fitted phased LME.
|
|
Sfull = tdcs_scenario_data([mergeKey '_full']);
|
|
T = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
|
|
Tp = T(T.day >= phase(1) & T.day <= phase(2), :);
|
|
Tp.dayp = Tp.day - phase(1);
|
|
Tp.tDCS = double(Tp.group == cfg.anchorHigh);
|
|
lme = fitlme(Tp, 'success ~ dayp*tDCS + (1|subject)');
|
|
cn = lme.CoefficientNames; be = lme.fixedEffects;
|
|
b0 = be(strcmp(cn,'(Intercept)')); bDay = be(strcmp(cn,'dayp'));
|
|
bT = be(strcmp(cn,'tDCS')); bInt = be(strcmp(cn,'dayp:tDCS'));
|
|
psi = covarianceParameters(lme); sSub = sqrt(psi{1}); sRes = sqrt(lme.MSE);
|
|
days = (phase(1):phase(2))' - phase(1);
|
|
|
|
fprintf('Power sim: truth=%s phase %d-%d | interaction=%.2f subjSD=%.2f resSD=%.2f | nrep=%d\n', ...
|
|
mergeKey, phase(1), phase(2), bInt, sSub, sRes, nrep);
|
|
|
|
rows = {};
|
|
for eMul = effMuls
|
|
bI = bInt * eMul;
|
|
fprintf('\n true interaction = %+.2f (%.0f%% of observed)\n', bI, eMul*100);
|
|
fprintf(' %-8s | per-animal power | LME power\n', 'N/group');
|
|
for N = Ns
|
|
sigPA = 0; sigL = 0;
|
|
for r = 1:nrep
|
|
tbl = localSim(N, days, b0, bDay, bT, bI, sSub, sRes);
|
|
sigPA = sigPA + (localPerAnimalP(tbl) < 0.05);
|
|
try
|
|
m = fitlme(tbl, 'success ~ dayp*tDCS + (1|subject)');
|
|
Cm = m.Coefficients;
|
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'dayp:tDCS')) < 0.05);
|
|
catch
|
|
end
|
|
end
|
|
pPA = sigPA / nrep; pL = sigL / nrep;
|
|
fprintf(' %-8d | %5.2f | %5.2f\n', N, pPA, pL);
|
|
rows(end+1, :) = {eMul, N, pPA, pL}; %#ok<AGROW>
|
|
end
|
|
end
|
|
PW = cell2table(rows, 'VariableNames', {'effMul','N','powerPerAnimal','powerLME'});
|
|
|
|
end
|
|
|
|
function tbl = localSim(N, days, b0, bDay, bT, bI, sSub, sRes)
|
|
nd = numel(days); rows = 2*N*nd;
|
|
subj = strings(rows,1); dayp = zeros(rows,1); tDCS = zeros(rows,1); success = zeros(rows,1);
|
|
k = 0;
|
|
for g = 0:1
|
|
for sIdx = 1:N
|
|
re = sSub * randn;
|
|
sid = sprintf('g%d_s%d', g, sIdx);
|
|
for d = 1:nd
|
|
k = k+1;
|
|
subj(k) = sid; dayp(k) = days(d); tDCS(k) = g;
|
|
success(k) = b0 + bDay*days(d) + bT*g + bI*days(d)*g + re + sRes*randn;
|
|
end
|
|
end
|
|
end
|
|
tbl = table(categorical(subj), dayp, tDCS, success, ...
|
|
'VariableNames', {'subject','dayp','tDCS','success'});
|
|
end
|
|
|
|
function p = localPerAnimalP(tbl)
|
|
subs = unique(tbl.subject); sl = zeros(numel(subs),1); gr = zeros(numel(subs),1);
|
|
for i = 1:numel(subs)
|
|
r = tbl.subject == subs(i);
|
|
c = polyfit(tbl.dayp(r), tbl.success(r), 1); sl(i) = c(1);
|
|
gr(i) = tbl.tDCS(find(r,1));
|
|
end
|
|
[~, p] = ttest2(sl(gr==1), sl(gr==0), 'Vartype', 'unequal');
|
|
end
|