Hi Darren,
no the PPMs (all brain voxels) are still black whatever value i choose for
the effectsize threshold and/or the posterior proability.
On Wed, 31 Jan 2007 07:33:50 -0600, d gitelman <d-
[log in to unmask]> wrote:
>Dear Martin
>
>This is just a warning, not an error. Do the maps look correct now- i.e.,
>not totally black and the thresholds appear reasonable. If so then it's
>fine.
>
>darren
>
>> -----Original Message-----
>> From: Martin Kronbichler [mailto:[log in to unmask]]
>> Sent: Wednesday, January 31, 2007 4:18 AM
>> To: Darren G
>> Cc: Martin Kronbichler
>> Subject: Re: spm5 bayesian problems
>>
>> 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')
>> %-#
>> >
>>
>>
|