Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/config/spm_run_bms_map.m
31022b3a13e6master

Last Change
Raw File

spm_run_bms_map.m
View Options

function out = spm_run_bms_map (job)
% Run Bayesian Model Selection Maps
% SPM job execution function
% takes a harvested job data structure and calls SPM functions to perform
% Bayesian Inference for Model Selection of Log. Evidence Maps
% Input:
% job - harvested job data structure (see matlabbatch help)
% Output:
% out - computation results, usually a struct variable.
%
%
% Bayesian Inference on Model Space:
%
% The Random-effects 'RFX' method is described in Stephan et al. [1]
% 'Bayesian Model Selection for Group Studies'.
% Output files (for each model):
% BMS.mat
% Exceedance Probability Maps (*epm.<ext>),
% Posterior Probability Maps (*ppm.<ext>),
% Dirichlet Paramters (alpha) Maps (*alpha.<ext>).
%
% The Fixed-effects 'FFX' method adds together the log-evidences over
% subjects/sessions for each group, then compares the group log-ev's.
% This is also known as the Group Bayes Factor (GBF) approach [2].
% Output files (for each model):
% BMS.mat
% Posterior Probability Maps (*ppm.<ext>).
%
% BMS contains:
% BMS.fname
% BMS.map.ffx(rfx).data
% BMS.map.ffx(rfx).ppm
% BMS.map.ffx(rfx).xppm - only for RFX (this is the expected posterior
% probability map ie. posterior mean)
% BMS.map.ffx(rfx).epm - only for RFX (optional) - this is the
% exceedance probability map
% BMS.map.ffx(rfx).alpha - only for RFX
%
% [1] Rosa et al., 2009, Bayesian Model Selection Maps for Group Studies,
% NeuroImage.
% [2] Stephan et al., 2009, Bayesian Model Selection for Group Studies,
% NeuroImage.
% [3] Penny et al., 2004, Comparing Dynamic Causal Models, NeuroImage.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Maria Joao Rosa
% $Id: spm_run_bms_map.m 5740 2013-11-13 12:00:04Z guillaume $
% Input
% -------------------------------------------------------------------------
direct = job.dir{1};
fname = fullfile(direct,'BMS.mat'); % Output filename (with full path)
mask = length(job.mask{1}); % Mask image
if mask
mask_image = spm_vol(job.mask); % Mask image Vol
end
nsamps = str2num(job.nsamp); % Number of samples (nmodels > 3)
do_maps = job.out_file;
do_ecp = do_maps > 0; % Compute Exceedance Probability
do_alpha = do_maps > 1; % Compute Alpha Parameters
% Nb. of subjects and models
% -------------------------------------------------------------------------
nsubjs = size(job.sess_map,2);
nmodels = size(job.sess_map{1}(1).mod_map,1);
nsess = size(job.sess_map{1},2);
nnames = size(job.mod_name,2);
names = job.mod_name;
% method
% -------------------------------------------------------------------------
if strcmp(job.method_maps,'FFX');
method = 'FFX';
else
method = 'RFX';
end
% Name models
% -------------------------------------------------------------------------
if nnames < nmodels
for nn=1:nmodels-nnames
names = [names, sprintf('m%d',nn)];
end
end
if size(unique(names(1:nmodels)),2) < nmodels,
id = 'Indentical names for different models!'; % Same name!
error(id);
end
if nsubjs == 1
method = 'FFX'; % If only 1 subject do FFX
end
% Sort out log-evidence images dimensions
% -------------------------------------------------------------------------
Vol_models(1,1) = spm_vol(job.sess_map{1}(1).mod_map(1));
first_vol = Vol_models(1,1);
M = first_vol{1}.mat;
DIM = first_vol{1}.dim(1:3)';
xdim = DIM(1);
ydim = DIM(2);
zdim = DIM(3);
[xords,yords] = ndgrid(1:xdim,1:ydim);
xords = xords(:)';
yords = yords(:)';
I = 1:xdim*ydim;
zords_init = ones(1,xdim*ydim);
% Setup images
% -------------------------------------------------------------------------
switch method
case 'FFX', % Fixed Effects
% Check if BMS.mat exists
if exist(fullfile(job.dir{1},'BMS.mat'),'file')
load(fname);
if isfield(BMS,'map') && isfield(BMS.map,'ffx')
str = { 'Warning: existing BMS.mat file has been over-written!'};
msgbox(str)
end
end
% Save BMS data
out.files{1} = fname;
% Create PPM image files for each model
model_ppm(1:nmodels) = struct(...
'fname', '',...
'dim', DIM',...
'dt', [spm_type('float32') spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'n', [1 1], ...
'descrip', '');
% Load Vols for all subjects/models
for i = 1:nmodels
model_ppm(i).fname = fullfile(direct,[sprintf('%s_model_ppm',names{i}) spm_file_ext]);
model_ppm(i).descrip = sprintf('PPM: %s model',names{i});
BMS.map.ffx.ppm{i} = model_ppm(i).fname;
for s = 1:nsubjs,
for se = 1:nsess,
nsessi = size(job.sess_map{s},2);
nmodelsi = size(job.sess_map{s}(se).mod_map,1);
if (nsess == nsessi && nmodels == nmodelsi)
Vol_models(s,i,se) = spm_vol(job.sess_map{s}(se).mod_map(i));
tmp = Vol_models(s,i,se);
else
msgbox('The number of sessions/models should be the same for all subjects!')
return
end
% Stop if log-ev images have different dimensions
if tmp{1}.dim(1)~=xdim || tmp{1}.dim(2)~=ydim || tmp{1}.dim(3)~=zdim
error('Log-evidence images must have the same dimensions!')
end
end
end
end
% Create files
model_ppm = spm_create_vol(model_ppm);
BMS.fname = fname;
% Save data and BMS
BMS.fname = fname;
BMS.map.ffx.data = job.sess_map;
save(out.files{1},'BMS', spm_get_defaults('mat.format'));
case 'RFX', % Random Effects
% Check if BMS.mat exists
if exist(fullfile(job.dir{1},'BMS.mat'),'file')
load(fname);
if isfield(BMS,'map') && isfield(BMS.map,'rfx')
str = { 'Warning: existing BMS.mat file has been over-written!'};
msgbox(str)
end
end
% BMS structure
out.files{1} = fname;
% Create PPM image files for each model
model_exp_r(1:nmodels) = struct(...
'fname', '',...
'dim', DIM',...
'dt', [spm_type('float32') spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'n', [1 1], ...
'descrip', '');
if do_ecp
% Create EPM image files for each model
model_xp(1:nmodels) = struct(...
'fname', '',...
'dim', DIM',...
'dt', [spm_type('float32') spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'n', [1 1], ...
'descrip', '');
end
if do_alpha
% Create alpha image files for each model
model_alpha(1:nmodels) = struct(...
'fname', '',...
'dim', DIM',...
'dt', [spm_type('float32') spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'n', [1 1], ...
'descrip', '');
end
% Load Vols for all subjects/models
for i = 1:nmodels
model_exp_r(i).fname = fullfile(direct,[sprintf('%s_model_xppm',names{i}) spm_file_ext]);
model_exp_r(i).descrip = sprintf('Exp_r: %s model',names{i});
BMS.map.rfx.ppm{i} = model_exp_r(i).fname;
if do_ecp
model_xp(i).fname = fullfile(direct,[sprintf('%s_model_epm',names{i}) spm_file_ext]);
model_xp(i).descrip = sprintf('XP: %s model',names{i});
BMS.map.rfx.epm{i} = model_xp(i).fname;
end
if do_alpha
model_alpha(i).fname = fullfile(direct,[sprintf('%s_model_alpha',names{i}) spm_file_ext]);
model_alpha(i).descrip = sprintf('Alpha: %s model',names{i});
BMS.map.rfx.alpha{i} = model_alpha(i).fname;
end
for s = 1:nsubjs,
for se = 1:nsess,
nsessi = size(job.sess_map{s},2);
nmodelsi = size(job.sess_map{s}(se).mod_map,1);
if (nsess == nsessi && nmodels == nmodelsi)
Vol_models(s,i,se) = spm_vol(job.sess_map{s}(se).mod_map(i));
tmp = Vol_models(s,i,se);
else
msgbox('The number of sessions/models should be the same for all subjects!')
return
end
% Stop if log-ev images have different dimensions
if tmp{1}.dim(1)~=xdim || tmp{1}.dim(2)~=ydim || tmp{1}.dim(3)~=zdim
error('Log-evidence images must have the same dimensions!')
end
end
end
end
% Create files
model_exp_r = spm_create_vol(model_exp_r);
if do_ecp, model_xp = spm_create_vol(model_xp); end
if do_alpha, model_alpha = spm_create_vol(model_alpha); end
% Save data and BMS
BMS.fname = fname;
BMS.map.rfx.data = job.sess_map;
save(out.files{1},'BMS', spm_get_defaults('mat.format'));
end
% Progress bar
% -------------------------------------------------------------------------
spm_progress_bar('Init',zdim,'BMS Maps (Inference)','Slices complete');
% Loop through image slices
% -------------------------------------------------------------------------
for z = 1:zdim,
spm_progress_bar('Set',z); % Update progress bar
j = repmat(NaN,xdim,ydim); % Init. image values
fprintf('%s%30s',repmat(sprintf('\b'),1,30),'Computing maps...')
str = sprintf('Slice %d out of %d',z,zdim); % Display slice nb.
fprintf('\r%-40s: %30s',str,' ')
zords = z*zords_init; % Slice z
xyz = [xords(I); yords(I); zords(I)]; % Slice coordinates
nVox = size(xyz,2); % Nb. of voxels per slice
if mask
% Voxels inside mask
mask_xyz = mask_image{1}.mat\M*[xyz(:,1:nVox);ones(1,nVox)];
gamma = spm_get_data(mask_image{1},mask_xyz);
b = find(gamma>0.5); % Voxels in the mask
else
b = 1:nVox; % All voxels
end
z_models = NaN(nsubjs,nmodels,nVox); % Data
z_models(1,1,:) = spm_get_data(first_vol{1},xyz); % Data: all subs/mods
non_nan = find(~isnan(z_models(1,1,:))); % Voxels ~NaN
% Find voxels ~NaN and sum sessions
% ---------------------------------------------------------------------
for s = 1:nsubjs,
for k = 1:nmodels,
sum_tmp_data = [];
for ns = 1:nsess,
tmp_data = Vol_models(s,k,ns);
sum_tmp_data = [sum_tmp_data; spm_get_data(tmp_data{1},xyz)];
end
z_models(s,k,:) = sum(sum_tmp_data,1);
non_nani = find(~isnan(z_models(s,k,:)));
non_nan = intersect(non_nan,non_nani);
end
end
% Voxels to be analysed
non_nan = intersect(non_nan,b);
Nvoxels = length(non_nan);
% Method
% ---------------------------------------------------------------------
switch method
% Fixed Effects
% ---------------------------------------------------------------
case 'FFX',
if Nvoxels > 0 % Slice with ~NaN voxels
zz = sum(z_models,1); % Sum all subjects/sessions
mz = mean(zz,2); % Get mean of all models
zzmean = zeros(1,nmodels,length(mz));
for jj = 1:nmodels
zzmean(1,jj,:) = mz;
end
if nmodels==1
% Process log Bayes factor image
zz = exp(zz); % Exponentiate log-bf values
% Calculate posterior probabiliy
pz = zeros(1,1,length(zz));
pz = zz./(1+zz);
j(non_nan) = pz(1,1,non_nan);
model_ppm(1) = spm_write_plane(model_ppm(k),j,z);
else
zz = zz-zzmean; % Subtract mean
zz = exp(zz); % Exponentiate log-ev values
tzz = sum(zz,2); % Sum exp(log-ev.)
% Calculate posterior probabiliy
pz = zeros(1,nmodels,length(zz));
for k = 1:nmodels,
pz(1,k,:) = zz(1,k,:)./tzz;
j(non_nan) = pz(1,k,non_nan);
model_ppm(k) = spm_write_plane(model_ppm(k),j,z);
end
end
else
% Nvoxels = 0
for k = 1:nmodels,
% Write NaN for slice z
model_ppm(k) = spm_write_plane(model_ppm(k),j,z);
end
end
% Fixed Effects
% ---------------------------------------------------------------
case 'RFX',
if Nvoxels > 0
% Initialise results
exp_r_total = zeros(Nvoxels,nmodels);
if do_ecp, xp_total = zeros(Nvoxels,nmodels); end
if do_alpha, alpha_total = zeros(Nvoxels,nmodels); end
% Do BMS in all voxels of slice z
for n = 1:Nvoxels,
lme = z_models(:,:,non_nan(n));
if nmodels==1
% Provide evidence for dummy null model
lme=0.5*[lme, -lme];
[alpha,exp_r,xp] = spm_BMS(lme,nsamps,0,0,do_ecp);
exp_r_total(n,:) = exp_r(1); % Cond. Expecta.
if do_ecp, xp_total(n,:) = xp(1); end % Exceeda. Prob.
if do_alpha, alpha_total(n,:) = alpha(1); end % Dirichlet par.
else
% Group BMS
[alpha,exp_r,xp] = spm_BMS(lme,nsamps,0,0,do_ecp);
exp_r_total(n,:) = exp_r; % Cond. Expecta.
if do_ecp, xp_total(n,:) = xp; end % Exceeda. Prob.
if do_alpha, alpha_total(n,:) = alpha; end % Dirichlet par.
end
end
% Write images
for i = 1:nmodels,
j(non_nan) = exp_r_total(:,i);
model_exp_r(i) = spm_write_plane(model_exp_r(i),j,z);
if do_ecp
j(non_nan) = xp_total(:,i);
model_xp(i) = spm_write_plane(model_xp(i),j,z);
end
if do_alpha
j(non_nan) = alpha_total(:,i);
model_alpha(i) = spm_write_plane(model_alpha(i),j,z);
end
end
else
% Write images when Nvoxels = 0
for i = 1:nmodels,
model_exp_r(i) = spm_write_plane(model_exp_r(i),j,z);
if do_ecp, model_xp(i) = spm_write_plane(model_xp(i),j,z); end
if do_alpha, model_alpha(i) = spm_write_plane(model_alpha(i),j,z); end
end
end
end
end % Loop over slices
% Clear progress bar
% -------------------------------------------------------------------------
spm_progress_bar('Clear');
disp('Done.');
c4science · Help