Wonderful ! Thanks a lot Peter and Adeel for all your help and time. It seems my analysis is working fine now.

Dear Sahil

Apologies for the delay in replying.

 

Thanks a lot for all the suggestions. So after I simplified the connectivity structure model and increased the number of iterations from 64 to 128, I am getting better connectivity parameters.

Attached figure shows the connectivity matrix A after I compared two groups: Controls +1 and patients -1.

 

I have some follow up questions:

(1). Is my interpretation that the connections whose values are between 0 and +1 are significantly stronger for controls than patients, and the connections whose values are between 0 and -1 are weaker in controls than patients? If that's correct, where can I find the details of statistical tests used here for publication purpose?

 

Nearly - a couple of minor tweaks to the way you stated the results. I’d avoid the word ‘significant’ when working with Bayesian statistics – you are able to state the probability for the effect of being a patient vs a control (get this by clicking on a bar). Significance is a concept from classical statistics. Also, the negative bars tell you that the connection is more negative in patients than controls. It could be that they are ‘weaker’ – i.e. closer to zero - but alternatively it could be that they are more inhibitory. These positive and negative effects of being a patient add to or subtract from the mean connectivity across groups– so use the GUI to have a look at that too. So you could say something like:

 

We took all subjects’ connectivity parameters to the group level and modelled them using a (Bayesian) GLM with regressors for the group mean and diagnosis per connection. We then used an automatic search to compare the evidence for this full GLM against hundreds of reduced GLMs, in which combinations of parameters had been switched off (Bayesian Model Reduction). We computed the average of the connectivity parameters of the best 256 models from this search (a Bayesian Model Average, BMA). The probability for each individual connectivity parameter was computed by comparing the evidence for all models (from the final 256) with the connection switched on, versus all models with the connection switched off. For convenience, we only show parameters exceeding 99% posterior probability.

 

The figure shows the BMA of the connectivity parameters. Positive bars indicate more positive connection strengths in patients than controls, and vice versa. Pink error bars are 90% confidence intervals. Between-region connections are in units of Hz, whereas self-connection parameters are unitless log scaling parameters, and multiply up or down the default inhibitory self-connection of -0.5Hz, according to the function:

 

A_hz = -0.5 * exp(A)

 

(2). I am not sure what is the unit of colorbar here, does this just represent the correlation coefficient?

 

No, there are no correlation coefficients here. This is in Hz  (except for the self-connections – see above) because all the connectivity parameters are rates of change (rate constants).

 

(3). This matrix, A, is generated at threshold settings: Free energy (with vs. without) at very strong evidence of pP > 0.99. I am not sure if these settings are correct and are commonly used. If so, is there any documentation available to understand these thresholds in more detail.

 

Perfect (although 99% is quite strict – one notch down may be appropriate). I have written a tutorial paper on all this and it’s with co-authors – looking forward to sharing it ASAP.

 

(4). Lastly, further I am interested in finding whether these particular significant connections (from matrix A) differ between pre- to post- treatment condition. How can I specify these particular connections for pre- and post-treatment data, rather than full connected model.

 

Perhaps you could just include treatment in your GLM, as you have diagnosis, and include all DCMs in the analysis?

 

Best

Peter

 

 

 

 

On Tue, Aug 21, 2018 at 2:03 AM, Zeidman, Peter <[log in to unmask]> wrote:

Dear Sahil

(1). I recommend going back to your DCMs per subject and looking at the A-matrix parameters (with pink error bars). Did many subjects get non-trivially large values? You can do this using spm_dcm_fmri_csd_results(DCM).

I have around 65 subjects. I was wondering if there is any script to check to see how many subjects have large values. I checked a few but I do not see any pink error bars in those subjects. Here I am attaching output window from one of those subjects.

 

No but you could write one to collate the largest between-region connection from each subject, e.g.

 

maxA = zeros(nsubjects,1);

for i = 1:nsubjects

   filename = sprint(‘DCM_subject%d.mat’,i);

   load(filename);

   A = DCM.Ep.A;

   maxA(i) = max(max(abs(A - diag(diag(A)))));

end

 

disp(maxA);

  

 

(2). For the group level model, after how many iterations did the PEB model converge or stop when you ran spm_dcm_peb?

PEB model converges after VL iteration 64:

VL Iteration 64      : F = 431180.57 dF: 36.0154  [+2.00]

 

64 iterations is the default limit for the PEB estimation. If the free energy is continuing to increase on each iteration when it reaches iteration 64 – i.e. if the dF (difference in free energy) is greater than 1 – then 64 isn’t enough iterations. Try increasing this limit. In the next version of SPM there’s a setting for this. For now, edit spm_dcm_peb.m and change the line around number 419 from:

 

for n = 1:64

 

to

 

for n = 1:128

 

Or you could try 256 iterations if that still isn’t enough. Make sure to write in your paper that you made this change.

 

(3). Perhaps you could share a DCM from a single subject so I can have a look at your timeseries?

Here, I am attaching 4 DCM files (dropbox link) - two for patients group and two for controls: https://www.dropbox.com/s/4fogxv8p03931sq/DCM_Files.zip?dl=0.

 

No obvious mistakes there. If the above doesn’t help, consider whether you could simplify the connectivity structure of the model at all to improve estimation efficiency, rather than having it fully connected?

 

Best

Peter

 

On Mon, Aug 13, 2018 at 9:02 AM, Zeidman, Peter <[log in to unmask]> wrote:

Dear Sahil

 

- After I include age and sex as covariates as we discussed above, it runs successfully. It gives 4 covariates: 1st for group mean, 2nd for group difference after controlling for age and sex, and if that's correct, then I am not sure what does 3rd and 4th covariate represent here?

 

The covariates in the GUI correspond to the order of covariates in your design matrix M.X. So if the 4 columns in your design matrix are:

 

1.       group mean

2.       group difference

3.       age

4.       sex

 

Then yes, the 4 covariates in the GUI will be in this order. As for any GLM, the parameters (betas) you get associated with each covariate will be the variance which can be attributed to that covariate after accounting for all the others.

 

- When I ran BMA = spm_dcm_peb_bmc(PEB(1)); and review my results using spm_dcm_peb_review(BMA,GCM), again I am getting very low posterior expectation values (~0.0001) (for all 4 covariates) in the connectivity matrix plot for A. On the other hand, output in the matlab window shows following:

13 models in Occams window:

            Model 1, p(m|Y)=0.00

            Model 2, p(m|Y)=0.00

            Model 3, p(m|Y)=0.00

            Model 4, p(m|Y)=0.10

            Model 5, p(m|Y)=0.00

            Model 6, p(m|Y)=0.03

            Model 7, p(m|Y)=0.01

            Model 8, p(m|Y)=0.00

            Model 9, p(m|Y)=0.00

            Model 10, p(m|Y)=0.00

            Model 11, p(m|Y)=0.00

            Model 12, p(m|Y)=0.00

            Model 13, p(m|Y)=0.87

 

 i.e model 13 with maximum posterior probability. I am not sure what does my results represent here, and where are my results in terms of winning model and winning parameters showing the group difference between controls and patients? Its also not clear from the wikipage that how to drop off parameters to test specific hypotheses.

 

I agree that it’s  not a good sign you only have small connectivity parameters. This shows that only a small amount of the variance could be ascribed to connectivity between neural populations. The issue is most likely at the first level – there’s either no effect for the DCM to explain, or the priors in the DCM aren’t suitable for your data. Less likely, there could be an issue at the group level – there could be inappropriate priors in the PEB model. I recommend going back to your DCMs per subject and looking at the A-matrix parameters (with pink error bars). Did many subjects get non-trivially large values? You can do this using spm_dcm_fmri_csd_results(DCM).

 

For the group level model, after how many iterations did the PEB model converge or stop when you ran spm_dcm_peb_bmc?

 

Regarding the model comparison output. The automatic search over reduced models will try hundreds of combinations of switching off connections to find the optimal combination. The best 256 models from the final iteration of the search will then be compared and their parameters averaged (Bayesian Model Average). The models in Occams window above are the best of the best – so clearly your BMA is mainly being driven by model 13. To see which connections were switched on or off in this model, re-run the search with the second output BMR:

 

[BMA,BMR]=spm_dcm_peb_bmc(PEB);

K = full(BMR.K);           % Matrix containing which parameters were turned OFF in each model

on = ~K;                           % Matrix containing which parameters were turned ON in each model

name = BMR.name;   % Names of parameters which varied across models

disp('Parameters switched ON in model 13: ');

disp(name(on(13,:)))

 

Basically, I am not sure how to proceed from here onwards. I specified full model containing 9 ROIs and my aim is to compare two groups and to identify the model and parameters which are stronger is one group compared to other.

I tried to find answers to all the above questions and interpretation of my current results from SPM discussion forum but couldn't find.

 

To summarise, if you find that your DCMs consistently have tiny A-matrix parameters across subjects (all less than 1/8), then it’s most likely that something is not right (or not typical) with the data, the preprocessing or the GLM analysis. I say that because the DCMs have priors which we find work in general. Perhaps you could share a DCM from a single subject so I can have a look at your timeseries? (Please zip it, as my outlook doesn’t allow .mat files.)

 

Best

Peter

 

 

Thanks a lot for your help and time.
Sahil

 

 

On Wed, Aug 8, 2018 at 11:43 PM, Adeel Razi <[log in to unmask]> wrote:

Dear Sahil,

 

1. yes you are right that since you have two covariates, the indexing of the parameters will change as you noted below. It would be easier if you use the GUI for plotting/visualization. And just to clarify what you are referring to as posterior probabilities are posterior expectations of the parameter estimates

2.  All good, only that you meant [1 1 25 1;1 1 22 1;1 -1 23 0;1 -1 26 1]? I would also advise to mean correct your continuous variables like age. Think of this as a usual design matrix, the group differences will take into account variables of no interest (age/gender).

 

Best wishes,

Adeel

 

On Wed, Aug 8, 2018 at 3:49 AM, Sahil Bajaj <[log in to unmask]> wrote:

Dear Dr. Razi,

 

Thank you so much for your reply. Yes, I redefined: M.X = [1 1;1 1;1 -1;1 -1]; to compare between two groups and I am getting additional covariate - covariate 2, as you described. So it seems its working fine. However, I have following doubts:

 

(1). Since posterior probability values are so low, so I set the limits of Ep as following:

 

Ep = reshape(BMA.Ep,9,9);

imagesc(Ep,[-0.02 0.02]), but I am getting the following error, because the dimensions of Ep are 162x1.

 

Ep = reshape(BMA.Ep,9,9);

Error using reshape

To RESHAPE the number of elements must not change.

 

I guess first 81 values are for covariate 1 and next 81 are for covariate 2. If so, then reshape(BMA.Ep(82:161),9,9); should be the correct adjusted Ep for 9x9 matrix for covariate 2. Correct?

 

(2). I am also interested in using age and sex as covariates. Previous post in SPM forum suggested to edit M.X matrix as M.X = [ones(N,1) age gender], which in my case will be [M.X] = [1 1;1 1;1 -1;1 -1; 25 1;22 1;23 0;26 1], where age and sex for subject 1 is: 25 1; for subject 2: 22 1; for subject 3: 23 0 and for subject 4: 26 1. Here 0 is for males and 1 for females.

 

Could you please just confirm if that's correct? If so, then again I am getting two covariates in output window - covariate 1 and covariate 2, does covariate 2 here represent the difference between two groups after eliminating the effect of age and sex? Or covariate 2 is just the group difference without accounting for age and sex, and there are some extra steps to include them as covariates? I am not sure what are those steps?

 

I really appreciate all your help and time.

 

Thanks,
Sahil

 

 

On Fri, Aug 3, 2018 at 3:07 PM, Adeel Razi <[log in to unmask]> wrote:

Dear Sahil,

I had a look at your data and the outputs. There are connectivity parameters in the middle panel ('Estimated Parameters') when you review your results. They are not flat lines if you look closely. The pink bars are 90% confidence intervals whereas the grey bars are the posterior expectations.

There is one covariate which represents the mean connectivity across groups as specified by the code:

PEB.M.X= ones(N,1)

If you are interested in group differences then add second column of 1 and -1 in the design matrix PEB.M.X for your two groups B and G. You can then view the group differences using the drop down menu by choosing covariate 2.

There are plenty of previous messages from Peter that provide lots of help in terms of interpretation of PEB outputs if you search through the archived messages in the mailing list.

I hope this helps.

Best wishes,

Adeel

 

 

On Sat, Aug 4, 2018 at 3:48 AM, Sahil Bajaj <[log in to unmask]> wrote:

Dear Drs. Adeel Razi and Peter Zeidman (and SPM/DCM experts),

 

I am trying to perform large-scale DCM analysis for resting-state fMRI data. I have 9 regions of interest - defined based on MNI coordinates from literature.

 

I am interested in comparing the connectivity between these 9 regions between two groups: B and G.  So I ran following steps (just to test, I used two subjects per group) using PEB wiki guidelines:

 

-------------------- Steps I ran:-------------------- 

% Two Groups: B and G

% DCM_DMN_B1.mat: Full model for subject 1 for group B

% DCM_DMN_B2.mat: Full model for subject 2 for group B

% DCM_DMN_G1.mat: Full model for subject 1 for group G

% DCM_DMN_G2.mat: Full model for subject 2 for group G

 

GCM = {'DCM_DMN_B1.mat';'DCM_DMN_B2.mat';'DCM_DMN_G1.mat';'DCM_DMN_G2.mat'}; 

GCM(:,1) = spm_dcm_fit(GCM(:,1));

 

GCM = spm_dcm_peb_fit(GCM);

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

 

% Specify PEB model settings

N = 4; % number of subjects

M = struct();

M.alpha = 1;

M.beta  = 16;

M.hE    = 0;

M.hC    = 1/16;

M.Q     = 'single';

 

% Specify design matrix for N subjects

M.X = ones(N,1);

 

% Choose field

field = {'A'};

 

% Estimate model

PEB     = spm_dcm_peb(GCM,M,field);

 

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

 

BMA = spm_dcm_peb_bmc(PEB(1));

spm_dcm_peb_review(BMA,GCM)

------------------------------------------------------------------------------ 

 

 

This runs successfully, but my output is showing a blank connectivity matrix (no connectivity between any pair of nodes among 9), even at no threshold. For your convenience, here I am attaching one *.mat file (DCM_DMN_B1.mat) and final output screen shot. Also, its showing that there is one covariate, I am not sure what is that, as I didn't define any covariate in above steps.

 

Could you please confirm if the above steps are correct and help me in figuring out what I am doing wrong here while comparing two groups: B and G? Also, any documentation to interpret the final results will be really useful.

 

Thank you so much !

 

Best,
Sahil