Experiments DB Dev
7b5c27df8c
Adds tdcs_paper_lme replicating the paper's exact formula/variable names (paper_* switch cases): mergeA2 reproduces the interaction F(1)=7.09 vs 7.12, p=0.009 vs 0.008, with a coverage-confound warning. Rewrites tdcs_power_sim to fit that same formula (interaction term canonicalized to day:stim). Adds a replication test. Suite 37/37. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
102 lines
4.2 KiB
Matlab
102 lines
4.2 KiB
Matlab
function PW = tdcs_power_sim(mergeKey, window, Ns, effMuls, nrep, cfg)
|
|
%TDCS_POWER_SIM Monte-Carlo power for the paper's stim x day interaction.
|
|
% Uses the paper's exact model -- fitlme(tbl, 'behavior ~ stim + day +
|
|
% stim:day + (1|rat)') -- throughout. The fitted model for MERGEKEY over the
|
|
% WINDOW = [lo hi] day range (Box-B2 = stim = 1 vs Box-A2 = stim = 0) is the
|
|
% ground truth (its fixed effects, rat random-intercept SD, residual SD);
|
|
% NREP datasets are simulated at each rats-per-group in NS, for each
|
|
% true-effect multiplier in EFFMULS (e.g. [1 0.5] = observed and half the
|
|
% observed stim:day slope). Each dataset is scored at alpha = 0.05 two ways:
|
|
% - per-animal (cluster-honest): per-rat behavior~day slope, Welch t (stim)
|
|
% - LME: the fitlme stim:day interaction p (observation-level DF)
|
|
% PW is a table (effMul, N, powerPerAnimal, powerLME); also printed. A fixed
|
|
% RNG seed makes the estimate reproducible.
|
|
%
|
|
% Defaults: WINDOW=[0 5] (early phase), 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(window); window = [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);
|
|
FORMULA = 'behavior ~ stim + day + stim:day + (1|rat)';
|
|
|
|
% Ground truth = fitted paper model over the window.
|
|
Sfull = tdcs_scenario_data([mergeKey '_full']);
|
|
A = Sfull(ismember(Sfull.group, {cfg.anchorLow, cfg.anchorHigh}), :);
|
|
A = A(A.day >= window(1) & A.day <= window(2), :);
|
|
tbl0 = table();
|
|
tbl0.behavior = A.success;
|
|
tbl0.day = A.day - window(1); % 0-based within the window
|
|
tbl0.stim = double(A.group == cfg.anchorHigh);
|
|
tbl0.rat = A.subject;
|
|
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 = (window(1):window(2))' - window(1);
|
|
|
|
fprintf('Power sim (paper formula): truth=%s window %d-%d | stim:day=%.2f ratSD=%.2f resSD=%.2f | nrep=%d\n', ...
|
|
mergeKey, window(1), window(2), bInt, sRat, sRes, nrep);
|
|
|
|
rows = {};
|
|
for eMul = effMuls
|
|
bI = bInt * eMul;
|
|
fprintf('\n true stim:day 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, bStim, bDay, bI, sRat, sRes);
|
|
sigPA = sigPA + (localPerAnimalP(tbl) < 0.05);
|
|
try
|
|
m = fitlme(tbl, FORMULA);
|
|
Cm = m.Coefficients;
|
|
sigL = sigL + (Cm.pValue(strcmp(Cm.Name,'day:stim')) < 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, 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 = localPerAnimalP(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
|