Skip to content

Instantly share code, notes, and snippets.

@jooh jooh/test_rfx_ffx.m
Created Nov 17, 2017

Embed
What would you like to do?
Matlab simulation comparing permutation tests to parametric RFX and FFX inference
% Explore correspondence between different types of permutation tests and
% parametric stats inference (FFX/RFX). We test two classes of permutation test:
% a within-subject permutation approach (Stelzer et al., 2013, NI), and the
% standard sign-flip permutation test (Nichols/Holmes, 2001, HBM).
%
% OBSERVATIONS:
% * RFX parametric and sign flip permutation p values are highly similar, as
% expected
% * FFX parametric and within-subject permutation p values are highly similar
% * using T as the test stat for within-subject permutation test seems to make
% the test behave very similarly to the RFX sign flip test. Intriguing, but this
% may not be a general result - this t stat is FFX and in the current similation
% within and between subject variance is the same, which is a much simpler case
% than what we expect in the real world (>within variance if modelling e.g.
% volumes, <within if modelling e.g. parameter estimates). Probably needs
% further testing for general use. In any case, I've never seen anyone do this -
% when this approach is used, it's always using mean as the test stat.
%
% 2017-04-28 J Carlin
%
% test_rfx_ffx(recompute)
function test_rfx_ffx(recompute)
if ~exist('recompute','var') || isempty(recompute)
recompute = false;
end
% data with a slight signal (mean effect > 0)
par.nobs = 100;
par.nsub = 10;
% predicted effect alternates between positive and negative
par.regressor = repmat([-1 1],[1, par.nobs/2]);
par.nsim = 100;
par.nsig = 5;
par.siglevels = logspace(-3,0,par.nsig);
par.nperm = 1000;
[fundir,funname,ext] = fileparts(mfilename('fullpath'));
simres = fullfile(fundir,[funname '_simres.mat']);
if recompute || ~exist(simres,'file')
logstr('running %d iterations\n',par.nsim*par.nsig*par.nperm);
parfor sim = 1:par.nsim
data = randn(par.nsub,par.nobs);
simres_p_rfx_para = NaN([1,par.nsig]);
simres_p_rfx_perm = NaN([1,par.nsig]);
simres_p_ffx_para = NaN([1,par.nsig]);
simres_p_ffx_m_perm = NaN([1,par.nsig]);
simres_p_ffx_t_perm = NaN([1,par.nsig]);
for sig = 1:par.nsig
sigdata = bsxfun(@plus,data,par.regressor*par.siglevels(sig));
[simres_p_rfx_para(sig),...
simres_p_rfx_perm(sig),...
simres_p_ffx_para(sig),...
simres_p_ffx_m_perm(sig),...
simres_p_ffx_t_perm(sig)] = runsim(sigdata,par);
end
res_p_rfx_para{sim} = simres_p_rfx_para;
res_p_rfx_perm{sim} = simres_p_rfx_perm;
res_p_ffx_para{sim} = simres_p_ffx_para;
res_p_ffx_m_perm{sim} = simres_p_ffx_m_perm;
res_p_ffx_t_perm{sim} = simres_p_ffx_t_perm;
end
logstr('saving to %s\n',simres);
save(simres,'res_p_rfx_para','res_p_rfx_perm','res_p_ffx_para',...
'res_p_ffx_m_perm','res_p_ffx_t_perm','par');
else
oldpar = par;
logstr('loading existing result from %s\n',simres);
load(simres);
assert(isequal(oldpar,par),'par changed relative since saved result was computed. Rerun?');
end
% now just summarise
means_p_rfx_para = matfun(@mean,res_p_rfx_para{:});
means_p_rfx_perm = matfun(@mean,res_p_rfx_perm{:});
means_p_ffx_para = matfun(@mean,res_p_ffx_para{:});
means_p_ffx_m_perm = matfun(@mean,res_p_ffx_m_perm{:});
means_p_ffx_t_perm = matfun(@mean,res_p_ffx_t_perm{:});
stdev = @(x,dim)std(x,[],dim);
std_p_rfx_para = matfun(stdev,res_p_rfx_para{:});
std_p_rfx_perm = matfun(stdev,res_p_rfx_perm{:});
std_p_ffx_para = matfun(stdev,res_p_ffx_para{:});
std_p_ffx_m_perm = matfun(stdev,res_p_ffx_m_perm{:});
std_p_ffx_t_perm = matfun(stdev,res_p_ffx_t_perm{:});
% and plot as a function of siglevel - one panel for RFX, one for FFX
F = figure(100);
clf(F);
% scatter with 2d error bars would be cool
subplot(1,2,1);
errorbar(repmat(par.siglevels',1,2),[means_p_rfx_para' means_p_rfx_perm'],...
[std_p_rfx_para',std_p_rfx_perm']);
title('RFX');
subplot(1,2,2);
errorbar(repmat(par.siglevels',1,2),[means_p_ffx_para' means_p_ffx_m_perm'],...
[std_p_ffx_para',std_p_ffx_m_perm']);
title('FFX');
F = figure(200);
clf(F);
subplot(1,2,1);
plot(horzcat(res_p_rfx_para{:}),horzcat(res_p_rfx_perm{:}),'o');
axis([0,1,0,1]);
set(gca,'dataaspectratio',[1,1,1],'xtick',[0,1],'ytick',[0,1]);
xlabel('parametric p');
ylabel('permutation p');
line([0,1],[0,1],'color',[0,0,0]);
title('RFX');
subplot(1,2,2);
plot(horzcat(res_p_ffx_para{:}),horzcat(res_p_ffx_m_perm{:}),'o');
axis([0,1,0,1]);
set(gca,'dataaspectratio',[1,1,1],'xtick',[0,1],'ytick',[0,1]);
xlabel('parametric p');
ylabel('permutation p');
line([0,1],[0,1],'color',[0,0,0]);
title('FFX');
printstandard(fullfile(fundir,funname));
keyboard;
function [rfx_p_para,rfx_p_perm,ffx_p_para,ffx_m_p,ffx_t_p] = runsim(data,par)
% generate single subject effects (parameter estimates from GLM)
mdata = arrayfun(@(x)par.regressor' \ data(x,:)',1:par.nsub);
[~,rfx_p_para,~,~] = ttest(mdata,0,'tail','right');
% parametric FFX p value
ffxdesign = repmat(par.regressor,[par.nsub,1]);
ffxmodel = GLM(ffxdesign(:),data(:));
ffx_p_para = pmap(ffxmodel,1,'right');
% permutation test by permuting observations within subjects
% nb these functions return true perm as first entry
ffx_inds = permuteindices(par.nobs,par.nperm);
rfx_inds = permflipindices(par.nsub,par.nperm);
ffx_nulldist_m = NaN([par.nperm,1]);
ffx_nulldist_t = NaN([par.nperm,1]);
rfx_nulldist_t = NaN([par.nperm,1]);
for p = 1:par.nperm
% FFX
% permute the design matrix (equivalent to permuting data) and re-calculate
% for each subject
permreg = par.regressor(ffx_inds(p,:));
mdata_null_ffx = arrayfun(@(x)permreg' \ data(x,:)',1:par.nsub);
ffx_nulldist_m(p) = mean(mdata_null_ffx);
[~,~,~,stats_null_ffx] = ttest(mdata_null_ffx,0,'tail','right');
ffx_nulldist_t(p) = stats_null_ffx.tstat;
% RFX
% just flip the sign and re-calculate test stat (nichols/holmes)
mdata_null_rfx = mdata;
mdata_null_rfx(rfx_inds(p,:)) = mdata_null_rfx(rfx_inds(p,:)) * -1;
[~,~,~,stats_null_rfx] = ttest(mdata_null_rfx,0,'tail','right');
rfx_nulldist_t(p) = stats_null_rfx.tstat;
end
% rfx p value is just percentile of trust stat in null dist
rfx_p_perm = sum(rfx_nulldist_t >= rfx_nulldist_t(1)) / par.nperm;
% nb most flavours use
ffx_m_p = sum(ffx_nulldist_m >= ffx_nulldist_m(1)) / par.nperm;
ffx_t_p = sum(ffx_nulldist_t >= ffx_nulldist_t(1)) / par.nperm;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.