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 n 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)
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.
_______________________________________________________________________
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 one or more
effects. The conjunction SPM is the minimum of the component SPMs
defined by the multiple contrasts. The distributional results used
for minimum fields 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 a
conjunction, 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.
By default the conjunction p-value is computed for one or more
effects (i.e. k > 1). This implies the null hypothesis is k = 0.
There may be instances where you want to infer two or more effects
(e.g. in cognitive conjunctions) or more. This can be tested with
a null hypothesis that k = 1. You will be asked to specify the number
of effects present under the null. The tabulated p-values correspond
to an upper bound on the false positive rate for inferring more than
this number of effects were present. For a conventional conjunction
analysis with consistent contrasts enter 0. To infer a conjunction of
effects enter n - 1, where n is the number of contrasts
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
_______________________________________________________________________
@(#)spm_getSPM.m 2.54 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 04/08/040001 function [SPM,xSPM] = spm_getSPM 0002 % computes a specified and thresholded SPM/PPM following parameter estimation 0003 % FORMAT [SPM,xSPM] = spm_getSPM; 0004 % 0005 % xSPM - structure containing SPM, distribution & filtering details 0006 % .swd - SPM working directory - directory containing current SPM.mat 0007 % .title - title for comparison (string) 0008 % .Z - minimum of n Statistics {filtered on u and k} 0009 % .n - conjunction number <= number of contrasts 0010 % .STAT - distribution {Z, T, X, F or P} 0011 % .df - degrees of freedom [df{interest}, df{residual}] 0012 % .STATstr - description string 0013 % .Ic - indices of contrasts (in SPM.xCon) 0014 % .Im - indices of masking contrasts (in xCon) 0015 % .pm - p-value for masking (uncorrected) 0016 % .Ex - flag for exclusive or inclusive masking 0017 % .u - height threshold 0018 % .k - extent threshold {voxels} 0019 % .XYZ - location of voxels {voxel coords} 0020 % .XYZmm - location of voxels {mm} 0021 % .S - search Volume {voxels} 0022 % .R - search Volume {resels} 0023 % .FWHM - smoothness {voxels} 0024 % .M - voxels -> mm matrix 0025 % .iM - mm -> voxels matrix 0026 % .VOX - voxel dimensions {mm} - column vector 0027 % .DIM - image dimensions {voxels} - column vector 0028 % .Vspm - Mapped statistic image(s) 0029 % .Ps - list of P values for voxels at SPM.xVol.XYZ (used by FDR) 0030 % 0031 % Required fields of SPM 0032 % 0033 % xVol - structure containing details of volume analysed 0034 % 0035 % xX - Design Matrix structure 0036 % - (see spm_spm.m for structure) 0037 % 0038 % xCon - Contrast definitions structure array 0039 % - (see also spm_FcUtil.m for structure, rules & handling) 0040 % .name - Contrast name 0041 % .STAT - Statistic indicator character ('T', 'F' or 'P') 0042 % .c - Contrast weights (column vector contrasts) 0043 % .X0 - Reduced design matrix data (spans design space under Ho) 0044 % Stored as coordinates in the orthogonal basis of xX.X from spm_sp 0045 % (Matrix in SPM99b) Extract using X0 = spm_FcUtil('X0',... 0046 % .iX0 - Indicates how contrast was specified: 0047 % If by columns for reduced design matrix then iX0 contains the 0048 % column indices. Otherwise, it's a string containing the 0049 % spm_FcUtil 'Set' action: Usually one of {'c','c+','X0'} 0050 % .X1o - Remaining design space data (X1o is orthogonal to X0) 0051 % Stored as coordinates in the orthogonal basis of xX.X from spm_sp 0052 % (Matrix in SPM99b) Extract using X1o = spm_FcUtil('X1o',... 0053 % .eidf - Effective interest degrees of freedom (numerator df) 0054 % - Or effect-size threshold for Posterior probability 0055 % .Vcon - Name of contrast (for 'T's) or ESS (for 'F's) image 0056 % .Vspm - Name of SPM image 0057 % 0058 % In addition, the xCon.mat file is updated. For newly evaluated 0059 % contrasts, SPM images (spmT_????.{img,hdr}) are written, along with 0060 % contrast (con_????.{img,hdr}) images for SPM{T}'s, or Extra 0061 % Sum-of-Squares images (ess_????.{img,hdr}) for SPM{F}'s. 0062 % 0063 % The contrast images are the weighted sum of the parameter images, 0064 % where the weights are the contrast weights, and are uniquely 0065 % estimable since contrasts are checked for estimability by the 0066 % contrast manager. These contrast images (for appropriate contrasts) 0067 % are suitable summary images of an effect at this level, and can be 0068 % used as input at a higher level when effecting a random effects 0069 % analysis. (Note that the ess_????.{img,hdr} and 0070 % SPM{T,F}_????.{img,hdr} images are not suitable input for a higher 0071 % level analysis.) See spm_RandFX.man for further details. 0072 % 0073 %_______________________________________________________________________ 0074 % 0075 % spm_getSPM prompts for an SPM and applies thresholds {u & k} 0076 % to a point list of voxel values (specified with their locations {XYZ}) 0077 % This allows the SPM be displayed and characterized in terms of regionally 0078 % significant effects by subsequent routines. 0079 % 0080 % For general linear model Y = XB + E with data Y, design matrix X, 0081 % parameter vector B, and (independent) errors E, a contrast c'B of the 0082 % parameters (with contrast weights c) is estimated by c'b, where b are 0083 % the parameter estimates given by b=pinv(X)*Y. 0084 % 0085 % Either single contrasts can be examined or conjunctions of different 0086 % contrasts. Contrasts are estimable linear combinations of the 0087 % parameters, and are specified using the SPM contrast manager 0088 % interface [spm_conman.m]. SPMs are generated for the null hypotheses 0089 % that the contrast is zero (or zero vector in the case of 0090 % F-contrasts). See the help for the contrast manager [spm_conman.m] 0091 % for a further details on contrasts and contrast specification. 0092 % 0093 % A conjunction assesses the conjoint expression of one or more 0094 % effects. The conjunction SPM is the minimum of the component SPMs 0095 % defined by the multiple contrasts. The distributional results used 0096 % for minimum fields require the SPMs to be identically distributed and 0097 % independent. Thus, all component SPMs must be either SPM{t}'s, or 0098 % SPM{F}'s with the same degrees of freedom. Independence is roughly 0099 % guaranteed for large degrees of freedom (and independent data) by 0100 % ensuring that the contrasts are "orthogonal". Note that it is *not* 0101 % the contrast weight vectors per se that are required to be 0102 % orthogonal, but the subspaces of the data space implied by the null 0103 % hypotheses defined by the contrasts (c'pinv(X)). Furthermore, 0104 % this assumes that the errors are i.i.d. (i.e. the estimates are 0105 % maximum likelihood or Gauss-Markov. This is the default in spm_spm). 0106 % To ensure approximate independence of the component SPMs in a 0107 % conjunction, non-orthogonal contrasts are serially orthogonalised 0108 % in the order specified, possibly generating new contrasts, such that the 0109 % second is orthogonal to the first, the third to the first two, and so on. 0110 % 0111 % By default the conjunction p-value is computed for one or more 0112 % effects (i.e. k > 1). This implies the null hypothesis is k = 0. 0113 % There may be instances where you want to infer two or more effects 0114 % (e.g. in cognitive conjunctions) or more. This can be tested with 0115 % a null hypothesis that k = 1. You will be asked to specify the number 0116 % of effects present under the null. The tabulated p-values correspond 0117 % to an upper bound on the false positive rate for inferring more than 0118 % this number of effects were present. For a conventional conjunction 0119 % analysis with consistent contrasts enter 0. To infer a conjunction of 0120 % effects enter n - 1, where n is the number of contrasts 0121 % 0122 % Masking simply eliminates voxels from the current contrast if they 0123 % do not survive an uncorrected p value (based on height) in one or 0124 % more further contrasts. No account is taken of this masking in the 0125 % statistical inference pertaining to the masked contrast. 0126 % 0127 % The SPM is subject to thresholding on the basis of height (u) and the 0128 % number of voxels comprising its clusters {k}. The height threshold is 0129 % specified as above in terms of an [un]corrected p value or 0130 % statistic. Clusters can also be thresholded on the basis of their 0131 % spatial extent. If you want to see all voxels simply enter 0. In this 0132 % instance the 'set-level' inference can be considered an 'omnibus test' 0133 % based on the number of clusters that obtain. 0134 % 0135 % BAYESIAN INFERENCE AND PPMS - POSTERIOR PROBABILITY MAPS 0136 % 0137 % If conditional estimates are available (and your contrast is a T 0138 % contrast) then you are asked whether the inference should be 'Bayesian' 0139 % or 'classical' (using GRF). If you choose Bayesian the contrasts are of 0140 % conditional (i.e. MAP) estimators and the inference image is a 0141 % posterior probability map (PPM). PPMs encode the probability that the 0142 % contrast exceeds a specified threshold. This threshold is stored in 0143 % the xCon.eidf. Subsequent plotting and tables will use the conditional 0144 % estimates and associated posterior or conditional probabilities. 0145 % 0146 % see spm_results_ui.m for further details of the SPM results section. 0147 % see also spm_contrasts.m 0148 %_______________________________________________________________________ 0149 % @(#)spm_getSPM.m 2.54 Andrew Holmes, Karl Friston & Jean-Baptiste Poline 04/08/04 0150 0151 0152 SCCSid = '2.54'; 0153 0154 %-GUI setup 0155 %----------------------------------------------------------------------- 0156 SPMid = spm('SFnBanner',mfilename,SCCSid); 0157 spm_help('!ContextHelp',mfilename) 0158 0159 %-Select SPM.mat & note SPM results directory 0160 %----------------------------------------------------------------------- 0161 swd = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H'); 0162 0163 %-Preliminaries... 0164 %======================================================================= 0165 0166 %-Load SPM.mat 0167 %----------------------------------------------------------------------- 0168 load(fullfile(swd,'SPM.mat')); 0169 SPM.swd = swd; 0170 0171 %-Get volumetric data from SPM.mat 0172 %----------------------------------------------------------------------- 0173 try 0174 if strcmp(spm('CheckModality'), 'EEG') & rank(full(SPM.xX.X)) == size(SPM.xX.X, 1) 0175 Vbeta = SPM.Vbeta; 0176 else 0177 xX = SPM.xX; %-Design definition structure 0178 XYZ = SPM.xVol.XYZ; %-XYZ coordinates 0179 S = SPM.xVol.S; %-search Volume {voxels} 0180 R = SPM.xVol.R; %-search Volume {resels} 0181 M = SPM.xVol.M(1:3,1:3); %-voxels to mm matrix 0182 VOX = sqrt(diag(M'*M))'; %-voxel dimensions 0183 end 0184 catch 0185 0186 % check the model has been estimated 0187 %--------------------------------------------------------------- 0188 str = { 'This model has not been estimated.';... 0189 'Would you like to estimate it now?'}; 0190 if spm_input(str,1,'bd','yes|no',[1,0],1) 0191 [SPM] = spm_spm(SPM); 0192 else 0193 return 0194 end 0195 end 0196 0197 0198 %-Contrast definitions 0199 %======================================================================= 0200 0201 %-Load contrast definitions (if available) 0202 %----------------------------------------------------------------------- 0203 try 0204 xCon = SPM.xCon; 0205 catch 0206 xCon = {}; 0207 end 0208 0209 0210 %======================================================================= 0211 % - 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 0212 %======================================================================= 0213 0214 %-Get contrasts 0215 %----------------------------------------------------------------------- 0216 [Ic,xCon] = spm_conman(SPM,'T&F',Inf,... 0217 ' Select contrasts...',' for conjunction',1); 0218 0219 0220 %-Enforce orthogonality of multiple contrasts for conjunction 0221 % (Orthogonality within subspace spanned by contrasts) 0222 %----------------------------------------------------------------------- 0223 if length(Ic) > 1 & ~spm_FcUtil('|_?',xCon(Ic), xX.xKXs) 0224 0225 0226 %-Successively orthogonalise 0227 %-NB: This loop is peculiarly controlled to account for the 0228 % possibility that Ic may shrink if some contrasts diasppear 0229 % on orthogonalisation (i.e. if there are colinearities) 0230 %------------------------------------------------------------------- 0231 i = 1; while(i < length(Ic)), i = i + 1; 0232 0233 %-Orthogonalise (subspace spanned by) contrast i wirit previous 0234 %--------------------------------------------------------------- 0235 oxCon = spm_FcUtil('|_',xCon(Ic(i)), xX.xKXs, xCon(Ic(1:i-1))); 0236 0237 %-See if this orthogonalised contrast has already been entered 0238 % or is colinear with a previous one. Define a new contrast if 0239 % neither is the case. 0240 %--------------------------------------------------------------- 0241 d = spm_FcUtil('In',oxCon,xX.xKXs,xCon); 0242 0243 if spm_FcUtil('0|[]',oxCon,xX.xKXs) 0244 0245 %-Contrast was colinear with a previous one - drop it 0246 %----------------------------------------------------------- 0247 Ic(i) = []; 0248 i = i - 1; 0249 0250 elseif any(d) 0251 0252 %-Contrast unchanged or already defined - note index 0253 %----------------------------------------------------------- 0254 Ic(i) = min(d); 0255 0256 else 0257 0258 %-Define orthogonalised contrast as new contrast 0259 %----------------------------------------------------------- 0260 oxCon.name = [xCon(Ic(i)).name,' (orth. w.r.t {',... 0261 sprintf('%d,',Ic(1:i-2)), sprintf('%d})',Ic(i-1))]; 0262 xCon = [xCon, oxCon]; 0263 Ic(i) = length(xCon); 0264 end 0265 0266 end % while... 0267 0268 SPM.xCon = xCon; 0269 end % if length(Ic)... 0270 0271 0272 %-Get contrasts for masking 0273 %----------------------------------------------------------------------- 0274 if spm_input('mask with other contrast(s)','+1','y/n',[1,0],2) 0275 [Im,xCon] = spm_conman(SPM,'T&F',-Inf,... 0276 'Select contrasts for masking...',' for masking',1); 0277 0278 %-Threshold for mask (uncorrected p-value) 0279 %--------------------------------------------------------------- 0280 pm = spm_input('uncorrected mask p-value','+1','r',0.05,1,[0,1]); 0281 0282 %-Inclusive or exclusive masking 0283 %--------------------------------------------------------------- 0284 Ex = spm_input('nature of mask','+1','b','inclusive|exclusive',[0,1]); 0285 else 0286 Im = []; 0287 pm = []; 0288 Ex = []; 0289 end 0290 0291 0292 %-Create/Get title string for comparison 0293 %----------------------------------------------------------------------- 0294 if length(Ic) == 1 0295 str = xCon(Ic).name; 0296 else 0297 str = [sprintf('contrasts {%d',Ic(1)),sprintf(',%d',Ic(2:end)),'}']; 0298 end 0299 if Ex 0300 mstr = 'masked [excl.] by'; 0301 else 0302 mstr = 'masked [incl.] by'; 0303 end 0304 if length(Im) == 1 0305 str = sprintf('%s (%s %s at p=%g)',str,mstr,xCon(Im).name,pm); 0306 0307 elseif ~isempty(Im) 0308 str = [sprintf('%s (%s {%d',str,mstr,Im(1)),... 0309 sprintf(',%d',Im(2:end)),... 0310 sprintf('} at p=%g)',pm)]; 0311 end 0312 titlestr = spm_input('title for comparison','+1','s',str); 0313 0314 0315 0316 %-Bayesian or classical Inference? 0317 %----------------------------------------------------------------------- 0318 if isfield(SPM,'PPM') & xCon(Ic(1)).STAT == 'T' 0319 0320 if length(Ic) == 1 & isempty(xCon(Ic).Vcon) 0321 0322 if spm_input('Inference',1,'b',{'Bayesian','classical'},[1 0]); 0323 0324 % set STAT to 'P' 0325 %--------------------------------------------------------------- 0326 xCon(Ic).STAT = 'P'; 0327 0328 %--------------------------------------------------------------- 0329 str = 'threshold {default: prior s.d.}'; 0330 if SPM.PPM.VB==1 0331 % For VB Gamma is stored explicitly 0332 Gamma=SPM.PPM.Gamma; 0333 else 0334 %-Get Bayesian threshold (Gamma) stored in xCon(Ic).eidf 0335 % The default is one conditional s.d. of the contrast 0336 Gamma = sqrt(xCon(Ic).c'*SPM.PPM.Cb*xCon(Ic).c); 0337 end 0338 xCon(Ic).eidf = spm_input(str,'+1','e',sprintf('%0.2f',Gamma)); 0339 0340 end 0341 end 0342 end 0343 0344 0345 %-Compute & store contrast parameters, contrast/ESS images, & SPM images 0346 %======================================================================= 0347 SPM.xCon = xCon; 0348 if size(SPM.xX.X, 1) > rank(full(SPM.xX.X)) 0349 SPM = spm_contrasts(SPM, unique([Ic, Im])); 0350 else 0351 SPM = spm_eeg_contrasts_conv(SPM, unique([Ic, Im])); 0352 xSPM = []; 0353 return; 0354 end 0355 0356 xCon = SPM.xCon; 0357 VspmSv = cat(1,xCon(Ic).Vspm); 0358 STAT = xCon(Ic(1)).STAT; 0359 n = length(Ic); 0360 0361 %-Check conjunctions - Must be same STAT w/ same df 0362 %----------------------------------------------------------------------- 0363 if (n > 1) & (any(diff(double(cat(1,xCon(Ic).STAT)))) | ... 0364 any(abs(diff(cat(1,xCon(Ic).eidf))) > 1)) 0365 error('illegal conjunction: can only conjoin SPMs of same STAT & df') 0366 end 0367 0368 %-Allow user to extend the null hypothesis for conjunctions 0369 %----------------------------------------------------------------------- 0370 if (n > 1) 0371 n = n - spm_input('effects under null ','!+1','w1','0',n - 1); 0372 end 0373 0374 0375 %-Degrees of Freedom and STAT string describing marginal distribution 0376 %----------------------------------------------------------------------- 0377 df = [xCon(Ic(1)).eidf xX.erdf]; 0378 if n > 1 0379 str = sprintf('^{%d}',n); 0380 else 0381 str = ''; 0382 end 0383 0384 switch STAT 0385 case 'T' 0386 STATstr = sprintf('%c%s_{%.0f}','T',str,df(2)); 0387 case 'F' 0388 STATstr = sprintf('%c%s_{%.0f,%.0f}','F',str,df(1),df(2)); 0389 case 'P' 0390 STATstr = sprintf('%s^{%0.2f}','PPM',df(1)); 0391 end 0392 0393 0394 %-Compute (unfiltered) SPM pointlist for masked conjunction requested 0395 %======================================================================= 0396 fprintf('\t%-32s: %30s\n','SPM computation','...initialising') %-# 0397 0398 0399 %-Compute conjunction as minimum of SPMs 0400 %----------------------------------------------------------------------- 0401 Z = Inf; 0402 for i = Ic 0403 Z = min(Z,spm_get_data(xCon(i).Vspm,XYZ)); 0404 end 0405 0406 % P values for False Discovery FDR rate computation (all search voxels) 0407 %======================================================================= 0408 switch STAT 0409 case 'T' 0410 Ps = (1 - spm_Tcdf(Z,df(2))).^n; 0411 case 'P' 0412 Ps = (1 - Z).^n; 0413 case 'F' 0414 Ps = (1 - spm_Fcdf(Z,df)).^n; 0415 end 0416 0417 0418 %-Compute mask and eliminate masked voxels 0419 %----------------------------------------------------------------------- 0420 for i = Im 0421 fprintf('%s%30s',sprintf('\b')*ones(1,30),'...masking') 0422 0423 Mask = spm_get_data(xCon(i).Vspm,XYZ); 0424 um = spm_u(pm,[xCon(i).eidf,xX.erdf],xCon(i).STAT); 0425 if Ex 0426 Q = Mask <= um; 0427 else 0428 Q = Mask > um; 0429 end 0430 XYZ = XYZ(:,Q); 0431 Z = Z(Q); 0432 if isempty(Q) 0433 fprintf('\n') %-# 0434 warning(sprintf('No voxels survive masking at p=%4.2f',pm)) 0435 break 0436 end 0437 end 0438 0439 %-clean up interface 0440 %----------------------------------------------------------------------- 0441 fprintf('\t%-32s: %30s\n','SPM computation','...done') %-# 0442 spm('Pointer','Arrow') 0443 0444 0445 0446 %======================================================================= 0447 % - H E I G H T & E X T E N T T H R E S H O L D S 0448 %======================================================================= 0449 0450 %-Height threshold - classical inference 0451 %----------------------------------------------------------------------- 0452 u = -Inf; 0453 k = 0; 0454 if STAT ~= 'P' 0455 0456 0457 %-Get height threshold 0458 %------------------------------------------------------------------- 0459 str = 'FWE|FDR|none'; 0460 % str = 'FWE|none'; % Use this line to disable FDR threshold 0461 switch spm_input('p value adjustment to control','+1','b',str,[],1) 0462 0463 0464 case 'FWE' % family-wise false positive rate 0465 %--------------------------------------------------------------- 0466 u = spm_input('p value (family-wise error)','+0','r',0.05,1,[0,1]); 0467 u = spm_uc(u,df,STAT,R,n,S); 0468 0469 case 'FDR' % False discovery rate 0470 %--------------------------------------------------------------- 0471 u = spm_input('p value (false discovery rate)','+0','r',0.05,1,[0,1]); 0472 u = spm_uc_FDR(u,df,STAT,n,VspmSv,0); 0473 0474 otherwise %-NB: no adjustment 0475 % p for conjunctions is p of the conjunction SPM 0476 %--------------------------------------------------------------- 0477 u = spm_input(['threshold {',STAT,' or p value}'],'+0','r',0.001,1); 0478 if u <= 1; u = spm_u(u^(1/n),df,STAT); end 0479 0480 end 0481 0482 %-Height threshold - Bayesian inference 0483 %----------------------------------------------------------------------- 0484 elseif STAT == 'P' 0485 0486 u = spm_input(['p value threshold for PPM'],'+0','r',.95,1); 0487 0488 end % (if STAT) 0489 0490 %-Calculate height threshold filtering 0491 %------------------------------------------------------------------- 0492 Q = find(Z > u); 0493 0494 %-Apply height threshold 0495 %------------------------------------------------------------------- 0496 Z = Z(:,Q); 0497 XYZ = XYZ(:,Q); 0498 if isempty(Q) 0499 warning(sprintf('No voxels survive height threshold u=%0.2g',u)) 0500 end 0501 0502 0503 %-Extent threshold (disallowed for conjunctions) 0504 %----------------------------------------------------------------------- 0505 if ~isempty(XYZ) & length(Ic) == 1 0506 0507 %-Get extent threshold [default = 0] 0508 %------------------------------------------------------------------- 0509 k = spm_input('& extent threshold {voxels}','+1','r',0,1,[0,Inf]); 0510 0511 %-Calculate extent threshold filtering 0512 %------------------------------------------------------------------- 0513 A = spm_clusters(XYZ); 0514 Q = []; 0515 for i = 1:max(A) 0516 j = find(A == i); 0517 if length(j) >= k; Q = [Q j]; end 0518 end 0519 0520 % ...eliminate voxels 0521 %------------------------------------------------------------------- 0522 Z = Z(:,Q); 0523 XYZ = XYZ(:,Q); 0524 if isempty(Q) 0525 warning(sprintf('No voxels survive extent threshold k=%0.2g',k)) 0526 end 0527 0528 else 0529 0530 k = 0; 0531 0532 end % (if ~isempty(XYZ)) 0533 0534 0535 %======================================================================= 0536 % - E N D 0537 %======================================================================= 0538 fprintf('\t%-32s: %30s\n','SPM computation','...done') %-# 0539 0540 %-Assemble output structures of unfiltered data 0541 %======================================================================= 0542 xSPM = struct('swd', swd,... 0543 'title', titlestr,... 0544 'Z', Z,... 0545 'n', n,... 0546 'STAT', STAT,... 0547 'df', df,... 0548 'STATstr', STATstr,... 0549 'Ic', Ic,... 0550 'Im', Im,... 0551 'pm', pm,... 0552 'Ex', Ex,... 0553 'u', u,... 0554 'k', k,... 0555 'XYZ', XYZ,... 0556 'XYZmm', SPM.xVol.M(1:3,:)*[XYZ; ones(1,size(XYZ,2))],... 0557 'S', SPM.xVol.S,... 0558 'R', SPM.xVol.R,... 0559 'FWHM', SPM.xVol.FWHM,... 0560 'M', SPM.xVol.M,... 0561 'iM', SPM.xVol.iM,... 0562 'DIM', SPM.xVol.DIM,... 0563 'VOX', VOX,... 0564 'Vspm', VspmSv,... 0565 'Ps', Ps); 0566 0567 % RESELS per voxel (density) if it exists 0568 %----------------------------------------------------------------------- 0569 if isfield(SPM,'VRpv'), xSPM.VRpv = SPM.VRpv; end