Dear Darren,
thanks for your fast response. However, the fix didn't work for me.
I tried to reestimate the design with you updated spm_spm_Bayes as well
as manually replacing "tmp = tmp + NaN*(~tmp); " with "tmp(~tmp) = NaN;" in
both lines
in my local spm_spm_Bayes. In both cases i get the same error when i test a
specific contrast
(3 conditions -1 0.5 0.5), however when i use a different contrast (e.g.,
condition 1 vs. implicit baseline: 1 0 0)
i don't get the error mesage and the results look reasonable. It seems as if
the error occurs whenever i try to contrast any of the modelled
conditions against one of the other modelled conditions, but not when i
estimate the design, Strangely, the same design gives reasonable results in
an other patients dataset and normally
first-level bayesian analysis works extremly well for my data-sets. You can
see the error returnded when i use dbstop below, to my limited understanding
the error starts in spm_Ncdf (the cum normal distribution..). I might add
that i use (still) MATLAB R13). Maybe you have any ideas left what's going
wrong?
Greetings,
Martin
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
In /opt/spm5/spm_input.m at line 1121
In /opt/spm5/spm_getSPM.m at line 376
In /opt/spm5/spm_results_ui.m at line 264
Warning in ==> /opt/spm5/spm_Ncdf.m
On line 80 ==> warning('Returning NaN for out of range arguments'), end
On Tue, 30 Jan 2007 10:03:24 -0600, d gitelman <[log in to unmask]>
wrote:
> 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') %-#
>
|