JiscMail Logo
Email discussion lists for the UK Education and Research communities

Help for SPM Archives


SPM Archives

SPM Archives


SPM@JISCMAIL.AC.UK


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

SPM Home

SPM Home

SPM  March 2007

SPM March 2007

Options

Subscribe or Unsubscribe

Subscribe or Unsubscribe

Log In

Log In

Get Password

Get Password

Subject:

Re: Question rephrased: is it possible to mask from a different statistical model

From:

Jan Gläscher <[log in to unmask]>

Reply-To:

Jan Gläscher <[log in to unmask]>

Date:

Thu, 1 Mar 2007 00:09:03 -0800

Content-Type:

multipart/mixed

Parts/Attachments:

Parts/Attachments

text/plain (88 lines) , spm_getSPM.m (881 lines)

Dear Chiara,

I ran into a similar problem and I modified the relevant SPM5 function 
(spm_getSPM) to include the feature that you are asking for. If you save 
the attached function in a directory in the Matlab path (before the SPM 
distribution), you can go ahead and press 'Results'. When you come to the
'mask with other contrast' section, and you press 'yes', then you can 
choose between 'current' or 'other' SPM analysis. Then you can select 
another SPM.mat and a masking contrast.

In addition (and unrelated to your question), my version of spm_getSPM also 
provides for an struct as an input parameter. In this struct you can 
specify the information about contrasts, thresholds etc. This way 
spm_getSPM will run non-interactively (e.g. in a batch mode) and create the 
usual xSPM struct with the list of voxel exceeding the contrast threshold. 
Please see the comments at the beginning of the function for more on this.

Hope this helps,
Jan

PS: If you need this for SPM2, let me know and I can email you a version 
that will run with SPM2.

Guillaume Sescousse wrote:
> Hi Chiara,
> I recently run into a similar problem, and the solution I found was to 
> use the PickAtlas toolbox (http://www.fmri.wfubmc.edu/download.htm) 
> which allows you to mask a T-map with any mask of your choice. The only 
> difference with the masking performed by SPM is that you have to create 
> your own mask from the desired contrast beforehand (using the "save" 
> option and then the "ImCalc" function to thresold and binarize your 
> saved image; you will find more helpful details in this post: 
> http://www.jiscmail.ac.uk/cgi-bin/webadmin?A2=ind04&L=SPM&P=R249758&D=0&I=-3&X=5562383E86C74C74). 
> 
> Hope it helps. All the best
> 
> Guillaume
> 
> Nosarti, Chiara a écrit :
>> Dear list
>> My question was asked in 1998 and at the time Karl stated that the
>> contrasts had to come from the same statistical model. Is this still the
>> case ?
>>
>> Best wishes
>> Chiara Nosarti
>>
>>
>> Date:         Tue, 22 Sep 1998 18:03:55 +0100
>> Reply-To:     [log in to unmask] (Karl Friston)
>> Sender:       spm
>> From:         [log in to unmask] (Karl Friston)
>> Subject:      Re: mask from another analysis? (was:Predicted Interval
>> Analysis
>>
>> Dear Ruth,
>>
>>  
>>> Dear Karl,      May I ask a further question on the subject of 
>>> masks.  Does the
>>>     
>> mask
>>  
>>> have to come from a contrast created in the very same analysis?  Or is
>>> it possible to "export" a mask from one analysis and use it as a mask
>>> with another set of scans?  Is there a way to save a mask as a file
>>>     
>> for
>>  
>>> later use?
>>>     
>>
>> The contrasts have to come from the same statistical model (i.e.
>> analysis).  This is not a mathematical requirement but is required by
>> the data reduction [to a point list] employed by SPM.
>>
>> All the best - Karl
>>
>>   
> 

-- 
Jan Gläscher, Ph.D.         Div. Humanities & Social Sciences
+1 (626) 395-4976 (office)  Caltech, Broad Center, M/C 114-96
+1 (626) 395-2000 (fax)     1200 California Blvd
[log in to unmask]    Pasadena, CA 91125



function [SPM,xSPM] = spm_getSPM(cfg) % computes a specified and thresholded SPM/PPM following parameter estimation % FORMAT [SPM,xSPM] = spm_getSPM; % % xSPM - structure containing SPM, distribution & filtering details % .swd - SPM working directory - directory containing current SPM.mat % .title - title for comparison (string) % .Z - minimum of Statistics {filtered on u and k} % .n - conjunction number <= number of contrasts % .STAT - distribution {Z, T, X, F or P} % .df - degrees of freedom [df{interest}, df{residual}] % .STATstr - description string % .Ic - indices of contrasts (in SPM.xCon) % .Im - indices of masking contrasts (in xCon) % .pm - p-value for masking (uncorrected) % .Ex - flag for exclusive or inclusive masking % .u - height threshold % .k - extent threshold {voxels} % .XYZ - location of voxels {voxel coords} % .XYZmm - location of voxels {mm} % .S - search Volume {voxels} % .R - search Volume {resels} % .FWHM - smoothness {voxels} % .M - voxels -> mm matrix % .iM - mm -> voxels matrix % .VOX - voxel dimensions {mm} - column vector % .DIM - image dimensions {voxels} - column vector % .Vspm - Mapped statistic image(s) % .Ps - list of P values for voxels at SPM.xVol.XYZ (used by FDR) % .thresDesc - description of height threshold (string) % % Required fields of SPM % % xVol - structure containing details of volume analysed % % xX - Design Matrix structure % - (see spm_spm.m for structure) % % xCon - Contrast definitions structure array % - (see also spm_FcUtil.m for structure, rules & handling) % .name - Contrast name % .STAT - Statistic indicator character ('T', 'F' or 'P') % .c - Contrast weights (column vector contrasts) % .X0 - Reduced design matrix data (spans design space under Ho) % Stored as coordinates in the orthogonal basis of xX.X from spm_sp % (Matrix in SPM99b) Extract using X0 = spm_FcUtil('X0',... % .iX0 - Indicates how contrast was specified: % If by columns for reduced design matrix then iX0 contains the % column indices. Otherwise, it's a string containing the % spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'} % .X1o - Remaining design space data (X1o is orthogonal to X0) % Stored as coordinates in the orthogonal basis of xX.X from spm_sp % (Matrix in SPM99b) Extract using X1o = spm_FcUtil('X1o',... % .eidf - Effective interest degrees of freedom (numerator df) % - Or effect-size threshold for Posterior probability % .Vcon - Name of contrast (for 'T's) or ESS (for 'F's) image % .Vspm - Name of SPM image % % In addition, the xCon.mat file is updated. For newly evaluated % contrasts, SPM images (spmT_????.{img,hdr}) are written, along with % contrast (con_????.{img,hdr}) images for SPM{T}'s, or Extra % Sum-of-Squares images (ess_????.{img,hdr}) for SPM{F}'s. % % The contrast images are the weighted sum of the parameter images, % where the weights are the contrast weights, and are uniquely % estimable since contrasts are checked for estimability by the % contrast manager. These contrast images (for appropriate contrasts) % are suitable summary images of an effect at this level, and can be % used as input at a higher level when effecting a random effects % analysis. (Note that the ess_????.{img,hdr} and % SPM{T,F}_????.{img,hdr} images are not suitable input for a higher % level analysis.) See spm_RandFX.man for further details. % %_______________________________________________________________________ % Non-interactive use of spm_getSPM % % The user may supply a config structure with the all the configuration % options. This will allow for a non-interactive creation of the % results-specific xSPM struct. This is useful for creating similar xSPM % structure for all subjects in a sample that can be then passed to % spm_regions for automatic data extraction (which can be also used % non-interactively). Additionally, the pop-up of windows may be during % the execution of spm_getSPM may be suppressed. % % The config structure has the following fields: % cfg.gui = flag for suppressing window output (1/0) % cfg.spmmatfile = full path and filename of the SPM.mat file % cfg.cnum = contrast number (can be determined in the contrast manager) % cfg.title = title for comparison % cfg.conjtype = different types of conjunctions: % 'conj','global','intermed' can be left empty of only % one contrast is specified in cfg.cnum % cfg.mask = flag to specify if the contrast should be masked with % another contrast (1/0) % cfg.omask = flag to specify if the masking contrast should be taken % from a different SPM.mat file (1/0) % cfg.mask_spmmatfile = full path and filename of SPM.mat from which the % masking contrast should be taken, can be left blank if % cfg.omask = 0; % cfg.mnum = number of masking contrast (can be determined from the % contrast manager of the SPM.mat which is used for masking % cfg.Pcon = p-value for thresholding the original contrast % cfg.Pthreshtype= type of p-value adjustment: 'FWE','FDR','none' % cfg.Extthresh = extent threshold (min. number of contiguous voxels) % cfg.Pmask = (uncorrected) p-value for masking contrast % cfg.masktype = 'exclusive' or 'inclusive' % cfg.ESthresh = only for PPMs: effect size threshold % cfg.PPMthresh = only for PPMs: posterior PPM threshold % cfg.plotES = only for PPMs: flag for plotting effect size (1) or % statistics (0) %_______________________________________________________________________ % % spm_getSPM prompts for an SPM and applies thresholds {u & k} % to a point list of voxel values (specified with their locations {XYZ}) % This allows the SPM be displayed and characterized in terms of regionally % significant effects by subsequent routines. % % For general linear model Y = XB + E with data Y, design matrix X, % parameter vector B, and (independent) errors E, a contrast c'B of the % parameters (with contrast weights c) is estimated by c'b, where b are % the parameter estimates given by b=pinv(X)*Y. % % Either single contrasts can be examined or conjunctions of different % contrasts. Contrasts are estimable linear combinations of the % parameters, and are specified using the SPM contrast manager % interface [spm_conman.m]. SPMs are generated for the null hypotheses % that the contrast is zero (or zero vector in the case of % F-contrasts). See the help for the contrast manager [spm_conman.m] % for a further details on contrasts and contrast specification. % % A conjunction assesses the conjoint expression of multiple effects. The % conjunction SPM is the minimum of the component SPMs defined by the % multiple contrasts. Inference on the minimum statistics can be % performed in different ways. Inference on the Conjunction Null (one or % more of the effects null) is accomplished by assessing the minimum as % if it were a single statistic; one rejects the conjunction null in % favor of the alternative that k=nc, that the number of active effects k % is equal to the number of contrasts nc. No assumptions are needed on % the dependence between the tests. % % Another approach is to make inference on the Global Null (all effects % null). Rejecting the Global Null of no (u=0) effects real implies an % alternative that k>0, that one or more effects are real. A third % Intermediate approach, is to use a null hypothesis of no more than u % effects are real. Rejecting the intermediate null that k<=u implies an % alternative that k>u, that more than u of the effects are real. % % The Global and Intermediate nulls use results for minimum fields which % require the SPMs to be identically distributed and independent. Thus, % all component SPMs must be either SPM{t}'s, or SPM{F}'s with the same % degrees of freedom. Independence is roughly guaranteed for large % degrees of freedom (and independent data) by ensuring that the % contrasts are "orthogonal". Note that it is *not* the contrast weight % vectors per se that are required to be orthogonal, but the subspaces of % the data space implied by the null hypotheses defined by the contrasts % (c'pinv(X)). Furthermore, this assumes that the errors are % i.i.d. (i.e. the estimates are maximum likelihood or Gauss-Markov. This % is the default in spm_spm). % % To ensure approximate independence of the component SPMs in the case of % the global or intermediate null, non-orthogonal contrasts are serially % orthogonalised in the order specified, possibly generating new % contrasts, such that the second is orthogonal to the first, the third % to the first two, and so on. Note that significant inference on the % global null only allows one to conclude that one or more of the effects % are real. Significant inference on the conjunction null allows one to % conclude that all of the effects are real. % % Masking simply eliminates voxels from the current contrast if they % do not survive an uncorrected p value (based on height) in one or % more further contrasts. No account is taken of this masking in the % statistical inference pertaining to the masked contrast. % % The SPM is subject to thresholding on the basis of height (u) and the % number of voxels comprising its clusters {k}. The height threshold is % specified as above in terms of an [un]corrected p value or % statistic. Clusters can also be thresholded on the basis of their % spatial extent. If you want to see all voxels simply enter 0. In this % instance the 'set-level' inference can be considered an 'omnibus test' % based on the number of clusters that obtain. % % BAYESIAN INFERENCE AND PPMS - POSTERIOR PROBABILITY MAPS % % If conditional estimates are available (and your contrast is a T % contrast) then you are asked whether the inference should be 'Bayesian' % or 'classical' (using GRF). If you choose Bayesian the contrasts are of % conditional (i.e. MAP) estimators and the inference image is a % posterior probability map (PPM). PPMs encode the probability that the % contrast exceeds a specified threshold. This threshold is stored in % the xCon.eidf. Subsequent plotting and tables will use the conditional % estimates and associated posterior or conditional probabilities. % % see spm_results_ui.m for further details of the SPM results section. % see also spm_contrasts.m %_______________________________________________________________________ % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience % Andrew Holmes, Karl Friston & Jean-Baptiste Poline % $Id: spm_getSPM.m 665 2006-10-24 11:51:52Z volkmar $ try cfg; catch cfg.gui = 1; end %-GUI setup %----------------------------------------------------------------------- if cfg.gui spm_help('!ContextHelp',mfilename) end %-Select SPM.mat & note SPM results directory %----------------------------------------------------------------------- try cfg.spmmatfile; if ~exist(cfg.spmmatfile); SPM = []; xSPM = []; return; end swd = spm_str_manip(cfg.spmmatfile,'H'); catch [spmmatfile, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat'); if ~sts, SPM = []; xSPM = []; return; end swd = spm_str_manip(spmmatfile,'H'); end %-Preliminaries... %======================================================================= %-Load SPM.mat %----------------------------------------------------------------------- try load(fullfile(swd,'SPM.mat')); catch error(['Cannot read ' fullfile(swd,'SPM.mat')]); end SPM.swd = swd; %-Get volumetric data from SPM.mat %======================================================================= % Dimensions: X: -68:68, Y: -100:72, Z: -42:82 DIM: [136; 172; 124] % adjust volumetrics if this is non-Talairach data %----------------------------------------------------------------------- global defaults % voxel-to-space mapping %----------------------------------------------------------------------- q = SPM.xVol.M - speye(4,4); % 3-D case %----------------------------------------------------------------------- if ~any(q(:))          % map x and y into anatomical space and make z a %     %-------------------------------------------------------------------     D = [136; 172; 100];     C = [68; 100; 0];     D = SPM.xVol.DIM./D;     C = D.*C;     iM = [D(1) 0 0 C(1);            0 D(2) 0 C(2);            0 0 D(3) C(3);            0 0 0 1];     SPM.xVol.iM = iM;     SPM.xVol.M = inv(iM);     % re-set units     %-------------------------------------------------------------------     defaults.units = {'mm' 'mm' '%'}; else     defaults.units = {'mm' 'mm' 'mm'}; end % 2-D case %----------------------------------------------------------------------- if strcmp(spm('CheckModality'), 'EEG') & SPM.xVol.DIM(3) == 1          defaults.units = {'mm' 'mm' ''};      end % get volumetrics %----------------------------------------------------------------------- try     if strcmp(spm('CheckModality'), 'EEG') & isfield(SPM.xX, 'fullrank')         Vbeta = SPM.Vbeta;     else         xX = SPM.xX; %-Design definition structure         XYZ = SPM.xVol.XYZ; %-XYZ coordinates         S = SPM.xVol.S; %-search Volume {voxels}         R = SPM.xVol.R; %-search Volume {resels}         M = SPM.xVol.M(1:3,1:3); %-voxels to mm matrix         VOX = sqrt(diag(M'*M))'; %-voxel dimensions     end catch      % check the model has been estimated %------------------------------------------------------------------- str = { 'This model has not been estimated.';... 'Would you like to estimate it now?'}; if spm_input(str,1,'bd','yes|no',[1,0],1) [SPM] = spm_spm(SPM); else return end end %-Contrast definitions %======================================================================= %-Load contrast definitions (if available) %----------------------------------------------------------------------- try xCon = SPM.xCon; catch xCon = {}; end %======================================================================= % - C O N T R A S T S , S P M C O M P U T A T I O N , M A S K I N G %======================================================================= %-Get contrasts %----------------------------------------------------------------------- try cfg.cnum; Ic = cfg.cnum; catch [Ic,xCon] = spm_conman(SPM,'T&F',Inf,... ' Select contrasts...',' for conjunction',1); end nc = length(Ic); % Number of contrasts %-Allow user to extend the null hypothesis for conjunctions % % n: conjunction number % u: Null hyp is k<=u effects real; Alt hyp is k>u effects real % (NB Here u is from Friston et al 2004 paper, not statistic thresh). % u n % Conjunction Null nc-1 1 | u = nc-n % Intermediate 1..nc-2 nc-u | #effects under null <= u % Global Null 0 nc | #effects under alt > u, >= u+1 %----------------------------------+------------------------------------- if (nc > 1) try cfg.conjtype; if strcmp(cfg.conjtype,'Conjunction') n = 1; elseif strcmp(cfg.conjtype,'Global') n = nc; elseif strcmp(cfg.conjtype,'Intermed') & nc > 2 n = NaN; else error('Unknown conjunction type in cfg.'); end catch if nc==2 But='Conjunction|Global'; Val=[1 nc]; else But='Conj''n|Intermed|Global'; Val=[1 NaN nc]; end n = spm_input('Null hyp. to assess?','+1','b',But,Val,1); if isnan(n) if nc==3, n = nc-1; else n = nc-spm_input('Effects under null ','0','n1','1',nc-1); end end end else n = 1; end %-Enforce orthogonality of multiple contrasts for conjunction % (Orthogonality within subspace spanned by contrasts) %----------------------------------------------------------------------- if nc>1 & n>1 & ~spm_FcUtil('|_?',xCon(Ic), xX.xKXs)     OrthWarn = 0;     %-Successively orthogonalise     %-NB: This loop is peculiarly controlled to account for the     % possibility that Ic may shrink if some contrasts diasppear     % on orthogonalisation (i.e. if there are colinearities)     %-------------------------------------------------------------------     i = 1; while(i < nc), i = i + 1;      %-Orthogonalise (subspace spanned by) contrast i wirit previous %--------------------------------------------------------------- oxCon = spm_FcUtil('|_',xCon(Ic(i)), xX.xKXs, xCon(Ic(1:i-1)));      %-See if this orthogonalised contrast has already been entered % or is colinear with a previous one. Define a new contrast if % neither is the case. %--------------------------------------------------------------- d = spm_FcUtil('In',oxCon,xX.xKXs,xCon); if spm_FcUtil('0|[]',oxCon,xX.xKXs) %-Contrast was colinear with a previous one - drop it %----------------------------------------------------------- Ic(i) = []; i = i - 1; elseif any(d) %-Contrast unchanged or already defined - note index %-----------------------------------------------------------             Ic(i) = min(d);         else             OrthWarn = OrthWarn + 1;             %-Define orthogonalised contrast as new contrast             %-----------------------------------------------------------             conlst = sprintf('%d,',Ic(1:i-1));             oxCon.name = sprintf('%s (orth. w.r.t {%s})', xCon(Ic(i)).name,...                                  conlst(1:end-1));             xCon = [xCon, oxCon];             Ic(i) = length(xCon);         end     end % while...          if OrthWarn       warning(sprintf(['Contrasts changed! %d contrasts orthogonalized ',... 'to allow conjunction inf.'], OrthWarn))     end     SPM.xCon = xCon; end % if nc>1... %-Get contrasts for masking %----------------------------------------------------------------------- try cfg.mask; catch cfg.mask = spm_input('mask with other contrast(s)','+1','y/n',[1,0],2); end if cfg.mask %- select analysis for masking contrasts %--------------------------------------------------------------- try cfg.omask; omask = cfg.omask; if omask otherSPM = load(cfg.mask_spmmatfile); mSPM = otherSPM.SPM; mxCon = mSPM.xCon; Im = cfg.mnum; else Im = cfg.mnum; end catch omask = spm_input('which analysis?','!+1','b','current|other',[0,1],1); if omask otherSPM = load(spm_select(1,'mat',... 'Select SPM.mat from other analysis',[],pwd,'^SPM.mat')); mSPM = otherSPM.SPM; [Im,mxCon] = spm_conman(mSPM,'T&F',-Inf,... 'Select contrasts for masking...',' for masking',1); else [Im,xCon] = spm_conman(SPM,'T&F',-Inf,... 'Select contrasts for masking...',' for masking',1); end end %-Threshold for mask (uncorrected p-value) %--------------------------------------------------------------- try cfg.Pmask; pm = cfg.Pmask; catch pm = spm_input('uncorrected mask p-value','+1','r',0.05,1,[0,1]); end %-Inclusive or exclusive masking %--------------------------------------------------------------- try cfg.masktype; if strcmp(cfg.masktype,'inclusive') Ex = 0; elseif strcmp(cfg.masktype,'exclusive') Ex = 1; else error('Unknown masktype in cfg.'); end catch Ex = spm_input('nature of mask','+1','b','inclusive|exclusive',[0,1]); end else Im = []; pm = []; Ex = []; end %-Create/Get title string for comparison %----------------------------------------------------------------------- if nc == 1 str = xCon(Ic).name; else str = [sprintf('contrasts {%d',Ic(1)),sprintf(',%d',Ic(2:end)),'}']; if n==nc str = [str ' (global null)']; elseif n==1 str = [str ' (conj. null)']; else str = [str sprintf(' (Ha: k>=%d)',(nc-n)+1)]; end end if Ex mstr = 'masked [excl.] by'; else mstr = 'masked [incl.] by'; end if length(Im) == 1 if exist('omask','var') & omask == 1 str = sprintf('%s (%s %s at p=%g)',str,mstr,mxCon(Im).name,pm); else str = sprintf('%s (%s %s at p=%g)',str,mstr,xCon(Im).name,pm); end elseif ~isempty(Im) str = [sprintf('%s (%s {%d',str,mstr,Im(1)),... sprintf(',%d',Im(2:end)),... sprintf('} at p=%g)',pm)]; end try cfg.title; titlestr = cfg.title; catch titlestr = spm_input('title for comparison','+1','s',str); end %-Bayesian or classical Inference? %----------------------------------------------------------------------- if isfield(SPM,'PPM')          % Make sure SPM.PPM.xCon field exists     if ~isfield(SPM.PPM,'xCon')         SPM.PPM.xCon=[];     end          % Set Bayesian con type     if length(SPM.PPM.xCon) < Ic        SPM.PPM.xCon(Ic).PSTAT = xCon(Ic).STAT;     end          % Make this one a Bayesian contrast     [xCon(Ic).STAT]=deal('P');          if all(strcmp([SPM.PPM.xCon(Ic).PSTAT],'T'))         % Simple contrast         str = 'Effect size threshold for PPM';         if isfield(SPM.PPM,'VB')             % For VB - set default effect size to zero             Gamma=0; try cfg.ppmthresh; xCon(ic).eidf = cfg.ppnthresh; catch xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma)); end         elseif nc == 1 & isempty(xCon(Ic).Vcon)             % con image not yet written             if spm_input('Inference',1,'b',{'Bayesian','classical'},[1 0]);                 %-Get Bayesian threshold (Gamma) stored in xCon(Ic).eidf                 % The default is one conditional s.d. of the contrast                 Gamma = sqrt(xCon(Ic).c'*SPM.PPM.Cb*xCon(Ic).c); try cfg.ESthresh; xCon(Ic).eidf = cfg.ESthresh catch xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma)); end             end         end     else         % Compound contrast using Chi^2 statistic         if ~isfield(xCon(Ic),'eidf')             xCon(Ic).eidf=0; % temporarily         end         if isempty(xCon(Ic).eidf)             xCon(Ic).eidf=0;         end     end end %-Compute & store contrast parameters, contrast/ESS images, & SPM images %======================================================================= SPM.xCon = xCon; if ~isfield(SPM.xX, 'fullrank')     if exist('omask','var') & omask == 1 SPM = spm_contrasts(SPM, unique([Ic])); %mSPM = spm_contrasts(mSPM,unique([Im])); else SPM = spm_contrasts(SPM, unique([Ic, Im])); end else     SPM = spm_eeg_contrasts_conv(SPM, unique([Ic, Im]));     xSPM = [];     return; end xCon = SPM.xCon; STAT = xCon(Ic(1)).STAT; VspmSv = cat(1,xCon(Ic).Vspm); %-Check conjunctions - Must be same STAT w/ same df %----------------------------------------------------------------------- if (nc > 1) & (any(diff(double(cat(1,xCon(Ic).STAT)))) | ... any(abs(diff(cat(1,xCon(Ic).eidf))) > 1)) error('illegal conjunction: can only conjoin SPMs of same STAT & df') end %-Degrees of Freedom and STAT string describing marginal distribution %----------------------------------------------------------------------- df = [xCon(Ic(1)).eidf xX.erdf]; if nc>1 if n>1 str = sprintf('^{%d \\{Ha:k\\geq%d\\}}',nc,(nc-n)+1); else str = sprintf('^{%d \\{Ha:k=%d\\}}',nc,(nc-n)+1); end else str = ''; end switch STAT case 'T' STATstr = sprintf('%c%s_{%.0f}','T',str,df(2)); case 'F' STATstr = sprintf('%c%s_{%.0f,%.0f}','F',str,df(1),df(2)); case 'P' STATstr = sprintf('%s^{%0.2f}','PPM',df(1)); end %-Compute (unfiltered) SPM pointlist for masked conjunction requested %======================================================================= fprintf('\t%-32s: %30s\n','SPM computation','...initialising') %-# %-Compute conjunction as minimum of SPMs %----------------------------------------------------------------------- Z = Inf; for i = Ic Z = min(Z,spm_get_data(xCon(i).Vspm,XYZ)); end % P values for False Discovery FDR rate computation (all search voxels) %======================================================================= switch STAT case 'T' Ps = (1 - spm_Tcdf(Z,df(2))).^n; case 'P' Ps = (1 - Z).^n; case 'F' Ps = (1 - spm_Fcdf(Z,df)).^n; end %-Compute mask and eliminate masked voxels %----------------------------------------------------------------------- for i = Im fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...masking') if exist('omask','var') & omask == 1 Vm = spm_vol(fullfile(mSPM.swd,mSPM.xCon(i).Vspm.fname)); Mask = spm_get_data(Vm,XYZ); um = spm_u(pm,[mSPM.xCon(i).eidf,mSPM.xX.erdf],mSPM.xCon(i).STAT); else Mask = spm_get_data(xCon(i).Vspm,XYZ); um = spm_u(pm,[xCon(i).eidf,xX.erdf],xCon(i).STAT); end if Ex Q = Mask <= um; else Q = Mask > um; end XYZ = XYZ(:,Q); Z = Z(Q); if isempty(Q) fprintf('\n') %-#      warning(sprintf('No voxels survive masking at p=%4.2f',pm)) break end end %-clean up interface %----------------------------------------------------------------------- fprintf('\t%-32s: %30s\n','SPM computation','...done') %-# try cfg.gui; if cfg.gui spm('Pointer','Arrow') end catch spm('Pointer','Arrow') end %======================================================================= % - H E I G H T & E X T E N T T H R E S H O L D S %======================================================================= %-Height threshold - classical inference %----------------------------------------------------------------------- u = -Inf; k = 0; if STAT ~= 'P'     %-Get height threshold     %-------------------------------------------------------------------     str = 'FWE|FDR|none';     % str = 'FWE|none'; % Use this line to disable FDR threshold try cfg.Pthreshtype; thresDesc = cfg.Pthreshtype; catch thresDesc = spm_input('p value adjustment to control','+1','b',str,[],1); end     switch thresDesc case 'FWE' % family-wise false positive rate     %--------------------------------------------------------------- try cfg.Pcon; u = cfg.Pcon; catch u = spm_input('p value (family-wise error)','+0','r',0.05,1,[0,1]); end     thresDesc = ['p<' num2str(u) ' (' thresDesc ')']; u = spm_uc(u,df,STAT,R,n,S); case 'FDR' % False discovery rate %--------------------------------------------------------------- try cfg.Pcon; u = cfg.Pcon; catch u = spm_input('p value (false discovery rate)','+0','r',0.05,1,[0,1]); end     thresDesc = ['p<' num2str(u) ' (' thresDesc ')']; u = spm_uc_FDR(u,df,STAT,n,VspmSv,0); otherwise %-NB: no adjustment % p for conjunctions is p of the conjunction SPM     %--------------------------------------------------------------- try cfg.Pcon; u = cfg.Pcon; catch u = spm_input(['threshold {',STAT,' or p value}'],'+0','r',0.001,1); end if u <= 1             thresDesc = ['p<' num2str(u) ' (unc.)'];             u = spm_u(u^(1/n),df,STAT);         else             thresDesc = [STAT '=' num2str(u) ];         end     end %-Height threshold - Bayesian inference %----------------------------------------------------------------------- elseif STAT == 'P'     u_default = 1- 1/SPM.xVol.S; try cfg.PPthresh; u = cfg.PPMthresh; catch u = spm_input(['Posterior probability threshold for PPM'],... '+0','r',u_default,1); end     thresDesc = ['P>' num2str(u) ' (PPM)']; end % (if STAT) %-Calculate height threshold filtering %------------------------------------------------------------------- Q = find(Z > u); %-Apply height threshold %------------------------------------------------------------------- Z = Z(:,Q); XYZ = XYZ(:,Q); if isempty(Q) warning(sprintf('No voxels survive height threshold u=%0.2g',u)) end %-Extent threshold (disallowed for conjunctions) %----------------------------------------------------------------------- if ~isempty(XYZ) & nc == 1     %-Get extent threshold [default = 0]     %------------------------------------------------------------------- try cfg.Extthresh; k = cfg.Extthresh; catch k = spm_input('& extent threshold {voxels}','+1','r',0,1,[0,Inf]); end     %-Calculate extent threshold filtering     %-------------------------------------------------------------------     A = spm_clusters(XYZ);     Q = [];     for i = 1:max(A)         j = find(A == i);         if length(j) >= k; Q = [Q j]; end     end     % ...eliminate voxels     %-------------------------------------------------------------------     Z = Z(:,Q);     XYZ = XYZ(:,Q);     if isempty(Q) warning(sprintf('No voxels survive extent threshold k=%0.2g',k))     end else     k = 0; end % (if ~isempty(XYZ)) %======================================================================= % - E N D %======================================================================= fprintf('\t%-32s: %30s\n','SPM computation','...done') %-# % For Bayesian inference provide (default) option to display contrast values if STAT=='P' try cfg.ESplot; pflag = cfg.ESplot; catch      pflag= spm_input('Plot effect-size/statistic',1,'b',{'Yes','No'},[1 0]); end if pflag         Z = spm_get_data(xCon(Ic).Vcon,XYZ);         %CC = spm_get_data(xCon(Ic).Vcon,XYZ);         %Z = CC(:,Q);     end end %-Assemble output structures of unfiltered data %======================================================================= xSPM = struct( ...         'swd', swd,... 'title', titlestr,... 'Z', Z,... 'n', n,... 'STAT', STAT,... 'df', df,... 'STATstr', STATstr,... 'Ic', Ic,... 'Im', Im,... 'pm', pm,... 'Ex', Ex,... 'u', u,... 'k', k,... 'XYZ', XYZ,... 'XYZmm', SPM.xVol.M(1:3,:)*[XYZ; ones(1,size(XYZ,2))],... 'S', SPM.xVol.S,... 'R', SPM.xVol.R,... 'FWHM', SPM.xVol.FWHM,... 'M', SPM.xVol.M,... 'iM', SPM.xVol.iM,... 'DIM', SPM.xVol.DIM,... 'VOX', VOX,... 'Vspm', VspmSv,... 'Ps', Ps,...         'thresDesc',thresDesc); % RESELS per voxel (density) if it exists %----------------------------------------------------------------------- if isfield(SPM,'VRpv'), xSPM.VRpv = SPM.VRpv; end

Top of Message | Previous Page | Permalink

JiscMail Tools


RSS Feeds and Sharing


Advanced Options


Archives

June 2017
May 2017
April 2017
March 2017
February 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
January 2016
December 2015
November 2015
October 2015
September 2015
August 2015
July 2015
June 2015
May 2015
April 2015
March 2015
February 2015
January 2015
December 2014
November 2014
October 2014
September 2014
August 2014
July 2014
June 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013
October 2013
September 2013
August 2013
July 2013
June 2013
May 2013
April 2013
March 2013
February 2013
January 2013
December 2012
November 2012
October 2012
September 2012
August 2012
July 2012
June 2012
May 2012
April 2012
March 2012
February 2012
January 2012
December 2011
November 2011
October 2011
September 2011
August 2011
July 2011
June 2011
May 2011
April 2011
March 2011
February 2011
January 2011
December 2010
November 2010
October 2010
September 2010
August 2010
July 2010
June 2010
May 2010
April 2010
March 2010
February 2010
January 2010
December 2009
November 2009
October 2009
September 2009
August 2009
July 2009
June 2009
May 2009
April 2009
March 2009
February 2009
January 2009
December 2008
November 2008
October 2008
September 2008
August 2008
July 2008
June 2008
May 2008
April 2008
March 2008
February 2008
January 2008
December 2007
November 2007
October 2007
September 2007
August 2007
July 2007
June 2007
May 2007
April 2007
March 2007
February 2007
January 2007
2006
2005
2004
2003
2002
2001
2000
1999
1998


WWW.JISCMAIL.AC.UK

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager