Dear experts,
Just a quick kind reminder regarding my previous email about large scale
DCMs, I would really appreciate your response on it.
Thanks,
Sahil
---------- Forwarded message ---------
From: Sahil Bajaj <[log in to unmask]>
Date: Wed, Aug 22, 2018 at 4:43 PM
Subject: Re: [SPM] Fwd: [SPM] Large-scale DCM (PEB) help and confirming the
steps
To:
Cc: [log in to unmask] <[log in to unmask]>
Dear Peter and Adeel,
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?
(2). I am not sure what is the unit of colorbar here, does this just
represent the correlation coefficient?
(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.
(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.
Thanks again for the help.
Best,
Sahil
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
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
