Matlab simulation comparing permutation tests to parametric RFX and FFX inference

test_rfx_ffx.m

% 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;