Print

---

Dear Hanneke,

I am copying this to the SPM helpline in case other people using DCM10 or DCM12 are encountering similar issues. I spoke to Klaas last week and he told me that several people were concerned about the differences between the results of model inversion using DCM8 and DCM10. It is probably useful to comment on the issue generally and then focus on the example you sent me.

The changes from DCM8 two DCM10 were driven by advances in two-state and stochastic dynamic causal modelling. These extensions and generalisations called for a revision and simplification of the prior densities (and parameterisation) of the models for fMRI used in DCM8. The idea was to make all DCMs for fMRI as internally compatible (and simple) as possible. I think it is important to emphasise that DCM is actually getting simpler with every revision and, in principle, easier to understand and tailor for specific applications. Having said this, the simplifications meant that with relatively noisy data and difficult inversions, DCM10 occasionally ‘flat-lined’. In other words, it supposed that the most likely explanation for the data was that they were caused entirely by observation noise. This led to the revision that Guillaume e-mailed about shortly after Christmas. This revision changed the priors – not on the parameters – but on the precision of observation noise – essentially telling DCM that there was signal to the explained. I will refer to this as DCM12.

In short, the continued development of DCM has entailed changes to the priors of both parameters and precisions and has, consequently, led to changes in posterior estimates and model selection. In one sense, this is a natural part of DCM development (and can be seen as informal model optimisation). On the other hand, I appreciate that this can be frustrating for people trying to publish reproducible results. As a general comment, the philosophy behind DCM is that it forces people to think mechanistically about how their data were caused – in terms of distributed neuronal processing. One could argue that the same issue is now emerging in terms of the prior beliefs about models in general and that people should be encouraged to try different priors – noting that every time you change the priors you implicitly change the model. Having said this, this is probably not terribly useful for practitioners of DCM, who would like to use canonical priors that work in most situations. Klaas and others are currently thinking about how to deal with this. Our internal experience here is that the current priors used in DCM12 provide more sensitive model (and parameter) inference than any previous version (the same appears to be true for the examples you sent me).

One thing that might help is to have some way of reassuring oneself that model inversion gives sensible results. Below, I enclose a short Matlab routine (that is compatible with early DCM.mat files), which provides some simple diagnostics and rules of thumb to ensure everything is okay (please see the help below for a fuller description and the attached PDF for an exemplar output).


In relation to the specific example you sent; I'll briefly describe history for people who are interested: you inverted the same dataset using DCM8, two versions of DCM10 and DCM12. As noted above, the two versions of DCM10 provided null results, with parameter estimates that were close their prior expectations, while DCM8 and DCM12 gave more interesting (if different) results. First, it may not be surprising that DCM10 flat-lined because this particular DCM is quite ambitious – the regions are very close together, which means that the data are highly correlated (the smallest correlation between regional responses was 0.84). Furthermore, you are allowed for full connectivity in both the A and B matrices. Finally, you used event related modulatory (bilinear) effects that were a subset of the driving events. Having said this, there are clearly interesting things going on as evidenced by the DCM8 and 12 inversions.

Looking at the results, they actually appear fairly consistent (although the bilinear effects focus on afferents to the first region in DCM12 and the second in DCM8). The posterior confidence and number of significant effects is greater with DCM12 than DCM8. The slight caveat here is the interesting self inhibition on the first region that is estimated to be positive (self excitation). It may well be the case that this region’s sensitivity becomes self excitatory during the transient appearance of the stimulus and you should believe the posterior estimates. Alternatively, your prior beliefs that regions should have self inhibition may be sufficiently strong to ignore this model. I increased the prior shrinkage variance on self connections (to 1/256) and fit the model – giving the results in the attached PDF. This had a free energy that was significantly better than either the original DCM8 or DCM12 inversions.

In summary, the fact that you were surprised by the results tells you immediately that they have violated your prior beliefs and that the specification of these beliefs in spm_dcm_fmri_priors was not sufficiently precise. When these prior beliefs were implemented the results seemed sensible and better than previous results. The outstanding problem is now how to operationalise this sort of model search and make it easy for people who are less experienced with DCM. I would value your thoughts and I hope that Klaas will comment upon this?

With very best wishes – Karl



function [DCM] = spm_dcm_fmri_check(P)

% post-hoc diagnostics for DCM (bilinear or nonlinear) of fMRI data

% FORMAT [DCM] = spm_dcm_fmri_check(DM)

% DCM - DCM structure or its filename

%

% This routine provides some diagnostics to ensure model inversion has

% converged. It plots the predicted and observed responses over all regions

% and provides the coefficient of determination – or percent variance

% explained. This should normally be above 10%. An abnormally low

% coefficient of determination is highlighted in red. Quantitatively, one

% would normally expect to see one or more extrinsic (between source)

% connections with the strength of 1/8 Hz or greater. If all the extrinsic

% posterior expectations are below this value, then this suggests a failure

% of convergence or that the data are very noisy (possibly due to using

% very small regions of interest to summarise regional responses). Finally,

% the posterior correlations among all parameters are shown in terms of a

% correlation matrix. The number of effective parameters estimated is

% reported in terms of the (KL) divergence between the posterior and

% prior densities over parameters. This is divided by the log of the

% number of observations, by appealing to the Bayesian information

% criterion. The divergence corresponds to complexity or Bayesian

% surprise. Normally, one would expect the posterior and prior to diverge

% in a non-trivial fashion.

%

% Posterior densities are shown as bars with 90% confidence intervals in

% pink. An informed model inversion would normally provide posterior

% densities with confidence intervals that are, for some connections,

% displaced from prior expectations (at or around zero).

%

% This routine is compatible with DCM8, DCM10 and DCM12 files.

%__________________________________________________________________________

% Copyright (C) 20012 Wellcome Trust Centre for Neuroimaging

% Karl Friston



% $Id: spm_dcm_fmri_check.m karl $



%-Load DCM structure

%--------------------------------------------------------------------------

if ~nargin

  uiopen('load');

elseif isstruct(P)

  DCM = P;

else

  load(P)

end



% Assemble diagnostics

%==========================================================================

% coefficient of determination (percent variance explained)

%--------------------------------------------------------------------------

PSS  = sum(sum(DCM.y.^2));

RSS  = sum(sum(DCM.R.^2));

D(1) = 100*PSS/(PSS + RSS);



% largest absolute posterior expectation (extrinsic connections)

%--------------------------------------------------------------------------

try

   A = DCM.Ep.A;

catch

   A = DCM.A;

end

D(2) = max(max(abs(A - diag(diag(A)))));



% complexity and effective number of parameters estimated

%--------------------------------------------------------------------------

qE = spm_vec(DCM.Ep);

pE = spm_vec(DCM.M.pE);

qC = DCM.Cp;

pC = DCM.M.pC;

k  = rank(full(pC));

pC = spm_inv(pC);



D(3) = (trace(pC*qC) + (pE - qE)'*pC*(pE - qE) - spm_logdet(qC*pC) - k)/2;

D(3) = D(3)/log(DCM.v);



% Plot summary of inversion

%==========================================================================

spm_figure('GetWin','DCM diagnostics'); clf



% plot predicted and observed regional responses

%--------------------------------------------------------------------------

subplot(2,1,1);

t = (1:DCM.v)*DCM.Y.dt;



plot(t,DCM.y,t,DCM.y + DCM.R,':');

str = sprintf('variance explained %0.0f%%', D(1));

str = {'responses and predictions',str};

if D(1) > 10

  title(str,'FontSize',16);

else

  title(str,'FontSize',16,'Color','r');

end

xlabel('time {seconds}');



% posterior densities over A parameters

%--------------------------------------------------------------------------

try

  i = spm_fieldindices(DCM.Ep,'A');

catch

  i = 1 + (1:DCM.n^2);

end

qE = spm_vec(DCM.Ep);

qC = DCM.Cp;

subplot(2,2,3)

spm_plot_ci(qE(i),qC(i,i)), hold on

str = sprintf('largest connection strength %0.2f', D(2));

str = {'intrinsic and extrinsic connections',str};

if D(2) > 1/8

  title(str,'FontSize',16);

else

  title(str,'FontSize',16,'Color','r');

end

xlabel('parameter}');

axis square



% posterior correlations among all parameters

%--------------------------------------------------------------------------

subplot(2,2,4)

imagesc(spm_cov2corr(DCM.Cp))

title('posterior correlations','FontSize',16)

str = sprintf('estimable parameters %0.0f', D(3));

str = {'posterior correlations',str};

if D(3) > 1

  title(str,'FontSize',16);

else

  title(str,'FontSize',16,'Color','r');

end

axis square