Dear Delin
What’s in your design matrix (M.X) ? It is expected that you have a column of ones (modelling the group mean), then your group regressor
of interest, then any other regressors of no interest. Also note that this works best with a small number of parameters – so you may want to include a few specific connections – e.g. {‘A(1,2)’,’A(3,2)’}.
If you have no clear winning model, then it’s fine to report that, and report the Bayesian Model Average.
Best
Peter
From:
孙得琳 [mailto:[log in to unmask]]
Sent: 03 February 2017 18:41
To: Zeidman, Peter <[log in to unmask]>
Cc: [log in to unmask]
Subject: Re: [SPM] Questions about fMRI DCM
Dear Peter,
Thank you so so so much for the detailed replies!
However, I ran "spm_dcm_loo(GCM(:,1),M,{'A'})" (also for 'B' and 'C') but still got errors after a few steps of iterations:
------------------------------------------------------------------
VL Iteration 1 : F = 31.78 dF: 0.0000 [-3.75]
VL Iteration 2 : F = 37.09 dF: 5.3100 [-3.50]
VL Iteration 3 : F = 44.22 dF: 7.1269 [-3.25]
VL Iteration 4 : F = 52.70 dF: 8.4854 [-3.00]
VL Iteration 5 : F = 59.44 dF: 6.7380 [-2.75]
VL Iteration 6 : F = 65.66 dF: 6.2218 [-2.50]
VL Iteration 7 : F = 71.58 dF: 5.9144 [-2.25]
VL Iteration 8 : F = 77.14 dF: 5.5650 [-2.00]
VL Iteration 9 : F = 82.26 dF: 5.1164 [-1.75]
VL Iteration 10 : F = 86.82 dF: 4.5637 [-1.50]
VL Iteration 11 : F = 90.75 dF: 3.9270 [-1.25]
VL Iteration 12 : F = 93.99 dF: 3.2379 [-1.00]
VL Iteration 13 : F = 96.52 dF: 2.5360 [-0.75]
VL Iteration 14 : F = 98.40 dF: 1.8743 [-0.50]
VL Iteration 15 : F = 99.71 dF: 1.3178 [-0.25]
VL Iteration 16 : F = 100.63 dF: 0.9148 [+0.00]
VL Iteration 17 : F = 101.29 dF: 0.6595 [+0.25]
VL Iteration 18 : F = 101.78 dF: 0.4947 [+0.50]
VL Iteration 19 : F = 102.14 dF: 0.3606 [+0.75]
VL Iteration 20 : F = 102.38 dF: 0.2343 [+1.00]
VL Iteration 21 : F = 102.50 dF: 0.1242 [+1.25]
VL Iteration 22 : F = 102.55 dF: 0.0451 [+1.50]
VL Iteration 23 : F = 102.55 dF: 0.0022 [+1.75]
VL Iteration 24 : F = 102.54 dF: 0.0022 [+0.75]
VL Iteration 25 : F = 102.55 dF: 0.0022 [-0.25]
VL Iteration 26 : F = 102.55 dF: 0.0022 [-1.25]
VL Iteration 27 : F = 102.55 dF: 0.0022 [-2.25]
VL Iteration 28 : F = 102.55 dF: 0.0022 [-3.25]
VL Iteration 29 : F = 102.55 dF: 0.0022 [-4.25]
VL Iteration 29 : F = 102.55 dF: 0.0022 [-4.25]
Index exceeds matrix dimensions.
Error in spm_dcm_ppd (line 88)
bC = var(X(:,iX))*4;
Error in spm_dcm_loo (line 82)
[Ep,Cp,P] = spm_dcm_ppd(DCM(i),DCM(j,1),M.X(i,:),M.X(j,:),field);
--------------------------------------------------------------------
May I know what's wrong it is? Thanks a lot!
BTW, if there is no winning model in my data (as you mentioned), can I just report the model with the highest posterior probability (about 80%)?
Bests,
Delin
On Fri, Feb 3, 2017 at 11:21 AM, Zeidman, Peter <[log in to unmask]> wrote:
Dear Delin
(1) for step 3 "Estimate a second level PEB (Parametric Empirical Bayes) model", could you please explain more about the meaning of "M" and its fields such as M.alpha, M.beta and M.hC?
Sure. Note that these are documented in the batch editor in recent releases of SPM12 (click SPM -> DCM -> Second level -> Specify / Estimate PEB and look at the help text at the bottom).. But I’m happy to recap. M stands for ‘model’ and it contains the priors for the PEB model. It also contains the design matrix, M.X.
M.alpha
=======
Each parameter in the PEB model represents a group-level effect on a connection. The parameters are normal distributions, which are specified by a mean and (co)variance. The mean of each parameter (PEB.Ep) is the estimated effect size and the variance (PEB.Cp) is the uncertainty. Each parameter also has a prior normal distribution, which is what you believe about the parameters before doing the analysis. If I remember correctly, the prior mean of each parameter is a small value close to zero, which says that in the absence of evidence, there’s no group level effect. The prior variance is set to a fraction of the priors in the first level DCM analysis, i.e.
prior PEB parameter variance = (1 / alpha) * (prior DCM parameter variance)
By default, alpha = 1. This says that your beliefs about the connections taking on a non-zero value are the same at the group level as at the single subject level.
M.beta and M.hE
=======
As well as estimating the group-level parameters, the PEB model also estimates the between-subject variability (random effects precision). The prior expectation of the random effects precision is set to:
prior RFX precision = 1 / (prior DCM parameter variance / M.beta)
By default, M.beta=16. In other words, we expect between-subject variability for each parameter to be 1/16 of our uncertainty about each parameter at the single-subject level. That may sound a bit confusing, but it’s just a quick way to choose some reasonable prior belief about between-subjects variability. This prior is scaled by a parameter, which is estimated from the data. The prior for this scaling parameter is a normal distribution with mean M.hE and variance M.hC. Altogether, the prior covariance of between-subjects random effects is:
exp(M.hE)*DCM{1}.M.pC/M.beta
By default, exp(M.hE) = 1. In the line above, DCM{1}.M.pC is the prior covariance matrix of the first subject’s DCM.
So in summary, M.alpha controls the prior belief about parameters, M.beta and M.hE control the prior belief about between-subjects variability.
(2) for step 4a "Compare the full PEB model to nested PEB models to test specific hypotheses", I ran "BMA = wpm_dcm_peb_bmc(PEB(1),GCM(1,:))" for field A but did not get any outputs. I guess maybe it is because the A field is the same across models (i.e. amygdala and PFC are always mutually connected). Am I right?
If there’s no variability in your A-matrix across models, and you only ran the PEB model on the A-matrix parameter, then I imagine you may not get output…
(3) for step 4a "Compare the full PEB model to nested PEB models to test specific hypotheses", I ran "BMA = wpm_dcm_peb_bmc(PEB(1),GCM(1,:))" for field B and C and got two figures (links shown below). Could you please tell me which parameters should be reported in paper? Can I state theta model 2 is the winning model for field B?
field C, https://dl.dropboxusercontent.com/u/100177332/4a%20Field%20C.png
You could write in your paper something like:
Having estimated a group level (PEB) model, we now tested our hypotheses by comparing the full PEB model to nested models where certain group effects had been switched off (fixed at their prior mean of zero). We did this separately for the B matrix (modulations) and the C matrix (driving inputs).
The top left of the figure shows our model space. Each column represents a model or hypotheses, and the rows indicate which group level parameters where switched on (white) or switched off (black) for each model. The middle left panels shows the posterior probability for each nested model, where model 1 is the full model. There was positive evidence for model 2 (posterior probability ~0.8). The second strongest model was model 3 (posterior probability XX). Because there was no single winning model, we performed Bayesian model averaging. Parameters B(1,1,2) and B(2,2,2) had negative values which we could be confident had moved away from their prior values (bottom right). The probability for each parameter was calculated by comparing the evidence for models which included each parameter against those which did not have that parameter. This is shown at the bottom left, and confirms that B(1,1,2) and B(2,2,2) had non-zero values with posterior probability > 0.95.
You can write something similar for the C-matrix results :-)
(3) for step 4a "Compare the full PEB model to nested PEB models to test specific hypotheses", about the command "BMA = wpm_dcm_peb_bmc(PEB(1),GCM(1,:))", does it mean that only the PEB of the first model (i.e. PEB(1)) is compared to GCM of all models of subject 1 (i.e. GCM(1,:))? should I also compare PEB(2), PEB(3) ... to GCM(2,:), GCM(3,:)...?
The idea is to run PEB only on the full model -i.e. PEB(1). Then you compare this PEB against nested PEBs. The software uses DCMs 2 to N to tell it which connections to switch off in each nested PEB. So DCMs 2 to N don’t even need to be estimated.
(5) for step 4b "Search over nested PEB models", I ran "BMA = wpm_dcm_peb_bmc(PEB(1))" and got outputs (links shown below). However, I am confused because there are 18 models (I only have 10 models [1 full connected model + 9 hypothesized models]). Did I make any mistake?
field C, https://dl.dropboxusercontent.com/u/100177332/4b%20field%20C.png
It doesn’t know about your hypothesized models. It only knows about your full model - PEB(1). The 18 models it shows have been discovered automatically – they’re different combinations of switching off conections.
(6) for step 5 "leave-one-out cross validation", I ran "wpm_dcm_loo(GCM,M,{'A'})" (also tried {'B'} and {'C'}) but got error messages:
Index exceeds matrix dimensions. What should I do?
Error in spm_dcm_loo (line 77)
[Ep,Cp,P] = spm_dcm_ppd(DCM(i),DCM(j,1),X(i,:),X(j,:),field);
Error in spm_dcm_loo (line 60)
[p,q,r] = spm_dcm_loo(DCM(:,i),X,field);
Try just running this on the full model for each subject:
spm_dcm_loo(GCM(:,1),M,{'A'})
(7) May I know what parameters I can use to do between-group comparisons and to correlate with individual measures (e.g. age, behavioral performance)?
One option is to put these in the PEB. It will then give you the influence of each covariate on each connection. You do this in the first step above – spm_dcm_peb – by including your covariates in the design matrix M.X. Make sure that the first column is the groups mean: M.X = [ones(N,1) age gender] etc.
If you need the DCM parameters for correlating against something else, externally from SPM, use the DCM connections informed by the group. i.e.
[PEB, GCM] = spm_dcm_peb…)
The second output of spm_dcm_peb are the DCMs informed by the group.
Best
Peter
Thank you so much!
Bests,
Delin
On Wed, Feb 1, 2017 at 4:57 AM, Zeidman, Peter <[log in to unmask]> wrote:
Dear Delin
I am having difficulty working out which figures come from which step of your analysis. Perhaps you have given them informative filenames – however I cannot see the filename, probably because you dragged or copy/pasted the images into the email software. Please you could re-send with these images as attachments, with filenames that identify which step of your analysis each one arises from?
Best
Peter
From: 孙得琳 [mailto:[log in to unmask]]
Sent: 31 January 2017 18:19
To: Zeidman, Peter <[log in to unmask]>
Cc: [log in to unmask]
Subject: Re: [SPM] Questions about fMRI DCM
Dear Peter,
Thank you very much! The revision works!
(1) The output are attached as "4A Field B" and "4A Field C" (There is no output for Field A: I guess maybe because the A filed are the same in all of my models, i.e. amygdala and PFC are always mutually connected). Could you please tell me which parameter is important to be reported in these figures? Can I say that model 2 is the winning model for Field B?
(2) I also run the section of "4b. Search over nested PEB models" (Please see the attached figures). However, I am confused because there are 18 models (I only have 10 models [1 flu connected model + 9 hypothesized models]).
(3) I met problems when running the section of "leave-one-out cross validation":
Index exceeds matrix dimensions.
Error in spm_dcm_loo (line 77)
[Ep,Cp,P] = spm_dcm_ppd(DCM(i),DCM(j,1),X(i,:),X(j,:),field);
Error in spm_dcm_loo (line 60)
[p,q,r] = spm_dcm_loo(DCM(:,i),X,field);
(4) I am unclear about the meaning of "M". May I ask for your explanation about M.alpha, M.beta ...?
Thank you so much! I am sorry for interrupting you with so many questions.
Bests,
Delin
On Tue, Jan 31, 2017 at 11:49 AM, Zeidman, Peter <[log in to unmask]> wrote:
Dear Delin
My apologies, I missed your previous email. That’s a strange error. Please could you try the following:
BMA = spm_dcm_peb_bmc(PEB(1), GCM(1,:))
If the problem still occurs, please could you send me all the code or the batch .mat file you are using to run this?
Regarding the figures you mention below, that’s the model comparison for a single subject. You can ignore that for your group analysis.
Best
Peter
From: 孙得琳 [mailto:[log in to unmask]]
Sent: 31 January 2017 16:23
To: Zeidman, Peter <[log in to unmask]>
Subject: Fwd: [SPM] Questions about fMRI DCM
Dear Peter,
I hope this email finds you well. I met a problem when running DCM PEB. May I ask for your help? The error messages and mid-way outputs please see below. Thank you very very much!
Bests,
Delin
---------- Forwarded message ----------
From: 孙得琳 <[log in to unmask]>
Date: Wed, Jan 25, 2017 at 2:10 PM
Subject: Re: [SPM] Questions about fMRI DCM
To: "Zeidman, Peter" <[log in to unmask]>
Cc: "[log in to unmask]" <[log in to unmask]>Dear Peter,
Thank you for the suggestions on PEB. I tried it after adding a fully connected model (including [A] intrinsic connectivities of amygdala->amygdala, PFC->PFC, amygdala->PFC and amygdala<-PFC; [B] task modulations on connectivity: amygdala->amygdala, PFC->PFC, amygdala->PFC and amygdala<-PFC; [C] input to both amygdala and PFC) as suggested. It works fine before the step of "Compare the full PEB model to nested PEB models to test specific hypotheses", but showed error message when I ran "BMA = spm_dcm_peb_bmc(PEB, GCM(1,:))". The error is as below:
Error using size
Too many output arguments.
Error in spm_dcm_peb_bmc (line 138)
[Np,Nx] = size(PEB.Ep);
The PEB.Ep contains one 4X1 double (the 1st item) and nine 4X1 sparse double.
May I know what's error it is?
Further, I got the output (see attached
) before the error. May I know where I can find the meaning of this picture?
Thank you very much!
Bests,
Delin
On Tue, Jan 10, 2017 at 6:55 AM, Zeidman, Peter <[log in to unmask]> wrote:
Dear Delin Sun
(1) Following SPM12 manual, I modeled two events: "neutral+emotional" and "emotional" in DCM's GLM. However, should I just model "neutral" and "emotional" seperately (as what I did in common fMRI GLM)?
Yes this is an ok way of modelling it - so neutral+emotional (which I'll call Task for convenience) will form the driving input to your DCM and emotional will form a modulatory input on one or more connections. Note that this isn't optimally efficient - the two conditions aren't independent and so there'll be covariance between the parameters - however it is the best you can do with this experimental design (ideally you would have at least 2 factors).
(2) Following SPM12 manual, I made the F contrast for "Effects of interest" (i.e. neutral+emotional) and two T contrasts (i.e. "neutral+emotional" and "emotional"). However, should I also get the T contrast for "emotional-neutral"?
These contrasts are simply there to help you choose ROIs (effects of interest has a special role in addition to this, for regressing out nuisance effects like motion). So the question of which contrasts you use depends on the hypotheses you wish to test. The "emotional" contrast will test for any effects of emotion trials relative to the unmodelled time in the GLM (i.e. the average of the neutral condition and inter-trial intervals). So that includes paying attention, viewing task instructions, etc. The "emotional-neutral" contrast will presumably control for things like viewing trial instructions / cues, so is probably a better contrast for identifying regions specifically modulated by emotion.
(3) I extracted time series from each ROI per session per subject. However, should I concatenate all sessions per subject before extracting time series, given that the number of blocks per session are unbalanced between emotional (1-2 blocks/session) and neutral (3-4 blocks/session) conditions?
In general I would concatenate sessions - instructions here https://en.wikibooks.org/wiki/SPM/Concatenation
(4) Following SPM12 manual, I defined the VOI of L amygdala by logically AND two images: [i1] the SPM results of "neutral+emotional" (height thresholded at p < 0.5) and [i2] the sphere covering the predefined ROI of L amygdala. On the other hand, I defined the VOI of vmPFC by logically AND two images: [i1] the SPM results of "emotional" (height thresholded at p < 0.5) and [i2] the sphere covering the predefined ROI of vmPFC. However, can I just extract signals from my predefined ROIs? Or, can I locate the VOI through using the SPM outputs corresponding to the contrast "emotional-neutral", which is more accurate than "neutral+emotional" or "emotional"?
Sorry I didn’t understand this question - perhaps you could rephrase it?
(5) I modified the "dcm_spm12_batch.m" in SPM official sample dataset of attention for my data. However, for the section of "Experimental inputs", I do not understand why there is a "33" in "DCM.U.u = [SPM.Sess(j).U(1).u(33:end,1) ..."
SPM adds 32 time bins before the first trial (for historical reasons). You should keep this as shown in the example script.
(6) for the section of "MODEL DEFINITION", may I know whether I am wrong to model as below:
DCM.a = [1 1; 1 1]; % the bi-directional intrinsic connections between L amygdala (ROI 1) and vmPFC (ROI 2)
DCM.b = zeros(2,2,2);
DCM.b(1,1,2) = 1;
DCM.b(2,2,2) = 1;
DCM.b(2,1,2) = 1; % task modulations on L amygdala, vmPFC and L amygdala --> vmPFC
DCM.c = [1 0;
0 1]; % input to both L amygdala and vmPFC
DCM.d = zeros(2,2,0); % There is no nonlinear modulations
The A and B matrices look correct, but the C-matrix does not look correct. I interpret your B-matrix as: emotion modulating the self-connection on region 1, the self-connection on region 2 and the connection from region 1 to region 2. The C-matrix is one row per region and one column per input. So your C-matrix says "Task is driving region 1 and emotion is driving region 2". Whereas, I think you meant to have Task driving both regions, which would be: [1 0; 1 0].
(7) For FFX method, is it correct that the best model is at least larger than the alternative models in 3 unit of log-evidence? On the other hand, I can't find log-evidence outputs in RFX method. How can I judge which model is the best using RFX method?
Correct. For the RFX, the convention is to look for a posterior probability or an exceedance probability of at least 95%.
(8) Am I right to extract the "mean DCM parameters per subject" from BMS.DCM.ffx.bma.mEps or BMS.DCM.rfx.bma.mEps for between-group comparison? If I understand it correctly, the DCM parameter, e.g. 0.123, of task modulation on region A-->B means that if there is a 1 unit neural activity increase in region A, there is a 0.123 unit neural activity increase in region B. Am I right? Is there any other DCM parameters or outputs for between-group comparison?
This procedure is fine, although you may wish to use the new Parametric Empirical Bayes (PEB) framework which is designed comparing groups of subjects. Temporary instructions at https://en.wikibooks.org/wiki/User:Peterz/sandbox .
Best
Peter
