SPM/Parametric Empirical Bayes (PEB)

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

Running PEB for DCM[edit]

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.

Specify your DCMs[edit]

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.

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

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 DCMs 2-N for each subject if applicable
if size(GCM,2) > 1
   GCM = spm_dcm_bmr(GCM);
end

% Alternatively, replace the above lines with this code 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);

% Write results
save('GCM_example.mat','GCM');

Estimate a second level PEB (Parametric Empirical Bayes) model[edit]

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     = 'all';

% 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');

Compare the full PEB model to nested PEB models to test specific hypotheses[edit]

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(1), 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).

Search over nested PEB models[edit]

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(1));

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).

Review results[edit]

To review the results of a BMA analysis (or a PEB model), use the PEB review function:

% Review results
spm_dcm_peb_review(BMA,GCM)

Note that in this line of code, we have provided the group DCM (GCM) to the review function, which provides useful information such as the names of the connections.

Interpreting the output[edit]

The model

To recap, the PEB model (or Bayesian Model Average of PEB models) is a general linear model:

Where is the vector of all connectivity parameters from all the individual subjects' DCMs. Design matrix has regressors (columns) for the effect of each covariate on each connection. For example, there may be a regressor for the group mean of a particular connection and another regressor for the difference in connection strength per year of the subjects' age. The corresponding parameters quantify the effect of each covariate on each connection and these are estimated from the data. The residual between-subject variability (random effects) form the vector .

The parameters

The estimated PEB parameters are shown as bar charts in the spm_dcm_peb_review GUI, with one bar chart per covariate. More specifically, when the PEB model is estimated, a multivariate normal probability density over the parameters is computed, with expected values stored in the field BMA.Ep and covariance matrix stored in the field BMA.Cp. The height of the bars correspond to the expected values (Ep) and the error bars are 90% confidence intervals, computed from the leading diagonal of the covariance matrix (Cp).

Some of the parameters in encode the commonalities over subjects. These are either the group mean of the DCM connections over subjects, or the baseline connectivity, depending on whether or not the covariates were mean-centred (see Example design matrices, below). If covariates (between-subject differences) are included in the design matrix, then there will be parameters encoding the differences between subjects due to each covariate. A positive value indicates the connection is positively associated with the covariate. Conversely, negative values indicate a negative relationship between the connection and the covariate.

Posterior probabilities

The scheme also returns a probability for each parameter. When performing an automatic search over reduced PEB models, this is computed by comparing the evidence for all models in which the particular connection was switched on, versus the evidence for all models in which the connection was switched off (out of the best 256 models from the search). The same approach to computing probabilities is used if pre-defined models are provided, rather than performing an automatic search, except probabilities will only be computed for parameters relating to the commonalities and the first group difference. The GUI provides the option to threshold parameters based on these probabilities. Note that this is intended to help focus attention on the most probable effects - there is no concept of 'statistical significance' with Bayesian statistics, and it is a good idea to report the posterior probability of all non-trivial effects.

Leave-one-out cross validation[edit]

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.

Example design matrices[edit]

Here we give some examples of how to define the between-subjects design matrix, M.X, for some typical experimental designs. All the same principles for the general linear model (GLM) also apply in the PEB scheme, the only restriction being that the PEB software expects the design matrix to start with a column of ones. For each type of experimental design, there are multiple options for how to encode the design matrix. Part of the decision of which design to use will be based on interpretability - for example, do you want parameters encoding the difference between groups, or the mean of each group separately? Additionally, different designs will have different degrees of statistical efficiency. For example, if the regressors are orthogonal (statistically independent), then efficiency will be maximised. If in doubt about which design is best, perform a Bayesian model comparison by trying different options and choosing the one which gives the most positive free energy PEB.F.

One between-subjects factor (two groups)[edit]

This can be modelled with two regressors, encoding 1) the group mean and 2) group difference relative to the mean:

Here, the rows are the four subjects and the columns are the regressors. The first two subjects are in group 1 (indicated by a value of -1 in the second regressor) and the second two subjects are in group 2 (indicated by value 1 in the second regressor). Note that the second regressor must have a mean of zero for this interpretation to hold. If the group difference regressor doesn't have mean of zero, you can set this manually:

X(:,2:end) = X(:,2:end) - mean(X(:,2:end))

The ensuing PEB model will have parameters encoding the group mean of each connection strength (due to the first regressor) as well as parameters encoding the deviation from the mean due to the group difference (thanks to regressor 2). For the group difference, positive estimated parameters indicate stronger connectivity in group 2 than group 1 and negative parameters indicate stronger connectivity in group 1 than group 2.

Alternatively, if the second regressor is not mean-centred, then the regressors are: 1) the connectivity of group 1, and 2) the difference between groups:

One between-subject factor (three groups)[edit]

Consider three groups of subjects, with two subjects per group. If a linear effect of group is hypothesised - i.e. that group 1 > group 2 > group 3 - this could be modelled using two regressors: 1) the overall mean and 2) the difference between groups 1 and 3:

Alternatively, if the objective is to test for the commonalities and between-group differences (without specifying a linear effect), three regressors can be used to encode: 1) the overall mean 2) the difference between the first and second groups, and 3) the difference between the second and third groups:

Or, to model each group separately, a group can be chosen to serve as the baseline (here, group 1), and the regressors encode: 1) the first group, 2) the additive effect of being in the second group relative to the first group, and 3) the additive effect of being in the third group relative to the first group:

Unbalanced design with two between-subject factors (three groups)[edit]

Consider an experiment with three groups: A) patients receiving a drug B) patients receiving a placebo, and C) a single control group. This is an example of an unbalanced design. Although the designs listed in the previous section could be used, it may be more intuitive to use the following three regressors: 1) baseline (controls), 2) the additive effect of being a patient, and 3) the difference between drug and placebo. With 2 subjects per group:

Note that these regressors are not orthogonal, so this will be a less efficient experimental design than a fully balanced factorial design, which is described next.

Balanced factorial design with two between-subject factors (four groups)[edit]

Consider an experiment with four groups of subjects in a balanced factorial design: 1) patients receiving a treatment, 2) patients receiving a placebo, 3) controls receiving a treatment and 4) controls receiving a placebo. The regressors will encode: 1) the overall mean, 2) the main effect of being a patient, 3) the main effect of treatment and 4) the interaction, i.e. the difference in drug effect between patients and controls. With two subjects per group:

Here, the final column - the interaction - is generated by element-wise multiplication of the two main effects (which must be mean-centred first). Note that this design gives rise to regressors that are orthogonal - making the factorial design optimally efficient.

Covariates[edit]

In addition to the regressors described above, it is common to have covariates, e.g. age or clinical scores. These can be added to the design matrix, and it is generally a good idea to mean-centre them (across all subjects) to enable the first regressor to be interpretable as the mean. Note that for every covariate added to the design matrix, one parameter is added per DCM connection, so you can end up with a lot of parameters. Therefore, it is a good idea to be conservative and keep the number of covariates to a minimum.

Hierarchical experimental designs[edit]

Consider an experiment with a factorial design at the between-run level. For example, we will imagine a study with two groups of subjects, where one group received a drug and the other received a placebo, and each subject was scanned twice: pre- and post-treatment. This is a balanced 2x2 design with cells or experimental conditions:

  • Group 1 pre-treatment (GCM1_pre)
  • Group 1 post-treatment (GCM1_post)
  • Group 2 pre-treatment (GCM2_pre)
  • Group 2 post-treatment (GCM2_post)

We'll assume that a separate DCM was fitted for each of these four experimental conditions for each subject. The names in brackets are example variable names to store the DCMs for each condition. Each variable is a cell array with one subject per row. Here are two ways to model a design of this sort using PEB.

Option 1: A 2-level design (first level: DCM, second level: PEB)[edit]

Create a PEB model in which the design matrix X has the following regressors:

  1. Overall mean
  2. Main effect of group (-1s for group 1, 1s for group 2, and then mean-corrected)
  3. Main effect of treatment (-1s for pre-treatment, 1s for post-treatment, and then mean-corrected)
  4. Interaction of group and treatment (the two mean-corrected main effects element-wise multiplied)

And fit this PEB model to all DCMs across subjects and time points (in the right order to match the regressors you create), e.g:

GCM = {GCM1_pre; GCM1_post; GCM2_pre; GCM2_post};
PEB = spm_dcm_peb(GCM, X);

An additional technical consideration arises if you wish to compare the evidence for pre-defined alternative PEB models, because the evidence for each model will only be compared in terms of the first two regressors (in this example, the overall mean and the main effect of group). To compare pre-defined models for all of these 4 four factors, simply re-run the PEB analysis with design matrix re-ordered so that the effect of interest is the second regressor. For example, to investigate which model best explains the effect of treatment, re-order the columns to: Mean, Treatment, Group, Interaction and then run spm_dcm_peb and spm_dcm_peb_bmc (or use the batch).

Option 2: A 3-level design (first level: DCM, second and third levels: PEB)[edit]

Alternatively, the same design can be modelled in a 3-level hierarchy, using the PEB-of-PEBs approach, where a PEB model is created for each group separately. For group 1, this will have a design matrix X1 with regressors:

  1. Mean of group 1
  2. Effect of treatment on group 1 (-1s for pre-treatment, 1s for post-treatment, and then mean-corrected)

And this will be fitted to all DCMs from group 1 in a single column vector:

GCM1 = {GCM1_pre; GCM1_post};
PEB1 = spm_dcm_peb(GCM1, X1);

For group 2, the PEB model will have a design matrix X2 with regressors:

  1. Mean of group 2
  2. Effect of treatment on group 2 (-1s for pre-treatment, 1s for post-treatment, and then mean-corrected)

And this will be fitted to the DCMs of group 2:

GCM2 = {GCM2_pre; GCM2_post};
PEB2 = spm_dcm_peb(GCM2, X2);

The ensuing PEB models PEB1 and PEB2 will have parameters encoding each of the the two effects (mean and treatment) on each DCM connection. Now, take the parameters of the PEB models up to the third level of the hierarchy, to identify commonalities and differences across groups. The design matrix X3 will contain regressors:

  1. Mean of all subjects
  2. Group difference (-1s for group 1, 1s for group 2, and then mean-corrected)

And this is fitted to the PEB parameters from the 2nd level:

PEBs = {PEB1; PEB2};
PEB3 = spm_dcm_peb(PEBs,X3);

The final model 3rd level PEB model, named PEB3, will include parameters for the group mean and group difference on each of the 2nd level effects. i.e. for every connection there will be a parameter for:

  • The overall mean connectivity
  • The main effect of treatment
  • The main effect of group
  • The interaction of group and treatment

Which option to use?[edit]

Option 1 is simpler, but option 2 might offer a better model of between subjects variability (random effects). To assess which is the better option, you can try both and compare the free energy. To do this, subtract the free energies of each model: PEB.F - PEB3.F . A positive number indicates option 1 is better, a negative number indicates option 2 is better. And make sure to feed back to the SPM mailing list on what you find - this will be useful for developing best practice.

Further reading[edit]