User:Peterz/sandbox

PEB

A common experimental aim is to test whether effective connectivity is different between groups of subjects, or is different according to a behavioural measure (e.g. test scores) within a group. One approach is to take DCM connectivity parameters and apply a classical test (e.g. t-tests, CVA). The disadvantage of this approach is that it throws away the estimated uncertainty (variance) about the connection strengths. Alternatively, one may construct a hierarchical model over the parameters - describing how group level effects constrain parameter estimates on a subject-by-subject basis. SPM12 includes a Parametric Empirical Bayes (PEB) model, which makes it possible to evaluate group effects and between-subjects variability on parameters.

Note: The GUI (batch) interface described here will be in the next public SPM release.

1. Specify your DCMs

Begin by specifying your DCMs in the normal way for each subject. For a tutorial on specifying DCMs, see the relevant chapter of the SPM manual.

2a. Assemble your DCMs in a group DCM file (GCM_name.mat) and estimate them (if not done already)

The next step is to assemble your DCMs into a cell array with one subject per row and one DCM per column, and estimate them. This can easily be done via the batch.

If your DCMs are not yet estimated:

• Click Batch from the main SPM window to open the batch editor.
• Click SPM -> DCM -> DCM estimation.
• Select your DCMs, either one model at a time ("Per model") or one subject at a time ("Per subject"), whichever is more convenient. If you have more than one DCM per subject, which is optional, ensure that the first model for each subject is a 'full' model containing all connections of interest. Any subsequent models should be 'nested' models, with some connections turned off, that you wish to compare.
• Select a directory and filename for the results to be stored.
• Check that you're happy with the estimation type (see the batch window for details on the options).
• Click the green play button.

If your DCMs are already estimated:

• Follow the steps above, but under estimation type, select "None (collate only)".

Here's how to achieve the same result using underlying Matlab functions:

% Collate DCMs into a GCM file
GCM = {'DCM_subject1_model1.mat','DCM_subject1_model2.mat';
       'DCM_subject2_model1.mat','DCM_subject2_model2.mat'};

% Fully estimate model 1
GCM(:,1) = spm_dcm_fit(GCM(:,1)); 

% Use Bayesian Model Reduction to rapidly estimated models 2-N
GCM = spm_dcm_bmr(GCM);

% Alternatively, replace the above lines with this to alternate between estimating
% DCMs and estimating group effects. This is slower, but can draw subjects out 
% of local optima towards the group mean.
% GCM = spm_dcm_peb_fit(GCM);

save('GCM_example.mat','GCM');

3. Estimate a second level PEB (Parametric Empirical Bayes) model

Having finished the first level analysis, we now create a second level (group) general linear model over the parameters:

• In the batch editor select SPM -> DCM -> Second level -> Specify / Estimate PEB.
• Give the analysis a name and select the GCM file from step 1.
• Leave DCM index set to 1 (we're going to build a group model using the full DCM from each subject).
• Under Covariates, add any between-subjects effects you'd like to model. (This is like specifying a design matrix in a regular second level SPM analysis). The first regressor will model the group mean and is added to the design matrix automatically. Some of the PEB functions assume that the first covariate you add will be of experimental interest, whereas subsequent covariates are nuisance variables of no interest. If you leave the covariates option blank, only the group mean will be modelled.
• Under Fields, select which fields of the DCM to model at the group level. We recommend keeping this to a small number of parameters. If, for instance, you're interested in the A-matrix and the B-matrix, then run a separate PEB analysis on each.

This will create a group level model named PEB_name.mat. It contains parameters representing each between-subjects effect on each DCM connection.

Here's the equivalent Matlab code using the underlying SPM functions:

% Specify PEB model settings (see batch editor for help on each setting)
M = struct();
M.alpha = 1;
M.beta  = 16;
M.hE    = 0;
M.hC    = 1/16;
M.Q     = 'single';

% Specify design matrix for N subjects. It should start with a constant column
M.X = ones(N,1);

% Choose field
field = {'A'};

% Estimate model
PEB     = spm_dcm_peb(GCM,M,field);

save('PEB_example.mat','PEB');

4a. Compare the full PEB model to nested PEB models to test specific hypotheses

The PEB model created above will have lots of parameters (the number of group effects multiplied by the number of DCM connections). To test hypotheses, you can compare this full PEB model to nested PEB models with certain parameters switched off (fixed at their prior mean of zero). For instance, you could switch off all group level effects, including the group mean, on the connection from Region R1 to Region R2. The difference in model evidence between the full PEB and this nested PEB, with R1->R2 parameters switched off, will tell you the importance of the R1->R2 at the group level.

To perform this analysis, the GCM file created above will need to contain a full DCM for each subject (column 1), as well as all nested DCMs you wish to test (columns 2-N). 'Nested' here means that one or more connections were switched off during model specification. These nested DCMs don't need to be estimated - they only serve to tell the software which group level PEB parameters to try switching off. To do this using the batch:

• In the batch editor select SPM -> DCM -> Second level -> Compare / Average PEB models.
• Select the PEB model you created in the previous step.
• For DCMs, select a GCM file specifying the models to be compared. This can be the GCM file you created earlier, or a new one.
• Click the green play button

Here's the equivalent Matlab code:

% Compare nested PEB models. Decide which connections to switch off based on the 
% structure of each DCM for subject 1.
BMA = spm_dcm_peb_bmc(PEB, GCM(1,:));

Two windows will be created. One titled BMC ("Bayesian Model Comparison") shows the results of comparing the full PEB model against the nested PEB models. The other window, titled "PEB - Review Parameters" shows the estimated connection strengths at the group level, averaged over PEBs (a Bayesian Model Average).

4b. Search over nested PEB models

Rather than compare specific hypotheses, as described above, you may wish to simply prune away any parameters from the PEB which don't contribute to the model evidence. This approach (previously referred to as post-hoc search) can be performed as follows.

• In the batch editor select SPM -> DCM -> Second level -> Search nested PEB models.
• Select the PEB file from earlier and the GCM file, which contains the first level DCMs.
• Press the green play button

Here's the equivalent Matlab code:

% Search over nested PEB models.
BMA = spm_dcm_peb_bmc(PEB);

Again, two windows will be created. One titled BMC ("Bayesian Model Comparison") shows the connections that were switched off - both for the group mean (commonalities across subjects) and first two between-subjects effects. The other window, titled "PEB - Review Parameters" shows the estimated group-level connection strengths, averaged over PEB models identified by the serach (a Bayesian Model Average).

5. Leave-one-out cross validation

Having identified one or more group effects in the results of the PEB analysis, you may wish to ask if the effect on a particular connection would be large enough to predict whether a new subject was in a particular group, or predict a continuous regressor, such as a clinical score. You can do this as follows.

• In the batch editor select SPM -> DCM -> Second level -> Predict (cross-validation).
• Specify the PEB model exactly as for step 3. This time, under Field, you may wish to include one or two connections you have previously found expresses a significant between-subjects effect.
• Press the green play button.

Alternatively, with Matlab code:

% Perform leave-one-out cross validation (GCM,M,field are as before)
spm_dcm_loo(GCM,M,field);

A PEB model will now be estimated while leaving out a subject, and will be used to predict the first between-subjects effect (after the constant column) in the design matrix, based on the specific connections chosen. The resulting plot shows the predicted group effect for each subject as well as the correlation between the predicted score and known score. If the effect to be predicted is binary (e.g. patient or control), then the bottom plots show on a subject-by-subject basis how confident one can be of the predicted group membership. If it's a continuous variable, like a clinical score, then the plot shows the prediction accuracy.