Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/packages/spm12/config/spm_run_dcm_bms.m
37b12c85c481master

Last Change
Raw File

spm_run_dcm_bms.m
View Options

function out = spm_run_dcm_bms(job)
% Compare DCMs on the basis of their log-evidences
% Four methods are available to identify the best among alternative models:
%
% (1) single subject BMS using Bayes factors
% (see Penny et al, NeuroImage, 2004)
% (2) fixed effects group BMS using group Bayes factors
% (see Stephan et al, NeuroImage, 2007)
% (3) random effects group BMS using exceedance probabilities
% (see Stephan et al, NeuroImage, 2009)
% (4) comparing model families
% (see Penny et al, PLOS-CB, 2010)
%
% Note: All functions use the negative free energy (F) as an approximation
% to the log model evidence.
%__________________________________________________________________________
% Copyright (C) 2009-2015 Wellcome Trust Centre for Neuroimaging
% CC Chen & Maria Joao Rosa
% $Id: spm_run_dcm_bms.m 7101 2017-06-08 12:36:13Z peter $
SVNid = '$Rev: 7101 $';
%-Say hello
%--------------------------------------------------------------------------
spm('FnBanner',mfilename,SVNid);
%-Input parameters
%--------------------------------------------------------------------------
fname = 'BMS.mat'; % Output filename
fname = fullfile(job.dir{1},fname); % Output filename (inc. directory)
ld_f = ~isempty(job.load_f{1}); % Log-evidence file
ld_msp = ~isempty(job.model_sp{1}); % Model space file
bma_do = isfield(job.bma,'bma_yes'); % Compute BMA
data_se = ~isempty(job.sess_dcm); % DCM files
out.bmsmat = {''};
% check whether a BMS.mat file exists in the selected directory
if exist(fname,'file')
str = {'Current directory contains a BMS file:',...
'Continuing will overwrite existing file !!'};
if spm_input(str,1,'bd','stop|continue',[1,0],1,mfilename);
fprintf('%-40s: %30s\n\n',...
'BMS stoped...existing BMS file in the selected directory',...
spm('time'));
Fgraph = spm_figure('GetWin','Graphics');
spm_clf(Fgraph);
return
end
end
%-Method
% -------------------------------------------------------------------------
if strcmp(job.method,'FFX');
method = 'FFX'; % Fixed-effects
else
method = 'RFX'; % Random-effects
end
%-Check DCM.mat files and BMA
%--------------------------------------------------------------------------
do_bma_famwin = 0;
do_bma_all = 0;
if bma_do
if ld_msp
fprintf('%-40s: ','Loading model space') %-#
load(job.model_sp{1});
load(subj(1).sess(1).model(1).fname)
fprintf('%30s\n','...done') %-#
else
if data_se
load(job.sess_dcm{1}(1).dcmmat{1})
else
spm('alert*','Please specify DCM.mat files or model_space.mat to do BMA.','Error');
return
end
end
nsamp = 1e4;
odds_ratio = 1/20;
if isfield(job.bma.bma_yes,'bma_famwin')
do_bma_famwin = 1;
else
if isfield(job.bma.bma_yes,'bma_all')
do_bma_all = 1;
else
bma_fam = double(job.bma.bma_yes.bma_part);
end
end
else
if ld_msp
fprintf('%-40s: ','Loading model space') %-#
load(job.model_sp{1},'subj');
fprintf('%30s\n','...done') %-#
if ~exist('subj','var')
spm('alert*',{'Incorrect model space file.','File must contain ''subj'' structure.'},'Error');
return
end
end
end
%-Initialise
%--------------------------------------------------------------------------
F = []; % Free energy
N = {}; % Models
%-Prepare the data
%--------------------------------------------------------------------------
if ld_f
data = job.load_f{1};
load(data,'F');
if ~exist('F','var')
spm('alert*',{'Incorrect log-evidence matrix file.','File must contain ''F'' matrix.'},'Error');
return
end
nm = size(F,2); % No of Models
ns = size(F,1); % No of Subjects
N = 1:nm;
subj = []; % subj structure
f_fname = data; % LogE file name
fname_msp = []; % Model space
else
f_fname = []; % LogE file name
% Verify dimensions of data and model space
if ld_msp
ns = length(subj); % No of Subjects
nsess = length(subj(1).sess); % No of sessions
nm = length(subj(1).sess(1).model); % No of Models
fname_msp = [job.model_sp{1}];
else
ns = size(job.sess_dcm,2); % No of Subjects
nsess = size(job.sess_dcm{1},2); % No of sessions
nm = size(job.sess_dcm{1}(1).dcmmat,1); % No of Models
fname_msp = fullfile(job.dir{1},'model_space.mat');
end
F = zeros(ns,nm);
N = cell(nm);
% Check if No of models > 2
if nm < 2
spm('alert*','Please select more than one file.','Error');
return
end
for k=1:ns
% Verify dimensions of data and model space for current session
if ld_msp
nsess_now = length(subj(k).sess);
nmodels = length(subj(k).sess(1).model);
else
disp(sprintf('Loading DCMs for subject %d', k));
nsess_now = size(job.sess_dcm{k},2);
nmodels = size(job.sess_dcm{k}(1).dcmmat,1);
end
if (nsess_now == nsess && nmodels== nm) % Check no of sess/mods
ID = zeros(nm,nsess);
for j=1:nm
F_sess = [];
for h = 1:nsess_now
if ~ld_msp
clear DCM
% Load DCM (model)
tmp = job.sess_dcm{k}(h).dcmmat{j};
if ~job.verify_id
st = warning('off','MATLAB:load:variableNotFound');
DCM.DCM = load(tmp, 'F', 'Ep', 'Cp');
warning(st);
if ~isfield(DCM.DCM, 'F')
DCM = load(tmp);
end
else
DCM = load(tmp);
end
% Free energy for sessions
F_sess = [F_sess,DCM.DCM.F];
% Create model space
subj(k).sess(h).model(j).fname = tmp;
subj(k).sess(h).model(j).F = DCM.DCM.F;
subj(k).sess(h).model(j).Ep = DCM.DCM.Ep;
subj(k).sess(h).model(j).Cp = DCM.DCM.Cp;
else
tmp = subj(k).sess(h).model(j).fname;
DCM = load(tmp);
F_sess = [F_sess,subj(k).sess(h).model(j).F];
end
% Data ID verification. At least for now we'll
% re-compute the IDs rather than use the ones stored
% with the DCM.
if ~exist('DCM','var') || ~isfield(DCM,'DCM'), job.verify_id = 0; end
if job.verify_id
disp(sprintf('Identifying data: model %d', j));
M = DCM.DCM.M;
str = spm_dcm_identify(DCM.DCM);
if strcmp(str,'fMRI') || strcmp(str,'fMRI_CSD')
Y = DCM.DCM.Y; % fMRI
else
Y = DCM.DCM.xY; % not fMRI
end
if isfield(M,'FS')
try
ID(j,h) = spm_data_id(feval(M.FS,Y.y,M));
catch
ID(j,h) = spm_data_id(feval(M.FS,Y.y));
end
else
ID(j,h) = spm_data_id(Y.y);
end
end
end
% Sum over sessions
F_mod = sum(F_sess);
F(k,j) = F_mod;
N{j} = sprintf('model%d',j);
end
if job.verify_id
failind = find(max(abs(diff(ID)),[],1) > eps);
if ~isempty(failind)
str_error = ['The models for subject ' num2str(k) ...
' session(s) ' num2str(failind) ' were not fitted to the same data.'];
spm('alert',str_error,'Error');
return
end
end
else
spm('alert*','The number of sessions/models should be the same for all subjects.','Error');
return
end
end
end
%-Bayesian Model Selection
%==========================================================================
% free energy
sumF = sum(F,1);
%-Family or model level
%--------------------------------------------------------------------------
if isfield(job.family_level,'family_file')
if ~isempty(job.family_level.family_file{1})
load(job.family_level.family_file{1});
do_family = 1;
family.infer = method;
nfam = size(family.names,2);
npart = length(unique(family.partition));
maxpart = max(family.partition);
m_indx = 1:nm;
if nfam ~= npart || npart == 1 || maxpart > npart
spm('alert*','Invalid family file.','Error');
return
end
else
do_family = 0;
end
else
if isempty(job.family_level.family)
do_family = 0;
else
do_family = 1;
nfam = size(job.family_level.family,2);
names_fam = {};
models_fam = [];
m_indx = [];
for f=1:nfam
names_fam = [names_fam,job.family_level.family(f).family_name];
m_indx = [m_indx,job.family_level.family(f).family_models'];
models_fam(job.family_level.family(f).family_models) = f;
end
family.names = names_fam;
family.partition = models_fam;
npart = length(unique(family.partition));
maxpart = max(family.partition);
nmodfam = length(m_indx);
if nfam ~= npart || npart == 1 || maxpart > npart || nmodfam > nm
spm('alert*','Invalid family.','Error');
return
end
family.infer = method;
end
end
%-Single subject BMS or 1st level (fixed effects) group BMS
%--------------------------------------------------------------------------
if strcmp(method,'FFX')
fprintf('%-40s: ','Computing FFX model/family posteriors') %-#
if ~do_family
family = [];
model.post = spm_api_bmc(sumF,N,[],[]);
else
Ffam = F(:,sort(m_indx));
[family,model] = spm_compare_families(Ffam,family);
P = spm_api_bmc(sumF,N,[],[],family);
end
fprintf('%30s\n','...done') %-#
if bma_do
if do_family
if do_bma_famwin
[fam_max,fam_max_i] = max(family.post);
indx = find(family.partition==fam_max_i);
post = model.post(indx);
else
if do_bma_all
post = model.post;
indx = 1:nm;
else
if bma_fam <= nfam && bma_fam > 0 && rem(bma_fam,1) == 0
indx = find(family.partition==bma_fam);
post = model.post(indx);
else
spm('alert*','Incorrect family for BMA.','Error');
return
end
end
end
else
post = model.post;
indx = 1:nm;
end
%-Bayesian Model Averaging
%------------------------------------------------------------------
fprintf('%-40s: ','FFX Bayesian model averaging') %-#
bma = spm_dcm_bma(post,indx,subj,nsamp,odds_ratio);
fprintf('%30s\n','...done') %-#
else
bma = [];
end
% if exist(fullfile(job.dir{1},'BMS.mat'),'file')
% load(fname);
% if isfield(BMS,'DCM') && isfield(BMS.DCM,'ffx')
% spm('alert','Existing BMS.mat file has been over-written!',...
% 'Warning');
% end
% end
BMS.DCM.ffx.data = fname_msp;
BMS.DCM.ffx.F_fname = f_fname;
BMS.DCM.ffx.F = F;
BMS.DCM.ffx.SF = sumF;
BMS.DCM.ffx.model = model;
BMS.DCM.ffx.family = family;
BMS.DCM.ffx.bma = bma;
% 2nd-level (random effects) BMS
%--------------------------------------------------------------------------
else
fprintf('%-40s: ','Computing RFX model/family posteriors') %-#
if ~do_family
[exp_r,xp,r_samp,g_post] = spm_BMS_gibbs(F);
model.g_post = g_post;
model.alpha = exp_r*(nm+ns);
model.exp_r = exp_r;
model.xp = xp;
family = [];
% Compute protected xp's
posterior.a=model.alpha;
posterior.r=g_post';
priors.a=ones(1,nm);
model.bor = spm_BMS_bor (F',posterior,priors);
% Compute protected exceedance probs - Eq 7 in Rigoux et al.
model.pxp=(1-model.bor)*xp+model.bor/nm;
else
Ffam = F(:,sort(m_indx));
[family,model] = spm_compare_families(Ffam,family);
end
fprintf('%30s\n','...done') %-#
%-Display result
if do_family
P = spm_api_bmc(sumF,N,model.exp_r,model.xp,family);
else
P = spm_api_bmc(sumF,N,model.exp_r,model.xp);
% Plot protected xp's
nm_tmp = length(N);
Fpxp = spm_figure('Create','Graphics','BMS: results ...');
figure(Fpxp);
subplot(2,1,1)
bar(1:length(N),model.pxp(:)')
set(gca,'XTick',1:length(N))
set(gca,'XTickLabel',1:nm_tmp)
ylabel('Protected Exceedance Probability','Fontsize',14)
xlabel('Models','Fontsize',14)
title(sprintf('Prob of Equal Model Frequencies (BOR) = %1.2f',model.bor),'Fontsize',14);
axis square
grid on
end
if bma_do
if do_family
if do_bma_famwin
[fam_max,fam_max_i] = max(family.exp_r);
indx = find(family.partition==fam_max_i);
post = model.g_post(:,indx);
else
if do_bma_all
post = model.g_post;
indx = 1:nm;
else
if bma_fam <= nfam && bma_fam > 0 && rem(bma_fam,1) == 0
indx = find(family.partition==bma_fam);
post = model.g_post(:,indx);
else
spm('alert','Incorrect family for BMA!','Error');
return
end
end
end
else
post = model.g_post;
indx = 1:nm;
end
%-Bayesian Model Averaging
%------------------------------------------------------------------
fprintf('%-40s: ','RFX Bayesian model averaging') %-#
bma = spm_dcm_bma(post,indx,subj,nsamp,odds_ratio);
fprintf('%30s\n','...done') %-#
else
bma = [];
end
% if exist(fullfile(job.dir{1},'BMS.mat'),'file')
% load(fname);
% if isfield(BMS,'DCM') && isfield(BMS.DCM,'rfx')
% spm('alert','Existing BMS.mat file has been over-written!',...
% 'Warning');
% end
% end
BMS.DCM.rfx.data = fname_msp;
BMS.DCM.rfx.F_fname = f_fname;
BMS.DCM.rfx.F = F;
BMS.DCM.rfx.SF = sumF;
BMS.DCM.rfx.model = model;
BMS.DCM.rfx.family = family;
BMS.DCM.rfx.bma = bma;
end
%-Data verification
%--------------------------------------------------------------------------
if job.verify_id
fprintf('%-40s: %30s\n','Data identity','...verified') %-#
else
fprintf('%-40s: %30s\n','Data identity','...not verified') %-#
end
%-Save BMS.mat
%--------------------------------------------------------------------------
fprintf('%-40s: %30s','Saving BMS.mat','...writing') %-#
save(fname,'BMS', spm_get_defaults('mat.format'));
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done') %-#
out.bmsmat = { fname };
%-Save model_space
%--------------------------------------------------------------------------
if ~ld_msp && data_se && ~ld_f
fprintf('%-40s: %30s','Saving model space','...writing') %-#
save(fname_msp,'subj', spm_get_defaults('mat.format'));
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done') %-#
out.model_space = { fname_msp };
end
c4science · Help