VB estimation of voxel-specific GLM-AR models with spatial regularisation
FORMAT [SPM] = spm_spm_vb(SPM)
This function implements a VB estimation scheme for voxel-specific GLM-AR
models. Both regression coefficients and AR coefficients are spatially
regularised. The algorithm is described in a series of papers:
Paper VB1: W. Penny, S. Kiebel and K. Friston (2003) Variational Bayesian
Inference for fMRI time series. NeuroImage 19, pp 727-741.
Paper VB2: W.Penny, N. Trujillo-Bareto and K. Friston (2004). Bayesian fMRI
time series analysis with spatial priors. Submitted to NeuroImage.
Paper VB3: W.Penny (2004). Bayesian analysis of single-subject fMRI data:
SPM implementation. Technical Report. WDIN, UCL.
Paper VB4: W.Penny et al. (2004) Beyond the voxel: Comparing spatially extended
fMRI models. Manuscript in preparation.
The space to be analysed is a 'Volume', 'Masked Volume' or 'Slices'
For 'Slices' the numbers of the slices to be analysed are then entered
________________________________________________________________________________
Required fields of SPM:
xY.VY - nScan x 1 struct array of mapped image volumes
Images must have the same orientation, voxel size and data type
- Any scaling should have already been applied via the image handle
scalefactors.
xX - Structure containing design matrix information
- Required fields are:
xX.X - Design matrix (raw, not temporally smoothed)
xX.name - cellstr of parameter names corresponding to columns
of design matrix
xM - Structure containing masking information, or a simple column vector
of thresholds corresponding to the images in VY.
- If a structure, the required fields are:
xM.TH - nVar x nScan matrix of analysis thresholds, one per image
xM.I - Implicit masking (0=>none, 1 => implicit zero/NaN mask)
xM.VM - struct array of mapped explicit mask image volumes
- (empty if no explicit masks)
- Explicit mask images are >0 for valid voxels to assess.
- Mask images can have any orientation, voxel size or data
type. They are interpolated using nearest neighbour
interpolation to the voxel locations of the data Y.
- Note that voxels with constant data (i.e. the same value across
scans) are also automatically masked out.
________________________________________________________________________________
spm_spm_vb adds the following fields to SPM:
SPM.VCbeta - Handles of posterior parameter estimates (Cbeta_????)
SPM.VPsd - Handles of SD of posterior parameter estimates (SDbeta_????)
SPM.VHp - Handle of standard deviation of the error (SDerror)
SPM.VAR - Handles of AR coefficient images (AR_????)
SPM.PPM - Posterior Probability Map data structure
.VB=1, tells later functions (spm_contrasts, spm_graph)
that parameters were estimated using VB
.AR_P, assumed AR model order
.priors, type of priors used (eg. 'Spatial-GMRF')
see spm_vb_set_priors.m
.compute_F, whether model evidence was computed
.Gamma, default effect size threshold (used in spm_getSPM)
.M_X, set to 1 if an explicit mask is used (otherwise 0)
The following parameters are set if the space
to be analysed chosen as 'Slices'
.AN_slices, numbers of slices analysed
.slice(z), further info about GLM-AR model at slice z
eg. slice(z).F is evidence for slice z
(if computed)
These fields are used later in spm_contrasts.m, spm_graph.m
________________________________________________________________________________
The following images are written to file:
mask.{img,hdr} - analysis mask image
8-bit (uint8) image of zero-s & one's indicating which voxels were
included in the analysis. This mask image is the intersection of the
explicit, implicit and threshold masks specified in the xM argument.
The XYZ matrix contains the voxel coordinates of all voxels in the
analysis mask. The mask image is included for reference, but is not
explicitly used by the results section.
Note mask.img is only written if the selected space is 'Volume' or
'Masked Volume' (ie not 'Slices')
Cbeta_????.{img,hdr}
These are 16-bit (float) images of the parameter posteriors. 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.
SDbeta_????.{img,hdr}
These are 16-bit (float) images of the standard deviation of parameter
posteriors.
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.
SDerror.{img,hdr}
This is a 16-bit (float) image of the standard deviation of the error.
Voxels outside the analysis mask (mask.img) are given value NaN.
AR_????.{img,hdr}
This is a 16-bit (float) image of the AR coefficient.
The image files are numbered according to the order of the corresponding
AR coefficient.
Voxels outside the analysis mask (mask.img) are given value NaN.
----------------
%W% Will Penny %E%0001 function [SPM] = spm_spm_vb(SPM) 0002 % VB estimation of voxel-specific GLM-AR models with spatial regularisation 0003 % FORMAT [SPM] = spm_spm_vb(SPM) 0004 % 0005 % This function implements a VB estimation scheme for voxel-specific GLM-AR 0006 % models. Both regression coefficients and AR coefficients are spatially 0007 % regularised. The algorithm is described in a series of papers: 0008 % 0009 % Paper VB1: W. Penny, S. Kiebel and K. Friston (2003) Variational Bayesian 0010 % Inference for fMRI time series. NeuroImage 19, pp 727-741. 0011 % 0012 % Paper VB2: W.Penny, N. Trujillo-Bareto and K. Friston (2004). Bayesian fMRI 0013 % time series analysis with spatial priors. Submitted to NeuroImage. 0014 % 0015 % Paper VB3: W.Penny (2004). Bayesian analysis of single-subject fMRI data: 0016 % SPM implementation. Technical Report. WDIN, UCL. 0017 % 0018 % Paper VB4: W.Penny et al. (2004) Beyond the voxel: Comparing spatially extended 0019 % fMRI models. Manuscript in preparation. 0020 % 0021 % The space to be analysed is a 'Volume', 'Masked Volume' or 'Slices' 0022 % For 'Slices' the numbers of the slices to be analysed are then entered 0023 % 0024 % ________________________________________________________________________________ 0025 % 0026 % Required fields of SPM: 0027 % 0028 % xY.VY - nScan x 1 struct array of mapped image volumes 0029 % Images must have the same orientation, voxel size and data type 0030 % - Any scaling should have already been applied via the image handle 0031 % scalefactors. 0032 % 0033 % xX - Structure containing design matrix information 0034 % - Required fields are: 0035 % xX.X - Design matrix (raw, not temporally smoothed) 0036 % xX.name - cellstr of parameter names corresponding to columns 0037 % of design matrix 0038 % 0039 % xM - Structure containing masking information, or a simple column vector 0040 % of thresholds corresponding to the images in VY. 0041 % - If a structure, the required fields are: 0042 % xM.TH - nVar x nScan matrix of analysis thresholds, one per image 0043 % xM.I - Implicit masking (0=>none, 1 => implicit zero/NaN mask) 0044 % xM.VM - struct array of mapped explicit mask image volumes 0045 % - (empty if no explicit masks) 0046 % - Explicit mask images are >0 for valid voxels to assess. 0047 % - Mask images can have any orientation, voxel size or data 0048 % type. They are interpolated using nearest neighbour 0049 % interpolation to the voxel locations of the data Y. 0050 % - Note that voxels with constant data (i.e. the same value across 0051 % scans) are also automatically masked out. 0052 % 0053 % ________________________________________________________________________________ 0054 % 0055 % spm_spm_vb adds the following fields to SPM: 0056 % 0057 % SPM.VCbeta - Handles of posterior parameter estimates (Cbeta_????) 0058 % SPM.VPsd - Handles of SD of posterior parameter estimates (SDbeta_????) 0059 % SPM.VHp - Handle of standard deviation of the error (SDerror) 0060 % SPM.VAR - Handles of AR coefficient images (AR_????) 0061 % 0062 % SPM.PPM - Posterior Probability Map data structure 0063 % 0064 % .VB=1, tells later functions (spm_contrasts, spm_graph) 0065 % that parameters were estimated using VB 0066 % 0067 % .AR_P, assumed AR model order 0068 % 0069 % .priors, type of priors used (eg. 'Spatial-GMRF') 0070 % see spm_vb_set_priors.m 0071 % 0072 % .compute_F, whether model evidence was computed 0073 % 0074 % .Gamma, default effect size threshold (used in spm_getSPM) 0075 % 0076 % .M_X, set to 1 if an explicit mask is used (otherwise 0) 0077 % 0078 % The following parameters are set if the space 0079 % to be analysed chosen as 'Slices' 0080 % .AN_slices, numbers of slices analysed 0081 % .slice(z), further info about GLM-AR model at slice z 0082 % eg. slice(z).F is evidence for slice z 0083 % (if computed) 0084 % These fields are used later in spm_contrasts.m, spm_graph.m 0085 % 0086 % ________________________________________________________________________________ 0087 % 0088 % The following images are written to file: 0089 % 0090 % mask.{img,hdr} - analysis mask image 0091 % 8-bit (uint8) image of zero-s & one's indicating which voxels were 0092 % included in the analysis. This mask image is the intersection of the 0093 % explicit, implicit and threshold masks specified in the xM argument. 0094 % The XYZ matrix contains the voxel coordinates of all voxels in the 0095 % analysis mask. The mask image is included for reference, but is not 0096 % explicitly used by the results section. 0097 % 0098 % Note mask.img is only written if the selected space is 'Volume' or 0099 % 'Masked Volume' (ie not 'Slices') 0100 % 0101 % Cbeta_????.{img,hdr} 0102 % These are 16-bit (float) images of the parameter posteriors. The image 0103 % files are numbered according to the corresponding column of the 0104 % design matrix. Voxels outside the analysis mask (mask.img) are given 0105 % value NaN. 0106 % 0107 % SDbeta_????.{img,hdr} 0108 % These are 16-bit (float) images of the standard deviation of parameter 0109 % posteriors. 0110 % The image files are numbered according to the corresponding column of the 0111 % design matrix. Voxels outside the analysis mask (mask.img) are given 0112 % value NaN. 0113 % 0114 % SDerror.{img,hdr} 0115 % This is a 16-bit (float) image of the standard deviation of the error. 0116 % Voxels outside the analysis mask (mask.img) are given value NaN. 0117 % 0118 % AR_????.{img,hdr} 0119 % This is a 16-bit (float) image of the AR coefficient. 0120 % The image files are numbered according to the order of the corresponding 0121 % AR coefficient. 0122 % Voxels outside the analysis mask (mask.img) are given value NaN. 0123 % 0124 % ---------------- 0125 % 0126 % %W% Will Penny %E% 0127 0128 % Let later functions know (eg. spm_contrasts) that 0129 % estimation was with VB 0130 PPM.VB=1; 0131 0132 %-Say hello 0133 %----------------------------------------------------------------------- 0134 Finter = spm('FigName','Stats: Bayesian Specification ...'); 0135 spm('Pointer','Arrow'); 0136 0137 %-Get SPM.mat if necessary 0138 %----------------------------------------------------------------------- 0139 if nargin ==0 0140 swd = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H'); 0141 load(fullfile(swd,'SPM.mat')); 0142 SPM.swd = swd; 0143 end 0144 0145 %-Change to SPM.swd if specified 0146 %----------------------------------------------------------------------- 0147 try 0148 cd(SPM.swd); 0149 end 0150 0151 %-Ensure data are assigned 0152 %----------------------------------------------------------------------- 0153 try 0154 SPM.xY.VY; 0155 catch 0156 helpdlg({ 'Please assign data to this design',... 0157 'Use fMRI under model specification'}); 0158 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') 0159 return 0160 end 0161 0162 %-Delete files from previous analyses 0163 %----------------------------------------------------------------------- 0164 if exist(fullfile('.','mask.img'),'file') == 2 0165 0166 str = {'Current directory contains SPM estimation files:',... 0167 'pwd = ',pwd,... 0168 'Existing results will be overwritten!'}; 0169 0170 abort = spm_input(str,1,'bd','stop|continue',[1,0],1); 0171 if abort 0172 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') 0173 return 0174 else 0175 str = sprintf('Overwriting old results\n\t (pwd = %s) ',pwd) 0176 warning(str) 0177 drawnow 0178 end 0179 end 0180 0181 files = { 'mask.???','ResMS.???','RVP.???',... 0182 'beta_????.???','con_????.???','ResI_????.???',... 0183 'ess_????.???', 'spm?_????.???'}; 0184 0185 for i=1:length(files) 0186 if any(files{i} == '*'|files{i} == '?' ) 0187 [j,null] = spm_list_files(pwd,files{i}); 0188 for i=1:size(j,1) 0189 spm_unlink(deblank(j(i,:))) 0190 end 0191 else 0192 spm_unlink(files{i}) 0193 end 0194 end 0195 0196 %======================================================================= 0197 % - A N A L Y S I S P R E L I M I N A R I E S 0198 %======================================================================= 0199 0200 0201 0202 %-Initialise output images 0203 %======================================================================= 0204 fprintf('%-40s: %30s','Output images','...initialising') %-# 0205 0206 VY = SPM.xY.VY; 0207 M = VY(1).mat; 0208 DIM = VY(1).dim(1:3)'; 0209 xdim = DIM(1); ydim = DIM(2); zdim = DIM(3); 0210 0211 %-Initialise XYZ matrix of in-mask voxel co-ordinates (real space) 0212 %----------------------------------------------------------------------- 0213 XYZ = zeros(3,xdim*ydim*zdim); 0214 0215 %-Intialise conditional estimate image files 0216 %----------------------------------------------------------------------- 0217 xX = SPM.xX; 0218 [nScan nBeta] = size(xX.X); 0219 Vbeta(1:nBeta) = deal(struct(... 0220 'fname', [],... 0221 'dim', [DIM',spm_type('float')],... 0222 'mat', M,... 0223 'pinfo', [1 0 0]',... 0224 'descrip', '')); 0225 for i = 1:nBeta 0226 Vbeta(i).fname = sprintf('Cbeta_%04d.img',i); 0227 Vbeta(i).descrip = sprintf('Posterior mean of beta (%04d) - %s',i,xX.name{i}); 0228 spm_unlink(Vbeta(i).fname) 0229 end 0230 Vbeta = spm_create_vol(Vbeta,'noopen'); 0231 0232 %-Intialise hyperparameter (AR 1..p and noise variance) image files 0233 %----------------------------------------------------------------------- 0234 0235 % Set number of AR coefficients 0236 str=' AR model order'; 0237 PPM.AR_P = spm_input(str,1,'e',0,[1 1]); 0238 0239 % Specify type of prior 0240 spm_input('Overall Priors ...',1,'d'); 0241 Ctype = { 0242 'Spatial - GMRF',... 0243 'Spatial - LORETA',... 0244 'Voxel - Shrinkage',... 0245 'Voxel - Uninformative'}; 0246 str = 'Select prior'; 0247 Sel = spm_input(str,2,'m',Ctype); 0248 PPM.priors.WA = Ctype{Sel}; 0249 0250 % Override option for prior on A - none 0251 PPM.priors.overA='none'; 0252 0253 % Compute evidence ? 0254 str = 'Compute evidence ?'; 0255 PPM.compute_F = spm_input(str,1,'b',{'yes','no'},[1 0]); 0256 0257 % Set number of VB iterations 0258 %str='Number of VB iterations'; 0259 %maxVBits = spm_input(str,1,'e',0,[1 1]); 0260 0261 % Initialise Error SD image 0262 VHp = deal(struct(... 0263 'fname', [],... 0264 'dim', [DIM',spm_type('double')],... 0265 'mat', M,... 0266 'pinfo', [1 0 0]',... 0267 'descrip', '')); 0268 0269 VHp.fname = sprintf('SDerror.img'); 0270 VHp.descrip = sprintf('Error SD'); 0271 spm_unlink(VHp.fname) 0272 VHp = spm_create_vol(VHp,'noopen'); 0273 0274 % Initialise Posterior SD images 0275 nPsd=size(SPM.xX.X,2); 0276 VPsd(1:nPsd) = deal(struct(... 0277 'fname', [],... 0278 'dim', [DIM',spm_type('double')],... 0279 'mat', M,... 0280 'pinfo', [1 0 0]',... 0281 'descrip', '')); 0282 for i = 1:nPsd 0283 VPsd(i).fname = sprintf('SDbeta_%04d.img',i); 0284 VPsd(i).descrip = sprintf('Posterior SD of beta (%04d)',i); 0285 spm_unlink(VPsd(i).fname) 0286 end 0287 VPsd = spm_create_vol(VPsd,'noopen'); 0288 0289 % Initialise AR images 0290 VAR(1:PPM.AR_P) = deal(struct(... 0291 'fname', [],... 0292 'dim', [DIM',spm_type('double')],... 0293 'mat', M,... 0294 'pinfo', [1 0 0]',... 0295 'descrip', '')); 0296 for i = 1:PPM.AR_P 0297 VAR(i).fname = sprintf('AR_%04d.img',i); 0298 VAR(i).descrip = sprintf('Autoregressive coefficient (%04d)',i); 0299 spm_unlink(VAR(i).fname) 0300 end 0301 VAR = spm_create_vol(VAR,'noopen'); 0302 0303 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...initialised'); %-# 0304 0305 % Set up masking details 0306 %-If xM is not a structure then assumme it's a vector of thresholds 0307 %----------------------------------------------------------------------- 0308 try 0309 xM = SPM.xM; 0310 catch 0311 xM = -ones(nScan,1)/0; 0312 end 0313 0314 % Analyse volume/masked volume/slices ? 0315 spm_input('Data to analyse...',1,'d') 0316 Ctype = { 0317 'Volume',... 0318 'Masked Volume',... 0319 'Slices'}; 0320 str = 'Select space'; 0321 Sel = spm_input(str,2,'m',Ctype); 0322 space_type = Ctype{Sel}; 0323 0324 switch space_type, 0325 case 'Masked Volume', 0326 PPM.M_X=1; 0327 PPM.AN_slices=[1:1:zdim]; 0328 case 'Volume', 0329 PPM.M_X=0; 0330 PPM.AN_slices=[1:1:zdim]; 0331 case 'Slices', 0332 PPM.M_X=0; 0333 PPM.AN_slices = spm_input(['Enter slice numbers eg. 3 14 2'],'+1'); 0334 end 0335 0336 0337 % Explicit masking ? 0338 %M_X = spm_input('explicitly mask images?','+1','y/n',[1,0],2); 0339 0340 if PPM.M_X, 0341 M_P = spm_get(Inf,'*.img',{'select mask images'}); 0342 else, 0343 M_P = {}; 0344 end 0345 0346 if ~isempty(M_P) 0347 xM.VM = spm_vol(char(M_P)); 0348 end 0349 0350 if ~isstruct(xM) 0351 xM = struct( 'T', [],... 0352 'TH', xM,... 0353 'I', 0,... 0354 'VM', {[]},... 0355 'xs', struct('Masking','analysis threshold')); 0356 end 0357 0358 %-Intialise the name of the new mask : current mask & conditions on voxels 0359 %----------------------------------------------------------------------- 0360 VM = struct( 'fname', 'mask.img',... 0361 'dim', [DIM',spm_type('uint8')],... 0362 'mat', M,... 0363 'pinfo', [1 0 0]',... 0364 'descrip', 'spm_spm:resultant analysis mask'); 0365 VM = spm_create_vol(VM); 0366 0367 0368 %======================================================================= 0369 % - 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 0370 %======================================================================= 0371 0372 %-Cycle to avoid memory problems (plane by plane) 0373 %======================================================================= 0374 spm_progress_bar('Init',100,'VB estimation',''); 0375 spm('Pointer','Watch') 0376 0377 %-maxMem is the maximum amount of data processed at a time (bytes) 0378 %----------------------------------------------------------------------- 0379 global defaults 0380 MAXMEM = defaults.stats.maxmem*10; 0381 blksz = ceil(MAXMEM/8/nScan); 0382 nbch = ceil(xdim*ydim/blksz); %-# blocks 0383 0384 SHp=0; % Sum of noise variance hyperparameters 0385 xords = [1:xdim]'*ones(1,ydim); xords = xords(:)'; % plane X coordinates 0386 yords = ones(xdim,1)*[1:ydim]; yords = yords(:)'; % plane Y coordinates 0387 S = 0; % Number of in-mask voxels 0388 s = 0; % Volume (voxels > UF) 0389 0390 % Initialise aspects of slice variables common to all slices 0391 slice_template = spm_vb_init_volume (SPM.xX.X,PPM.AR_P); 0392 0393 index=1; 0394 for z = 1:zdim, 0395 0396 % current plane-specific parameters 0397 %------------------------------------------------------------------- 0398 zords = z*ones(xdim*ydim,1)'; %-plane Z coordinates 0399 CrBl = []; %-conditional parameter estimates 0400 CrHp = []; % VB hyperparameter estimates 0401 CrPsd = []; % 0402 CrAR = []; % 0403 Q = []; %-in mask indices for this plane 0404 0405 analyse_slice = length(find((PPM.AN_slices-z)==0)); 0406 if analyse_slice 0407 for bch = 1:nbch %-loop over bunches of lines (planks) 0408 0409 %-# Print progress information in command window 0410 %--------------------------------------------------------------- 0411 str = sprintf('Plane %3d/%-3d, block %3d/%-3d',z,zdim,bch,nbch); 0412 fprintf('\r%-40s: %30s',str,' ') %-# 0413 0414 %-construct list of voxels in this block 0415 %--------------------------------------------------------------- 0416 I = [1:blksz] + (bch - 1)*blksz; 0417 I = I(I <= xdim*ydim); %-truncate 0418 xyz = [xords(I); yords(I); zords(I)]; %-voxel coordinates 0419 nVox = size(xyz,2); 0420 0421 %-Get data & construct analysis mask 0422 %=============================================================== 0423 fprintf('%s%30s',sprintf('\b')*ones(1,30),'...read & mask data')%-# 0424 Cm = logical(ones(1,nVox)); %-current mask 0425 0426 %-Compute explicit mask 0427 % (note that these may not have same orientations) 0428 %--------------------------------------------------------------- 0429 for i = 1:length(xM.VM) 0430 0431 %-Coordinates in mask image 0432 %------------------------------------------------------- 0433 j = xM.VM(i).mat\M*[xyz;ones(1,nVox)]; 0434 0435 %-Load mask image within current mask & update mask 0436 %------------------------------------------------------- 0437 Cm(Cm) = spm_get_data(xM.VM(i),j(:,Cm)) > 0; 0438 end 0439 0440 %-Get the data in mask, compute threshold & implicit masks 0441 %--------------------------------------------------------------- 0442 Y = zeros(nScan,nVox); 0443 for i = 1:nScan 0444 0445 %-Load data in mask 0446 %------------------------------------------------------- 0447 if ~any(Cm), break, end %-Break if empty mask 0448 Y(i,Cm) = spm_get_data(VY(i),xyz(:,Cm)); 0449 0450 Cm(Cm) = Y(i,Cm) > xM.TH(i); %-Threshold (& NaN) mask 0451 if xM.I & ~YNaNrep & xM.TH(i) < 0 %-Use implicit mask 0452 Cm(Cm) = abs(Y(i,Cm)) > eps; 0453 end 0454 end 0455 0456 %-Mask out voxels where data is constant 0457 %--------------------------------------------------------------- 0458 Cm(Cm) = any(diff(Y(:,Cm),1)); 0459 Y = Y(:,Cm); %-Data within mask 0460 CrS = sum(Cm); %- Number of current voxels 0461 0462 %-Conditional estimates (per partition, per voxel) 0463 %--------------------------------------------------------------- 0464 beta = zeros(nBeta,CrS); 0465 Hp = zeros(1, CrS); 0466 AR = zeros(PPM.AR_P, CrS); 0467 Psd = zeros(nPsd, CrS); 0468 0469 if CrS 0470 0471 % - Assume single session 0472 vxyz = spm_vb_neighbors(xyz(:,Cm)'); 0473 0474 slice = slice_template; 0475 slice.verbose=1; 0476 slice.update_w=1; 0477 slice.update_lambda=1; 0478 slice.update_F=0; 0479 0480 slice=spm_vb_set_priors(slice,PPM.priors,vxyz); 0481 slice=spm_vb_glmar(Y,slice); 0482 0483 if PPM.AR_P > 0 0484 AR(1:PPM.AR_P,:)=slice.ap_mean; 0485 end 0486 0487 if PPM.compute_F 0488 slice.F=spm_vb_F(Y,slice); 0489 PPM.slice(z).F=slice.F; 0490 end 0491 0492 beta = slice.wk_mean; 0493 Hp(1,:) = sqrt(1./slice.mean_lambda'); 0494 Psd = slice.w_dev; 0495 0496 % Get slice-wise Taylor approximation to posterior correlation 0497 slice = spm_vb_taylor_R (Y,slice); 0498 PPM.slice(z).mean=slice.mean; 0499 PPM.slice(z).elapsed_seconds=slice.elapsed_seconds; 0500 0501 0502 0503 end % if CrS 0504 0505 0506 %-Append new inmask voxel locations and volumes 0507 %--------------------------------------------------------------- 0508 XYZ(:,S + [1:CrS]) = xyz(:,Cm); %-InMask XYZ voxel coords 0509 Q = [Q I(Cm)]; %-InMask XYZ voxel indices 0510 S = S + CrS; %-Volume analysed (voxels) 0511 0512 % Sum of noise variance hyperparameters 0513 SHp = SHp + sum(Hp(1,:)); 0514 0515 %-Save for current plane in memory as we go along 0516 %--------------------------------------------------------------- 0517 CrBl = [CrBl beta]; 0518 CrHp = [CrHp Hp]; 0519 CrPsd = [CrPsd Psd]; 0520 CrAR = [CrAR AR]; 0521 0522 end % (bch) 0523 end % (if analyse_slice) 0524 0525 %-Write Mask image 0526 %------------------------------------------------------------------- 0527 j = sparse(xdim,ydim); 0528 if length(Q), j(Q) = 1; end 0529 VM = spm_write_plane(VM, j, z); 0530 0531 %-Write conditional beta images 0532 %------------------------------------------------------------------- 0533 j = NaN*ones(xdim,ydim); 0534 for i = 1:nBeta 0535 if length(Q), j(Q) = CrBl(i,:); end 0536 Vbeta(i) = spm_write_plane(Vbeta(i),j,z); 0537 end 0538 0539 %-Write SD error images 0540 %------------------------------------------------------------------- 0541 j = NaN*ones(xdim,ydim); 0542 if length(Q), j(Q) = CrHp(1,:); end 0543 VHp = spm_write_plane(VHp,j,z); 0544 0545 0546 %-Write posterior standard-deviation of beta images 0547 %------------------------------------------------------------------- 0548 j = NaN*ones(xdim,ydim); 0549 for i = 1:nPsd 0550 if length(Q), j(Q) = CrPsd(i,:); end 0551 VPsd(i) = spm_write_plane(VPsd(i),j,z); 0552 end 0553 0554 %-Write AR images 0555 %------------------------------------------------------------------- 0556 j = NaN*ones(xdim,ydim); 0557 for i = 1:PPM.AR_P 0558 if length(Q), j(Q) = CrAR(i,:); end 0559 VAR(i) = spm_write_plane(VAR(i),j,z); 0560 end 0561 0562 %-Report progress 0563 %------------------------------------------------------------------- 0564 spm_progress_bar('Set',100*(z - 1)/zdim); 0565 0566 0567 end % (for z = 1:zdim) 0568 0569 fprintf('\n') %-# 0570 spm_progress_bar('Clear') 0571 spm('Pointer','Arrow'); 0572 %======================================================================= 0573 % - P O S T E S T I M A T I O N 0574 %======================================================================= 0575 0576 if S == 0, warning('No inmask voxels - empty analysis!'), end 0577 0578 0579 %-"close" written image files, updating scalefactor information 0580 %======================================================================= 0581 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...closing files') %-# 0582 Vbeta = spm_close_vol(Vbeta); 0583 VHp = spm_close_vol(VHp); 0584 VPsd = spm_close_vol(VPsd); 0585 VAR = spm_close_vol(VAR); 0586 VM = spm_close_vol(VM); 0587 0588 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done') %-# 0589 0590 %-Create 1st contrast for 'effects of interest' (all if not specified) 0591 %======================================================================= 0592 Fcname = 'effects of interest'; 0593 try 0594 iX0 = [xX.iB xX.iG]; 0595 catch 0596 iX0 = []; 0597 end 0598 0599 xX.xKXs = spm_sp('Set',spm_filter(xX.K,xX.X)); % ** Not Whitened ** 0600 xX.erdf = size(xX.X,1); % Just set to number of scans so, when 0601 % we assess the results, spm_getSPM is happy 0602 xX.W= eye(size(xX.X,1)); % Set whitening matrix to identity - 0603 % we must set it to keep spm_graph happy 0604 0605 xCon = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs); 0606 0607 %-Compute scaled design matrix for display purposes 0608 %----------------------------------------------------------------------- 0609 xX.nKX = spm_DesMtx('sca',xX.xKXs.X,xX.name); 0610 0611 %-Save remaining results files and analysis parameters 0612 %======================================================================= 0613 fprintf('%-40s: %30s','Saving results','...writing') %-# 0614 0615 0616 0617 %-place fields in SPM 0618 %----------------------------------------------------------------------- 0619 SPM.xVol.XYZ = XYZ(:,1:S); %-InMask XYZ coords (voxels) 0620 SPM.xVol.M = M; %-voxels -> mm 0621 SPM.xVol.iM = inv(M); %-mm -> voxels 0622 SPM.xVol.DIM = DIM; %-image dimensions 0623 SPM.xVol.S = S; %-Volume (voxels) 0624 SPM.xVol.R = 100; % Set R - number of RESELS - to arbitrary value 0625 % as, if R not set, SPM will think model has not 0626 % been estimated 0627 SPM.xVol.FWHM = 10; % Set to arbitrary value so spm_getSPM is happy 0628 0629 SPM.VCbeta = Vbeta; % Filenames - parameters 0630 SPM.VHp = VHp; % Filenames - hyperparameters 0631 SPM.VPsd = VPsd; % Filenames - hyperparameters 0632 SPM.VAR = VAR; % Filenames - hyperparameters 0633 SPM.VM = VM; %-Filehandle - Mask 0634 0635 %PPM.Cb = 1; % Default threshold for effect size (1 per cent) 0636 0637 PPM.Gamma = 1; % Default threshold for effect size (1 per cent) 0638 SPM.PPM = PPM; % PPM structure 0639 0640 SPM.xX = xX; %-design structure 0641 SPM.xM = xM; %-mask structure 0642 0643 SPM.xCon = xCon; %-contrast structure 0644 0645 save SPM SPM 0646 0647 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done') %-# 0648 0649 0650 %======================================================================= 0651 %- E N D: Cleanup GUI 0652 %======================================================================= 0653 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') 0654 fprintf('%-40s: %30s\n','Completed',spm('time')) %-# 0655 fprintf('...use the results section for assessment\n\n') %-# 0656 0657