Conditional parameter estimation of a General Linear Model
FORMAT [SPM] = spm_spm_Bayes(SPM)
_______________________________________________________________________
spm_spm_Bayes returns to voxels identified by spm_spm (ML parameter
estimation) to get conditional parameter estimates and ReML hyper-
parameter estimates. These estimates use prior covariances, on the
parameters, from emprical Bayes. These PEB prior variances come from
the hierarchical model that obtains by considering voxels as providing
a second level. Put simply, the variance in parameters, over voxels,
is used as a prior variance from the point of view of any one voxel.
The error covariance hyperparameters are re-estimated in the light of
these priors. The approach adopted is essentially a fully Bayesian
analysis at each voxel, using emprical Bayesian prior variance
estimators over voxels.
Each separable partition (i.e. session) is assigned its own
hyperparameter but within session covariance components are lumped
together, using their relative expectations over voxels. This makes
things much more computationally efficient and avoids inefficient
voxel-specific multiple hyperparameter estimates.
spm_spm_Bayes adds the following feilds to SPM:
----------------
SPM.PPM.l = session-specific hyperparameter means
SPM.PPM.Cb = empirical prior parameter covariances
SPM.PPM.C = conditional covariances of parameters
SPM.PPM.dC{i} = dC/dl;
SPM.PPM.ddC{i} = ddC/dldl
The derivatives are used to compute the conditional variance of various
contrasts in spm_getSPM, using a Taylor expansion about the hyperparameter
means.
----------------
SPM.VCbeta - Handles of conditional parameter estimates
SPM.VHp - Handles of hyperparameter estimates
----------------
Cbeta_????.{img,hdr} - conditional parameter images
These are 16-bit (float) images of the conditional estimates. The image
files are numbered according to the corresponding column of the
design matrix. Voxels outside the analysis mask (mask.img) are given
value NaN.
----------------
CHp_????.{img,hdr} - error covariance hyperparamter images
This is a 32-bit (double) image of the ReML error variance estimate.
for each separable partition (Session). Voxels outside the analysis
mask are given value NaN.
_______________________________________________________________________
@(#)spm_spm_Bayes.m 2.7 Karl Friston 03/03/120001 function [SPM] = spm_spm_Bayes(SPM) 0002 % Conditional parameter estimation of a General Linear Model 0003 % FORMAT [SPM] = spm_spm_Bayes(SPM) 0004 %_______________________________________________________________________ 0005 % spm_spm_Bayes returns to voxels identified by spm_spm (ML parameter 0006 % estimation) to get conditional parameter estimates and ReML hyper- 0007 % parameter estimates. These estimates use prior covariances, on the 0008 % parameters, from emprical Bayes. These PEB prior variances come from 0009 % the hierarchical model that obtains by considering voxels as providing 0010 % a second level. Put simply, the variance in parameters, over voxels, 0011 % is used as a prior variance from the point of view of any one voxel. 0012 % The error covariance hyperparameters are re-estimated in the light of 0013 % these priors. The approach adopted is essentially a fully Bayesian 0014 % analysis at each voxel, using emprical Bayesian prior variance 0015 % estimators over voxels. 0016 % 0017 % Each separable partition (i.e. session) is assigned its own 0018 % hyperparameter but within session covariance components are lumped 0019 % together, using their relative expectations over voxels. This makes 0020 % things much more computationally efficient and avoids inefficient 0021 % voxel-specific multiple hyperparameter estimates. 0022 % 0023 % spm_spm_Bayes adds the following feilds to SPM: 0024 % 0025 % ---------------- 0026 % 0027 % 0028 % SPM.PPM.l = session-specific hyperparameter means 0029 % SPM.PPM.Cb = empirical prior parameter covariances 0030 % SPM.PPM.C = conditional covariances of parameters 0031 % SPM.PPM.dC{i} = dC/dl; 0032 % SPM.PPM.ddC{i} = ddC/dldl 0033 % 0034 % The derivatives are used to compute the conditional variance of various 0035 % contrasts in spm_getSPM, using a Taylor expansion about the hyperparameter 0036 % means. 0037 % 0038 % 0039 % ---------------- 0040 % 0041 % SPM.VCbeta - Handles of conditional parameter estimates 0042 % SPM.VHp - Handles of hyperparameter estimates 0043 % 0044 % ---------------- 0045 % 0046 % Cbeta_????.{img,hdr} - conditional parameter images 0047 % These are 16-bit (float) images of the conditional estimates. The image 0048 % files are numbered according to the corresponding column of the 0049 % design matrix. Voxels outside the analysis mask (mask.img) are given 0050 % value NaN. 0051 % 0052 % ---------------- 0053 % 0054 % CHp_????.{img,hdr} - error covariance hyperparamter images 0055 % This is a 32-bit (double) image of the ReML error variance estimate. 0056 % for each separable partition (Session). Voxels outside the analysis 0057 % mask are given value NaN. 0058 % 0059 %_______________________________________________________________________ 0060 % @(#)spm_spm_Bayes.m 2.7 Karl Friston 03/03/12 0061 0062 0063 %-Say hello 0064 %----------------------------------------------------------------------- 0065 Finter = spm('FigName','Stats: Bayesian estimation...'); 0066 0067 %-Select SPM.mat & change directory 0068 %----------------------------------------------------------------------- 0069 if ~nargin 0070 swd = spm_str_manip(spm_get(1,'SPM.mat','Select SPM.mat'),'H'); 0071 load(fullfile(swd,'SPM.mat')) 0072 cd(swd) 0073 end 0074 0075 % Single subject fMRI ? 0076 str = 'Single Subject fMRI ?'; 0077 ss_fMRI = spm_input(str,1,'b',{'yes','no'},[1 0]); 0078 if ss_fMRI 0079 % Use VB-GLM-AR algorithm 0080 spm_spm_vb(SPM); 0081 return 0082 end 0083 0084 try 0085 M = SPM.xVol.M; 0086 DIM = SPM.xVol.DIM; 0087 xdim = DIM(1); ydim = DIM(2); zdim = DIM(3); 0088 XYZ = SPM.xVol.XYZ; 0089 catch 0090 helpdlg({ 'Please do a ML estimation first',... 0091 'This identifies the voxels to analyse'}); 0092 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') 0093 return 0094 end 0095 0096 0097 %======================================================================= 0098 % - A N A L Y S I S P R E L I M I N A R I E S 0099 %======================================================================= 0100 0101 %-Initialise output images 0102 %======================================================================= 0103 fprintf('%-40s: %30s','Output images','...initialising') %-# 0104 0105 %-Intialise oonditional estimate image files 0106 %----------------------------------------------------------------------- 0107 xX = SPM.xX; 0108 [nScan nBeta] = size(xX.X); 0109 Vbeta(1:nBeta) = deal(struct(... 0110 'fname', [],... 0111 'dim', [DIM',spm_type('float')],... 0112 'mat', M,... 0113 'pinfo', [1 0 0]',... 0114 'descrip', '')); 0115 for i = 1:nBeta 0116 Vbeta(i).fname = sprintf('Cbeta_%04d.img',i); 0117 Vbeta(i).descrip = sprintf('Cond. beta (%04d) - %s',i,xX.name{i}); 0118 spm_unlink(Vbeta(i).fname) 0119 end 0120 Vbeta = spm_create_vol(Vbeta,'noopen'); 0121 0122 %-Intialise ReML hyperparameter image files 0123 %----------------------------------------------------------------------- 0124 try 0125 nHp = length(SPM.nscan); 0126 catch 0127 nHp = nScan; 0128 SPM.nscan = nScan; 0129 end 0130 0131 VHp(1:nHp) = deal(struct(... 0132 'fname', [],... 0133 'dim', [DIM',spm_type('double')],... 0134 'mat', M,... 0135 'pinfo', [1 0 0]',... 0136 'descrip', '')); 0137 for i = 1:nHp 0138 VHp(i).fname = sprintf('Hp_%04d.img',i); 0139 VHp(i).descrip = sprintf('Hyperparameter (%04d)',i); 0140 spm_unlink(VHp(i).fname) 0141 end 0142 VHp = spm_create_vol(VHp,'noopen'); 0143 0144 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...initialised') %-# 0145 0146 0147 %======================================================================= 0148 % - E M P I R I C A L B A Y E S F O R P R I O R V A R I A N C E 0149 %======================================================================= 0150 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...estimatng priors') %-# 0151 0152 % get row u{i} and column v{i}/v0{i} indices for separable designs 0153 %---------------------------------------------------------------------- 0154 s = nHp; 0155 if isfield(SPM,'Sess') 0156 for i = 1:s 0157 u{i} = SPM.Sess(i).row; 0158 v{i} = SPM.Sess(i).col; 0159 v0{i} = xX.iB(i); 0160 end 0161 else 0162 u{1} = [1:nScan]; 0163 v{1} = [xX.iH xX.iC]; 0164 v0{1} = [xX.iB xX.iG]; 0165 end 0166 0167 % cycle over separarable partitions 0168 %----------------------------------------------------------------------- 0169 for i = 1:s 0170 0171 % Get design X and confounds X0 0172 %--------------------------------------------------------------- 0173 fprintf('%-30s- %i\n',' ReML Session',i); 0174 X = xX.X(u{i}, v{i}); 0175 X0 = xX.X(u{i},v0{i}); 0176 [m n] = size(X); 0177 0178 % add confound in 'filter' 0179 %--------------------------------------------------------------- 0180 if isstruct(xX.K) 0181 X0 = full([X0 xX.K(i).X0]); 0182 end 0183 0184 % orthogonalize X w.r.t. X0 0185 %--------------------------------------------------------------- 0186 X = X - X0*(pinv(X0)*X); 0187 0188 % covariance components induced by parameter variations {Q} 0189 %--------------------------------------------------------------- 0190 for j = 1:n 0191 Q{j} = X*sparse(j,j,1,n,n)*X'; 0192 end 0193 0194 % covariance components induced by error non-sphericity {V} 0195 %--------------------------------------------------------------- 0196 Q{n + 1} = SPM.xVi.V(u{i},u{i}); 0197 0198 % ReML covariance component estimation 0199 %--------------------------------------------------------------- 0200 [C h W] = spm_reml(SPM.xVi.CY,X0,Q); 0201 0202 % check for negative variance components 0203 %--------------------------------------------------------------- 0204 h = abs(h); 0205 0206 % 2-level model for this partition using prior variances sP(i) 0207 % treat confounds as fixed (i.e. infinite prior variance) 0208 %--------------------------------------------------------------- 0209 n0 = size(X0,2); 0210 Cb = blkdiag(diag(h(1:n)),speye(n0,n0)*1e8); 0211 P{1}.X = [X X0]; 0212 P{1}.C = {SPM.xVi.V}; 0213 P{2}.X = sparse(size(P{1}.X,2),1); 0214 P{2}.C = Cb; 0215 0216 sP(i).P = P; 0217 sP(i).u = u{:}; 0218 sP(i).v = [v{:} v0{:}]; 0219 end 0220 0221 0222 %======================================================================= 0223 % - 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 0224 %======================================================================= 0225 0226 %-Cycle to avoid memory problems (plane by plane) 0227 %======================================================================= 0228 spm_progress_bar('Init',100,'Bayesian estimation',''); 0229 helpdlg('This could take some time'); 0230 spm('Pointer','Watch') 0231 0232 %-maxMem is the maximum amount of data processed at a time (bytes) 0233 %----------------------------------------------------------------------- 0234 global defaults 0235 MAXMEM = defaults.stats.maxmem; 0236 blksz = ceil(MAXMEM/8/nScan); 0237 SHp = 0; % sum of hyperparameters 0238 for z = 1:zdim 0239 0240 % current plane-specific parameters 0241 %------------------------------------------------------------------- 0242 U = find(XYZ(3,:) == z); 0243 nbch = ceil(length(U)/blksz); 0244 CrBl = zeros(nBeta,length(U)); %-conditional parameter estimates 0245 CrHp = zeros(nHp, length(U)); %-ReML hyperparameter estimates 0246 for bch = 1:nbch %-loop over bunches of lines (planks) 0247 0248 %-construct list of voxels in this block 0249 %--------------------------------------------------------------- 0250 I = [1:blksz] + (bch - 1)*blksz; 0251 I = I(I <= length(U)); 0252 xyz = XYZ(:,U(I)); 0253 nVox = size(xyz,2); 0254 0255 %-Get response variable 0256 %--------------------------------------------------------------- 0257 Y = spm_get_data(SPM.xY.VY,xyz); 0258 0259 %-Conditional estimates (per partition, per voxel) 0260 %--------------------------------------------------------------- 0261 beta = zeros(nBeta,nVox); 0262 Hp = zeros(nHp, nVox); 0263 for j = 1:s 0264 P = sP(j).P; 0265 u = sP(j).u; 0266 v = sP(j).v; 0267 for i = 1:nVox 0268 [C P] = spm_PEB(Y(u,i),P); 0269 beta(v,i) = C{2}.E(1:length(v)); 0270 Hp(j,i) = C{1}.h; 0271 end 0272 end 0273 0274 %-Save for current plane in memory as we go along 0275 %--------------------------------------------------------------- 0276 CrBl(:,I) = beta; 0277 CrHp(:,I) = Hp; 0278 SHp = SHp + sum(Hp,2); 0279 0280 end % (bch) 0281 0282 0283 %-write out plane data to image files 0284 %=================================================================== 0285 0286 %-Write conditional beta images 0287 %------------------------------------------------------------------- 0288 for i = 1:nBeta 0289 tmp = sparse(XYZ(1,U),XYZ(2,U),CrBl(i,:),xdim,ydim); 0290 tmp = tmp + NaN*(~tmp); 0291 Vbeta(i) = spm_write_plane(Vbeta(i),tmp,z); 0292 end 0293 0294 %-Write hyperparameter images 0295 %------------------------------------------------------------------- 0296 for i = 1:nHp 0297 tmp = sparse(XYZ(1,U),XYZ(2,U),CrHp(i,:),xdim,ydim); 0298 tmp = tmp + NaN*(~tmp); 0299 VHp(i) = spm_write_plane(VHp(i),tmp,z); 0300 end 0301 0302 0303 %-Report progress 0304 %------------------------------------------------------------------- 0305 spm_progress_bar('Set',100*(z - 1)/zdim); 0306 0307 0308 end % (for z = 1:zdim) 0309 fprintf('\n') %-# 0310 spm_progress_bar('Clear') 0311 0312 %======================================================================= 0313 % - P O S T E S T I M A T I O N 0314 %======================================================================= 0315 0316 % Taylor expansion for conditional covariance 0317 %----------------------------------------------------------------------- 0318 fprintf('%-40s: %30s\n','Non-sphericity','...REML estimation') %-# 0319 0320 % expansion point (mean hyperparameters) 0321 %----------------------------------------------------------------------- 0322 l = SHp/SPM.xVol.S; 0323 0324 % change in conditional coavriance w.r.t. hyperparameters 0325 %----------------------------------------------------------------------- 0326 n = size(xX.X,2); 0327 PPM.l = l; 0328 for i = 1:s 0329 PPM.dC{i} = sparse(n,n); 0330 PPM.ddC{i} = sparse(n,n); 0331 end 0332 for i = 1:s 0333 0334 P = sP(i).P; 0335 u = sP(i).u; 0336 v = sP(i).v; 0337 0338 % derivatives of conditional covariance w.r.t. hyperparameters 0339 %--------------------------------------------------------------- 0340 d = P{1}.X'*inv(P{1}.C{1})*P{1}.X; 0341 Cby = inv(d/l(i) + inv(P{2}.C)); 0342 d = d*Cby; 0343 dC = Cby*d/(l(i)^2); 0344 ddC = 2*(dC/(l(i)^2) - Cby/(l(i)^3))*d; 0345 0346 % place in output structure 0347 %--------------------------------------------------------------- 0348 j = 1:length(v); 0349 PPM.Cb(v,v) = P{2}.C(j,j); 0350 PPM.Cby(v,v) = Cby(j,j); 0351 PPM.dC{i}(v,v) = dC(j,j); 0352 PPM.ddC{i}(v,v) = ddC(j,j); 0353 end 0354 0355 0356 %-"close" written image files, updating scalefactor information 0357 %======================================================================= 0358 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...closing files') %-# 0359 Vbeta = spm_close_vol(Vbeta); 0360 VHp = spm_close_vol(VHp); 0361 0362 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done') %-# 0363 0364 0365 %-Save remaining results files and analysis parameters 0366 %======================================================================= 0367 fprintf('%-40s: %30s','Saving results','...writing') %-# 0368 0369 %-Save analysis parameters in SPM.mat file 0370 %----------------------------------------------------------------------- 0371 SPM.VCbeta = Vbeta; % Filenames - parameters 0372 SPM.VHp = VHp; % Filenames - hyperparameters 0373 SPM.PPM = PPM; % PPM structure 0374 0375 save SPM SPM 0376 0377 fprintf('%s%30s\n',sprintf('\b')*ones(1,30),'...done') %-# 0378 0379 0380 %======================================================================= 0381 %- E N D: Cleanup GUI 0382 %======================================================================= 0383 spm('FigName','Stats: done',Finter); spm('Pointer','Arrow') 0384 fprintf('%-40s: %30s\n','Completed',spm('time')) %-# 0385 fprintf('...use the results section for assessment\n\n') %-#