Dear Martin
Greetings. The problem is in spm_spm_Bayes.m when one uses Matlab R14 or
greater. If you go to approx lines 294 and 301 this code :
tmp = tmp + NaN*(~tmp);
produces the NaN's.
I have updated this in the FIL version of the file and hopefully it will
make its way to wider distribution at a future update (or it may have
already).
In any case replace the above with this line in both places.
tmp(~tmp) = NaN;
You will have to re-estimate the dataset, but then all should be fine.
Attached is also an updated version of spm_spm_Bayes.m
Regards,
Darren
> -----Original Message-----
> From: SPM (Statistical Parametric Mapping)
> [mailto:[log in to unmask]] On Behalf Of Martin Kronbichler
> Sent: Tuesday, January 30, 2007 4:12 AM
> To: [log in to unmask]
> Subject: [SPM] spm5 bayesian problems
>
> Dear SPM Experts,
>
> i have recently started to use variational bayesian first
> level analyses for analysing presurgical language and motor
> experiments because i had the feeling that the posterior
> probalities are easier to communicate and easier to
> understand (especially when making inferences on how likely
> increased activatity for speech or motor is in regions near
> to the lesions).
> Generally, the bayesian analyses worked extremly well.
> However i have stumpled across a strange result i do not
> understand when analysing a block-design language experiment
> (2 blockedc onditions: simple tones vs. auditory words
> presented in 2 sessions). Whatever effect size threshold
> (even completerly unrealistic values like > 1000) or
> posterior probality i choose the PPM shows the whole brains
> as exhibiting posterior prob. of 1 (i.e. the whole brain is
> "activated"). SPM5 displays the following error message when
> is estimate the contrasts:
>
> spm{P} image 1 : ...computing
> Warning: Returning NaN for out of range arguments
> > In /opt/spm5/spm_Ncdf.m at line 80
> In /opt/spm5/spm_contrasts.m at line 214
> In /opt/spm5/spm_getSPM.m at line 474
> In /opt/spm5/spm_results_ui.m at line 264
> ...written spmP_0001.img
> SPM computation :
> ...initialising
> SPM computation :
> ...done
> SPM computation :
> ...done
>
> I have tried to remove and add motion parameters, analyse
> only one session, use different basis functions etc... but
> whatever i do the problem persits.
> Inspection of the files generated by SPM5 shows that the beta
> and SDbeta images look reasonable. The spmP and con_SD images
> however contain only NaNs..
> I might add that a conventional analysis results in
> stastically highly signficant voxels...
> Anybody has any idea what is causing this problem?
>
> Greetings from Salzburg,
>
> Martin
>
>
> Martin Kronbichler, M.Sc.
> ---------------------------------------------
> Department of Psychology
>
>
> Center for Neurocognitive Research
> University of Salzburg
>
> Hellbrunnerstr.34, A 5020 Salzburg, Austria
> e-mail:[log in to unmask]
> Tel.:+43/(0)/662/8044-5162
>
> Fax:+43/(0)/662/8044-5126
>
> ---------------------------------------------
> Department of Neurology
> Center for Neurocognitive Research Christian-Doppler Clinic,
> Paracelsus Private Medical University Ignaz Harrerstr.79, A
> 5020 Salzburg, Austria e-mail:[log in to unmask]
> Tel.:+43/(0)/662/4483-3966
> ---------------------------------------------
>
function [SPM] = spm_spm_Bayes(SPM)
% Conditional parameter estimation of a General Linear Model
% FORMAT [SPM] = spm_spm_Bayes(SPM)
%_______________________________________________________________________
%
% For single subject fMRI analysis this function switches to
% using voxel-wise GLM-AR models that are spatially regularised
% using the VB framework. This is implemented using spm_spm_vb.m.
% One need not have previously analysed the data using classical estimation.
% For more on this option please see the help in that file.
% Otherwise, read on.
%
% spm_spm_Bayes returns to voxels identified by spm_spm (ML parameter
% estimation) to get conditional parameter estimates and ReML hyper-
% parameter estimates. These estimates use prior covariances, on the
% parameters, from emprical Bayes. These PEB prior variances come from
% the hierarchical model that obtains by considering voxels as providing
% a second level. Put simply, the variance in parameters, over voxels,
% is used as a prior variance from the point of view of any one voxel.
% The error covariance hyperparameters are re-estimated in the light of
% these priors. The approach adopted is essentially a fully Bayesian
% analysis at each voxel, using emprical Bayesian prior variance
% estimators over voxels.
%
% Each separable partition (i.e. session) is assigned its own
% hyperparameter but within session covariance components are lumped
% together, using their relative expectations over voxels. This makes
% things much more computationally efficient and avoids inefficient
% voxel-specific multiple hyperparameter estimates.
%
% spm_spm_Bayes adds the following fields to SPM:
%
% ----------------
%
%
% SPM.PPM.l = session-specific hyperparameter means
% SPM.PPM.Cb = empirical prior parameter covariances
% SPM.PPM.C = conditional covariances of parameters
% SPM.PPM.dC{i} = dC/dl;
% SPM.PPM.ddC{i} = ddC/dldl
%
% The derivatives are used to compute the conditional variance of various
% contrasts in spm_getSPM, using a Taylor expansion about the hyperparameter
% means.
%
%
% ----------------
%
% SPM.VCbeta - Handles of conditional parameter estimates
% SPM.VHp - Handles of hyperparameter estimates
%
% ----------------
%
% Cbeta_????.{img,hdr} - conditional parameter images
% These are 16-bit (float) images of the conditional estimates. The image
% files are numbered according to the corresponding column of the
% design matrix. Voxels outside the analysis mask (mask.img) are given
% value NaN.
%
% ----------------
%
% CHp_????.{img,hdr} - error covariance hyperparamter images
% This is a 32-bit (double) image of the ReML error variance estimate.
% for each separable partition (Session). Voxels outside the analysis
% mask are given value NaN.
%
%_______________________________________________________________________
% Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
% Karl Friston
% $Id: spm_spm_Bayes.m 587 2006-08-07 04:38:22Z Darren $
%-Say hello
%-----------------------------------------------------------------------
Finter = spm('FigName','Stats: Bayesian estimation...');
%-Select SPM.mat & change directory
%-----------------------------------------------------------------------
if ~nargin
swd = spm_str_manip(spm_select(1,'^SPM\.mat$','Select SPM.mat'),'H');
load(fullfile(swd,'SPM.mat'))
cd(swd)
end
try
M = SPM.xVol.M;
DIM = SPM.xVol.DIM;
xdim = DIM(1); ydim = DIM(2); zdim = DIM(3);
XYZ = SPM.xVol.XYZ;
catch
helpdlg({ 'Please do a ML estimation first',...
'This identifies the voxels to analyse'});
spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
return
end
%=======================================================================
% - A N A L Y S I S P R E L I M I N A R I E S
%=======================================================================
%-Initialise output images
%=======================================================================
fprintf('%-40s: %30s','Output images','...initialising') %-#
%-Intialise oonditional estimate image files
%-----------------------------------------------------------------------
xX = SPM.xX;
[nScan nBeta] = size(xX.X);
Vbeta(1:nBeta) = deal(struct(...
'fname', [],...
'dim', DIM',...
'dt', [spm_type('float32'), spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'descrip', ''));
for i = 1:nBeta
Vbeta(i).fname = sprintf('Cbeta_%04d.img',i);
Vbeta(i).descrip = sprintf('Cond. beta (%04d) - %s',i,xX.name{i});
spm_unlink(Vbeta(i).fname)
end
Vbeta = spm_create_vol(Vbeta);
%-Intialise ReML hyperparameter image files
%-----------------------------------------------------------------------
try
nHp = length(SPM.nscan);
catch
nHp = nScan;
SPM.nscan = nScan;
end
VHp(1:nHp) = deal(struct(...
'fname', [],...
'dim', DIM',...
'dt', [spm_type('float64'), spm_platform('bigend')],...
'mat', M,...
'pinfo', [1 0 0]',...
'descrip', ''));
for i = 1:nHp
VHp(i).fname = sprintf('Hp_%04d.img',i);
VHp(i).descrip = sprintf('Hyperparameter (%04d)',i);
spm_unlink(VHp(i).fname)
end
VHp = spm_create_vol(VHp);
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised') %-#
%=======================================================================
% - E M P I R I C A L B A Y E S F O R P R I O R V A R I A N C E
%=======================================================================
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...estimating priors') %-#
% get row u{i} and column v{i}/v0{i} indices for separable designs
%----------------------------------------------------------------------
s = nHp;
if isfield(SPM,'Sess')
for i = 1:s
u{i} = SPM.Sess(i).row;
v{i} = SPM.Sess(i).col;
v0{i} = xX.iB(i);
end
else
u{1} = [1:nScan];
v{1} = [xX.iH xX.iC];
v0{1} = [xX.iB xX.iG];
end
% cycle over separarable partitions
%-----------------------------------------------------------------------
for i = 1:s
% Get design X and confounds X0
%---------------------------------------------------------------
fprintf('%-30s- %i\n',' ReML Session',i);
X = xX.X(u{i}, v{i});
X0 = xX.X(u{i},v0{i});
[m n] = size(X);
% add confound in 'filter'
%---------------------------------------------------------------
if isstruct(xX.K)
X0 = full([X0 xX.K(i).X0]);
end
% orthogonalize X w.r.t. X0
%---------------------------------------------------------------
X = X - X0*(pinv(X0)*X);
% covariance components induced by parameter variations {Q}
%---------------------------------------------------------------
for j = 1:n
Q{j} = X*sparse(j,j,1,n,n)*X';
end
% covariance components induced by error non-sphericity {V}
%---------------------------------------------------------------
Q{n + 1} = SPM.xVi.V(u{i},u{i});
% ReML covariance component estimation
%---------------------------------------------------------------
[C h] = spm_reml(SPM.xVi.CY,X0,Q);
% check for negative variance components
%---------------------------------------------------------------
h = abs(h);
% 2-level model for this partition using prior variances sP(i)
% treat confounds as fixed (i.e. infinite prior variance)
%---------------------------------------------------------------
n0 = size(X0,2);
Cb = blkdiag(diag(h(1:n)),speye(n0,n0)*1e8);
P{1}.X = [X X0];
P{1}.C = {SPM.xVi.V};
P{2}.X = sparse(size(P{1}.X,2),1);
P{2}.C = Cb;
sP(i).P = P;
sP(i).u = u{:};
sP(i).v = [v{:} v0{:}];
end
%=======================================================================
% - F I T M O D E L & W R I T E P A R A M E T E R I M A G E S
%=======================================================================
%-Cycle to avoid memory problems (plane by plane)
%=======================================================================
spm_progress_bar('Init',100,'Bayesian estimation','');
helpdlg('This could take some time');
spm('Pointer','Watch')
%-maxMem is the maximum amount of data processed at a time (bytes)
%-----------------------------------------------------------------------
global defaults
MAXMEM = defaults.stats.maxmem;
blksz = ceil(MAXMEM/8/nScan);
SHp = 0; % sum of hyperparameters
for z = 1:zdim
% current plane-specific parameters
%-------------------------------------------------------------------
U = find(XYZ(3,:) == z);
nbch = ceil(length(U)/blksz);
CrBl = zeros(nBeta,length(U)); %-conditional parameter estimates
CrHp = zeros(nHp, length(U)); %-ReML hyperparameter estimates
for bch = 1:nbch %-loop over bunches of lines (planks)
%-construct list of voxels in this block
%---------------------------------------------------------------
I = [1:blksz] + (bch - 1)*blksz;
I = I(I <= length(U));
xyz = XYZ(:,U(I));
nVox = size(xyz,2);
%-Get response variable
%---------------------------------------------------------------
Y = spm_get_data(SPM.xY.VY,xyz);
%-Conditional estimates (per partition, per voxel)
%---------------------------------------------------------------
beta = zeros(nBeta,nVox);
Hp = zeros(nHp, nVox);
for j = 1:s
P = sP(j).P;
u = sP(j).u;
v = sP(j).v;
for i = 1:nVox
[C P] = spm_PEB(Y(u,i),P);
beta(v,i) = C{2}.E(1:length(v));
Hp(j,i) = C{1}.h;
end
end
%-Save for current plane in memory as we go along
%---------------------------------------------------------------
CrBl(:,I) = beta;
CrHp(:,I) = Hp;
SHp = SHp + sum(Hp,2);
end % (bch)
%-write out plane data to image files
%===================================================================
%-Write conditional beta images
%-------------------------------------------------------------------
for i = 1:nBeta
tmp = sparse(XYZ(1,U),XYZ(2,U),CrBl(i,:),xdim,ydim);
% The following line (and the same code on line 306 below) results in
% images with all NaN's (in Matlab R14 or greater). It has been
% commented out by DRG, and replaced by the line: tmp(~tmp) = NaN;
% tmp = tmp + NaN*(~tmp);
tmp(~tmp) = NaN;
Vbeta(i) = spm_write_plane(Vbeta(i),tmp,z);
end
%-Write hyperparameter images
%-------------------------------------------------------------------
for i = 1:nHp
tmp = sparse(XYZ(1,U),XYZ(2,U),CrHp(i,:),xdim,ydim);
% tmp = tmp + NaN*(~tmp);
tmp(~tmp) = NaN;
VHp(i) = spm_write_plane(VHp(i),tmp,z);
end
%-Report progress
%-------------------------------------------------------------------
spm_progress_bar('Set',100*(z - 1)/zdim);
end % (for z = 1:zdim)
fprintf('\n') %-#
spm_progress_bar('Clear')
%=======================================================================
% - P O S T E S T I M A T I O N
%=======================================================================
% Taylor expansion for conditional covariance
%-----------------------------------------------------------------------
fprintf('%-40s: %30s\n','Non-sphericity','...REML estimation') %-#
% expansion point (mean hyperparameters)
%-----------------------------------------------------------------------
l = SHp/SPM.xVol.S;
% change in conditional coavriance w.r.t. hyperparameters
%-----------------------------------------------------------------------
n = size(xX.X,2);
PPM.l = l;
for i = 1:s
PPM.dC{i} = sparse(n,n);
PPM.ddC{i} = sparse(n,n);
end
for i = 1:s
P = sP(i).P;
u = sP(i).u;
v = sP(i).v;
% derivatives of conditional covariance w.r.t. hyperparameters
%---------------------------------------------------------------
if spm_matlab_version_chk('6.5.1') < 0
% MATLAB 6.5.0 R13 has a bug in dealing with sparse matrices so revert to full matrices
fprintf('%-40s','MATLAB 6.5.0 R13: must use full matrices - this could cause memory problems...');
d = full(P{1}.X'*inv(P{1}.C{1})*P{1}.X);
Cby = full(inv(d/l(i) + inv(P{2}.C)));
d = full(d*Cby);
dC = full(Cby*d/(l(i)^2));
ddC = full(2*(dC/(l(i)^2) - Cby/(l(i)^3))*d);
else
% all other MATLAB versions: use sparse matrices
d = P{1}.X'*inv(P{1}.C{1})*P{1}.X;
Cby = inv(d/l(i) + inv(P{2}.C));
d = d*Cby;
dC = Cby*d/(l(i)^2);
ddC = 2*(dC/(l(i)^2) - Cby/(l(i)^3))*d;
end
% place in output structure
%---------------------------------------------------------------
if spm_matlab_version_chk('6.5.1') < 0,
% MATLAB 6.5.0 R13: revert to full matrices
j = 1:length(v);
PPM.Cb(v,v) = full(P{2}.C(j,j));
PPM.Cby(v,v) = full(Cby(j,j));
PPM.dC{i}(v,v) = full(dC(j,j));
PPM.ddC{i}(v,v) = full(ddC(j,j));
else
% all other MATLAB versions: use sparse matrices
j = 1:length(v);
PPM.Cb(v,v) = P{2}.C(j,j);
PPM.Cby(v,v) = Cby(j,j);
PPM.dC{i}(v,v) = dC(j,j);
PPM.ddC{i}(v,v) = ddC(j,j);
end
end
%-Save remaining results files and analysis parameters
%=======================================================================
fprintf('%-40s: %30s','Saving results','...writing') %-#
%-Save analysis parameters in SPM.mat file
%-----------------------------------------------------------------------
SPM.VCbeta = Vbeta; % Filenames - parameters
SPM.VHp = VHp; % Filenames - hyperparameters
SPM.PPM = PPM; % PPM structure
if spm_matlab_version_chk('7') >=0
save('SPM', 'SPM', '-V6');
else
save('SPM', 'SPM');
end;
fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done') %-#
%=======================================================================
%- E N D: Cleanup GUI
%=======================================================================
spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
fprintf('%-40s: %30s\n','Completed',spm('time')) %-#
fprintf('...use the results section for assessment\n\n') %-#
|