User:Peterz

From Wikibooks, open books for an open world
Jump to: navigation, search

Draft Notes ahead of contributions[edit]

Getting Bayesian: Priors and Posteriors[edit]

All parameters in DCM, such as connections between regions, have prior values (beliefs prior to seeing the data) and posterior values (estimated parameter values after seeing the data).

Priors[edit]

If you take a look at a DCM, you will find the priors in DCM.M.pE (prior mean of each parameter) and DCM.M.pC (prior covariance of each parameter). Whereas DCM.M.pE is simply laid out with matrices for A, B, C etc, DCM.M.pC is an nxn matrix, where n is the number of parameters. To make DCM.M.pC intelligible, you can convert this to the format with separate matrices for A,B,C etc. The code is as follows:

Cp_matrix = spm_unvec(diag(DCM.M.pC), DCM.Ep);

If you'd like to know the priors on each connection in the A-matrix, you can then simply do:

full(Cp_matrix.A)

The priors for DCM are set by the script spm_dcm_fmri_priors.m . If you wish to change any of them, however, don't alter this script. Rather, set your new priors in DCM.options:

% These two override spm_dcm_fmri_priors
DCM.options.pE = my_pE; % Prior means
DCM.options.pC = my_pC; % Prior covariance
 
% These are used only if the above two are not set:
DCM.options.precision = x; % exp(x) is the prior covariance on DCM.a entries
DCM.options.decay = y; % y is multiplied by the prior means

Post-hoc DCM[edit]

A walk-through of the code[edit]

1. Bringing the number of free parameters to less than 16[edit]

Each iteration of the greedy search begins by listing the indices of A,B and C-matrix parameters with non-zero prior variance. The resulting list of parameters, k, are those which post-hoc DCM will switch on and off while searching the model space:

k = spm_fieldindices(DCM.Ep,'A','B','C');
k = k(C(k));
nparam = length(k);

If there are more than 16 parameters, those with the least evidence are removed. To prepare for this, priors and posteriors of the parameters are collected into vectors - prior variables begin with 'p', posteriors begin with 'q'. Any parameters with zero prior covariance are removed from all these variables, performed by conducting an SVD on the prior covariance matrix:

U     = spm_svd(pC);
qE    = U'*spm_vec(qE); % Posterior means stripped of those with zero prior covariance

The code then loops over each parameter in k (those with non-zero prior variance), and creates an updated prior variance vector with one parameter switched off, by having its variance set to zero. This is then used to generate a reduced prior covariance matrix (using the matrix U from the SVD above, rather than pc):

r      = C; r(k(i)) = 0;
R      = U(r,:)'*U(r,:); % The updated prior covariance is later given by R*pC*R

Armed with this reduced prior covariance matrix without the ith parameter, the full model can now be compared to the reduced model:

Z(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);

After repeating this for each parameter, only the 8 parameters in k which when removed give the strongest evidence are retained.

2. Searching all possible reduced models[edit]

Safe in the knowledge that the number of free parameters (listed in variable k) is not excessive, a list of all possible prior parameter vectors is generated, each with a different parameter switched off:

K = spm_perm_mtx(length(k));

Models with each of these prior covariances are then evaluated, just as described above:

for i = 1:length(K)
  r      = C; r(k(K(i,:))) = 0;
  R      = U(r,:)'*U(r,:);
  G(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);            
  ...
end

The results of the search across models are then collated (variable S). If models from multiple subjects are given, their log evidence (free energy) is pooled by adding the free energy of each. The evidence is then made relative to that of the full model (the last entry in S):

S     = sum(G,2);
S     = S - S(end);

And converted to a posterior probability for each model:

p     = exp(S - max(S));
p     = p/sum(p);

If we started with more than 16 parameters, and if any parameters were removed, then all the steps described are repeated.

Conjunction[edit]

Conjunction analysis allows the experimenter to ask "do two or more conditions activate the same voxels"? You may find two different kinds of conjunction analysis used:

  • Global null. The null hypothesis is that all conditions are null in a given voxel. Therefore, a single condition having an effect would be sufficient to disprove the null hypothesis. This is implemented in SPM when pressing the 'global' button by taking the minimum of two statistical images.
  • Conjunction null Given two conditions, A and B, the conjunction null hypothesis is that there is no activation for condition A OR no activation for condition B. So to disprove it, both condition A and B would have to be show an effect.