123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482 |
- function [SPM] = spm_spm_Bayes(SPM)
- % 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 empirical 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 empirical Bayesian prior variance estimators over
- % voxels.
- %
- % Each separable partition (i.e. session) is assigned its own
- % hyperparameter but within session covariance components are lumped
- % together, using their relative expectations over voxels. This makes
- % things much more computationally efficient and avoids inefficient
- % voxel-specific multiple hyperparameter estimates.
- %
- % spm_spm_Bayes adds the following fields to SPM:
- %
- % ----------------
- %
- %
- % SPM.PPM.l = session-specific hyperparameter means
- % SPM.PPM.Cb = empirical prior parameter covariances
- % SPM.PPM.C = conditional covariances of parameters
- % SPM.PPM.dC{i} = dC/dl;
- % SPM.PPM.ddC{i} = ddC/dldl
- %
- % The derivatives are used to compute the conditional variance of various
- % contrasts in spm_getSPM, using a Taylor expansion about the hyperparameter
- % means.
- %
- %
- % ----------------
- %
- % SPM.VCbeta - Handles of conditional parameter estimates
- % SPM.VHp - Handles of hyperparameter estimates
- %
- % ----------------
- %
- % Cbeta_????.<ext> - conditional parameter images
- % These are 32-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.<ext>) are given
- % value NaN.
- %
- % ----------------
- %
- % CHp_????.<ext> - error covariance hyperparameter images
- % This is a 64-bit (double) image of the ReML error variance estimate.
- % for each separable partition (Session). Voxels outside the analysis
- % mask are given value NaN.
- %__________________________________________________________________________
- %
- % For single subject fMRI analysis there is an alternative function
- % using voxel-wise GLM-AR models that are spatially regularised
- % using the VB framework. This is implemented using spm_spm_vb.m.
- %__________________________________________________________________________
- % Copyright (C) 2002-2013 Wellcome Trust Centre for Neuroimaging
- % Karl Friston
- % $Id: spm_spm_Bayes.m 5354 2013-03-25 18:50:54Z guillaume $
- %-Say hello
- %--------------------------------------------------------------------------
- Finter = spm('FigName','Stats: Bayesian estimation...');
- %-Select SPM.mat & change directory
- %--------------------------------------------------------------------------
- if ~nargin
- [Pf, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
- if ~sts, return; end
- swd = spm_file(Pf,'fpath');
- load(fullfile(swd,'SPM.mat'))
- cd(swd)
- end
- try
- M = SPM.xVol.M;
- DIM = SPM.xVol.DIM;
- xdim = DIM(1); ydim = DIM(2); zdim = DIM(3);
- XYZ = SPM.xVol.XYZ;
- catch
- helpdlg({ 'Please do a ML estimation first.',...
- 'This identifies the voxels to analyse.'});
- spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
- return
- end
- %==========================================================================
- % - A N A L Y S I S P R E L I M I N A R I E S
- %==========================================================================
- %-Initialise output images
- %==========================================================================
- fprintf('%-40s: %30s','Output images','...initialising') %-#
- %-Initialise conditional estimate image files
- %--------------------------------------------------------------------------
- xX = SPM.xX;
- [nScan,nBeta] = size(xX.X);
- Vbeta(1:nBeta) = deal(struct(...
- 'fname', [],...
- 'dim', DIM',...
- 'dt', [spm_type('float32'), spm_platform('bigend')],...
- 'mat', M,...
- 'pinfo', [1 0 0]',...
- 'descrip', ''));
- for i = 1:nBeta
- Vbeta(i).fname = [sprintf('Cbeta_%04d',i) spm_file_ext];
- Vbeta(i).descrip = sprintf('Cond. beta (%04d) - %s',i,xX.name{i});
- spm_unlink(Vbeta(i).fname)
- end
- Vbeta = spm_create_vol(Vbeta);
- %-Initialise ReML hyperparameter image files
- %--------------------------------------------------------------------------
- try
- nHp = length(SPM.nscan);
- catch
- nHp = nScan;
- SPM.nscan = nScan;
- end
- VHp(1:nHp) = deal(struct(...
- 'fname', [],...
- 'dim', DIM',...
- 'dt', [spm_type('float64'), spm_platform('bigend')],...
- 'mat', M,...
- 'pinfo', [1 0 0]',...
- 'descrip', ''));
- for i = 1:nHp
- VHp(i).fname = [sprintf('Hp_%04d',i) spm_file_ext];
- VHp(i).descrip = sprintf('Hyperparameter (%04d)',i);
- spm_unlink(VHp(i).fname)
- end
- VHp = spm_create_vol(VHp);
- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...initialised') %-#
- %==========================================================================
- % - A V E R A G E S A M P L E C O V A R I A N C E M A T R I X
- %==========================================================================
- fprintf('%s%30s',repmat(sprintf('\b'),1,30),'...estimating CY') %-#
- CY = 0; % <(Y - <Y>) * (Y - <Y>)'>
- EY = 0; % <Y> for ReML
- nScan = size(xX.X,1);
- xVi = SPM.xVi;
- %-Compute Hsqr and F-threshold under i.i.d.
- %--------------------------------------------------------------------------
- xX.xKXs = spm_sp('Set',spm_filter(xX.K,xX.W*xX.X));
- xX.xKXs.X = full(xX.xKXs.X);
- xX.pKX = spm_sp('x-',xX.xKXs);
- if isfield(xVi,'Fcontrast')
- Fcname = 'User-specified contrast';
- xCon = spm_FcUtil('Set',Fcname,'F','c',xVi.Fcontrast,xX.xKXs);
- else
- Fcname = 'effects of interest';
- iX0 = [xX.iB xX.iG];
- xCon = spm_FcUtil('Set',Fcname,'F','iX0',iX0,xX.xKXs);
- end
- if ~isempty(xCon(1).c)
- X1o = spm_FcUtil('X1o', xCon(1),xX.xKXs);
- Hsqr = spm_FcUtil('Hsqr',xCon(1),xX.xKXs);
- trMV = spm_SpUtil('trMV',X1o);
- else
- % Force all voxels to enter non-sphericity
- trMV = 1;
- Hsqr = Inf;
- end
- trRV = spm_SpUtil('trRV',xX.xKXs);
- %-Threshold for voxels entering non-sphericity estimates
- %--------------------------------------------------------------------------
- try
- modality = lower(spm_get_defaults('modality'));
- UFp = spm_get_defaults(['stats.' modality '.ufp']);
- catch
- UFp = 0.001;
- end
- xVi.UFp = UFp;
- UF = spm_invFcdf(1 - UFp,[trMV,trRV]);
- %-Split data into chunks
- %--------------------------------------------------------------------------
- VY = SPM.xY.VY;
- mask = logical(spm_read_vols(SPM.VM));
- chunksize = floor(spm_get_defaults('stats.maxmem') / 8 / nScan);
- nbchunks = ceil(prod(DIM) / chunksize);
- chunks = min(cumsum([1 repmat(chunksize,1,nbchunks)]),prod(DIM)+1);
- for i=1:nbchunks
- chunk = chunks(i):chunks(i+1)-1;
-
- %-Get data & construct analysis mask
- %----------------------------------------------------------------------
- Y = zeros(nScan,numel(chunk));
- cmask = mask(chunk);
- for j=1:nScan
- if ~any(cmask), break, end %-Break if empty mask
- Y(j,cmask) = spm_data_read(VY(j),chunk(cmask));%-Read chunk of data
- end
- mask(chunk) = cmask;
- if ~any(cmask), continue, end
- Y = Y(:,cmask); %-Data within mask
- %-Remove filter confounds
- %----------------------------------------------------------------------
- KWY = spm_filter(xX.K,xX.W*Y);
-
- %-Ordinary Least Squares estimation
- %----------------------------------------------------------------------
- beta = xX.pKX*KWY; %-Parameter estimates
- if any(cmask)
- res = spm_sp('r',xX.xKXs,KWY); %-Residuals
- else
- res = zeros(nScan,0);
- end
- ResSS = sum(res.^2); %-Residual SSQ
- clear res
-
- %-F-threshold & accumulate spatially whitened Y*Y'
- %----------------------------------------------------------------------
- j = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV;
- if nnz(j)
- Y = Y(:,j);
- CY = CY + Y*Y';
- EY = EY + sum(Y,2);
- end
-
- end
- %-average sample covariance and mean of Y (over voxels)
- %--------------------------------------------------------------------------
- S = nnz(mask);
- CY = CY/S;
- EY = EY/S;
- CY = CY - EY*EY';
- SPM.xVi.CY = CY;
- clear CY EY
- %==========================================================================
- % - E M P I R I C A L B A Y E S F O R P R I O R V A R I A N C E
- %==========================================================================
- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...estimating priors') %-#
- % get row u{i} and column v{i}/v0{i} indices for separable designs
- %--------------------------------------------------------------------------
- s = nHp;
- if isfield(SPM,'Sess')
- for i = 1:s
- u{i} = SPM.Sess(i).row;
- v{i} = SPM.Sess(i).col;
- v0{i} = xX.iB(i);
- end
- else
- u{1} = [1:nScan];
- v{1} = [xX.iH xX.iC];
- v0{1} = [xX.iB xX.iG];
- end
- % cycle over separarable partitions
- %--------------------------------------------------------------------------
- for i = 1:s
- % Get design X and confounds X0
- %----------------------------------------------------------------------
- fprintf('%-30s\n',sprintf(' ReML Session %i',i)); %-#
- X = xX.X(u{i}, v{i});
- X0 = xX.X(u{i},v0{i});
- [m,n] = size(X);
- % add confound in 'filter'
- %----------------------------------------------------------------------
- if isstruct(xX.K)
- X0 = full([X0 xX.K(i).X0]);
- end
- % orthogonalize X w.r.t. X0
- %----------------------------------------------------------------------
- X = X - X0*(pinv(X0)*X);
- % covariance components induced by parameter variations {Q}
- %----------------------------------------------------------------------
- for j = 1:n
- Q{j} = X*sparse(j,j,1,n,n)*X';
- end
- % covariance components induced by error non-sphericity {V}
- %----------------------------------------------------------------------
- Q{n + 1} = SPM.xVi.V(u{i},u{i});
- % ReML covariance component estimation
- %----------------------------------------------------------------------
- [C,h] = spm_reml(SPM.xVi.CY,X0,Q);
- % check for negative variance components
- %----------------------------------------------------------------------
- h = abs(h);
- % 2-level model for this partition using prior variances sP(i)
- % treat confounds as fixed (i.e. infinite prior variance)
- %----------------------------------------------------------------------
- n0 = size(X0,2);
- Cb = blkdiag(diag(h(1:n)),speye(n0,n0)*1e8);
- P{1}.X = [X X0];
- P{1}.C = {SPM.xVi.V};
- P{2}.X = sparse(size(P{1}.X,2),1);
- P{2}.C = Cb;
- sP(i).P = P;
- sP(i).u = u{:};
- sP(i).v = [v{:} v0{:}];
- end
- %==========================================================================
- % - F I T M O D E L & W R I T E P A R A M E T E R I M A G E S
- %==========================================================================
- %-Cycle to avoid memory problems (plane by plane)
- %==========================================================================
- spm_progress_bar('Init',100,'Bayesian estimation','');
- spm('Pointer','Watch')
- %-maxMem is the maximum amount of data processed at a time (bytes)
- %--------------------------------------------------------------------------
- MAXMEM = spm_get_defaults('stats.maxmem');
- blksz = ceil(MAXMEM/8/nScan);
- SHp = 0; % sum of hyperparameters
- for z = 1:zdim
- % current plane-specific parameters
- %----------------------------------------------------------------------
- U = find(XYZ(3,:) == z);
- nbch = ceil(length(U)/blksz);
- CrBl = zeros(nBeta,length(U)); %-conditional parameter estimates
- CrHp = zeros(nHp, length(U)); %-ReML hyperparameter estimates
- for bch = 1:nbch %-loop over bunches of lines (planks)
-
- %-construct list of voxels in this block
- %------------------------------------------------------------------
- I = [1:blksz] + (bch - 1)*blksz;
- I = I(I <= length(U));
- xyz = XYZ(:,U(I));
- nVox = size(xyz,2);
-
- %-Get response variable
- %------------------------------------------------------------------
- Y = spm_get_data(SPM.xY.VY,xyz);
-
- %-Conditional estimates (per partition, per voxel)
- %------------------------------------------------------------------
- beta = zeros(nBeta,nVox);
- Hp = zeros(nHp, nVox);
- for j = 1:s
- P = sP(j).P;
- u = sP(j).u;
- v = sP(j).v;
- for i = 1:nVox
- C = spm_PEB(Y(u,i),P);
- beta(v,i) = C{2}.E(1:length(v));
- Hp(j,i) = C{1}.h;
- end
- end
-
- %-Save for current plane in memory as we go along
- %------------------------------------------------------------------
- CrBl(:,I) = beta;
- CrHp(:,I) = Hp;
- SHp = SHp + sum(Hp,2);
-
- end % (bch)
- %-write out plane data to image files
- %======================================================================
- %-Write conditional beta images
- %----------------------------------------------------------------------
- for i = 1:nBeta
- tmp = sparse(XYZ(1,U),XYZ(2,U),CrBl(i,:),xdim,ydim);
- tmp(~tmp) = NaN;
- Vbeta(i) = spm_write_plane(Vbeta(i),tmp,z);
- end
- %-Write hyperparameter images
- %----------------------------------------------------------------------
- for i = 1:nHp
- tmp = sparse(XYZ(1,U),XYZ(2,U),CrHp(i,:),xdim,ydim);
- tmp(~tmp) = NaN;
- VHp(i) = spm_write_plane(VHp(i),tmp,z);
- end
- %-Report progress
- %----------------------------------------------------------------------
- spm_progress_bar('Set',100*(z - 1)/zdim);
- end % (for z = 1:zdim)
- fprintf('\n') %-#
- spm_progress_bar('Clear')
- %==========================================================================
- % - P O S T E S T I M A T I O N
- %==========================================================================
- % Taylor expansion for conditional covariance
- %--------------------------------------------------------------------------
- fprintf('%-40s: %30s\n','Non-sphericity','...REML estimation') %-#
- % expansion point (mean hyperparameters)
- %--------------------------------------------------------------------------
- l = SHp/SPM.xVol.S;
- % change in conditional coavriance w.r.t. hyperparameters
- %--------------------------------------------------------------------------
- n = size(xX.X,2);
- PPM.l = l;
- for i = 1:s
- PPM.dC{i} = sparse(n,n);
- PPM.ddC{i} = sparse(n,n);
- end
- for i = 1:s
- P = sP(i).P;
- u = sP(i).u;
- v = sP(i).v;
- % derivatives of conditional covariance w.r.t. hyperparameters
- %----------------------------------------------------------------------
- d = P{1}.X'*inv(P{1}.C{1})*P{1}.X;
- Cby = inv(d/l(i) + inv(P{2}.C));
- d = d*Cby;
- dC = Cby*d/(l(i)^2);
- ddC = 2*(dC/(l(i)^2) - Cby/(l(i)^3))*d;
- % place in output structure
- %----------------------------------------------------------------------
- j = 1:length(v);
- PPM.Cb(v,v) = P{2}.C(j,j);
- PPM.Cby(v,v) = Cby(j,j);
- PPM.dC{i}(v,v) = dC(j,j);
- PPM.ddC{i}(v,v) = ddC(j,j);
-
- end
- %-Save remaining results files and analysis parameters
- %==========================================================================
- fprintf('%-40s: %30s','Saving results','...writing') %-#
- %-Save analysis parameters in SPM.mat file
- %--------------------------------------------------------------------------
- SPM.VCbeta = Vbeta; % Filenames - parameters
- SPM.VHp = VHp; % Filenames - hyperparameters
- SPM.PPM = PPM; % PPM structure
- save('SPM.mat', 'SPM', spm_get_defaults('mat.format'));
- fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done') %-#
- %==========================================================================
- %- E N D: Cleanup GUI
- %==========================================================================
- spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
- fprintf('%-40s: %30s\n','Completed',spm('time')) %-#
|