Dear Carsten,
Your analysis is absolutely correct. The covariance components returned by spm_Ce correspond to:
% C{i} = AR(a) - a*dAR(a)/da;
% C{i + 1} = AR(a) + a*dAR(a)/da;
with a = 0.2. In other words, they are the first order Taylor approximations to an AR process with a = 0 and a = 0.4 (using an expansion point of a = 0.2). As you note, these correspond to a white noise process and an AR(1) with a = 0.4. This is a fairly arbitrary parameterisation but is sufficient for most fMRI applications.
With very best wishes – Karl
-----Original Message-----
From: Carsten Allefeld [mailto:[log in to unmask]]
Sent: 17 October 2016 14:45
To: Friston, Karl <[log in to unmask]>
Subject: Fwd: estimation of serial correlations
Dear Karl,
have you seen my email below? I sent it to you and the SPM list. If you can find the time, I would very much appreciate some insight.
Thank you very much!
Carsten
----- Forwarded Message -----
Dear Karl,
thank you for your answer.
> The model used for serial correlations in standard MRI timeseries is
> usually described as "an AR model plus white noise". Practically, this
> means there are two covariance components (C{1} and C{2} in your code
> below), where the first is the identity or white component and the
> second is an AR component.
> This is simply a covariance component based upon an AR(1) process with
> an AR coefficient of 0.2 (or user-specified).
I understand, but as far as I can tell, SPM12 doesn't seem to be actually doing this.
According to your description (consistent with the paper), the components should be the identity matrix (truncated to the first 6 x 6 entries)
>> eye(6)
ans =
1 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0
0 0 0 1 0 0
0 0 0 0 1 0
0 0 0 0 0 1
and the autocorrelation matrix of an AR(1) process with parameter 0.2
>> toeplitz(0.2 .^ (0 : 5))
ans =
1.0000 0.2000 0.0400 0.0080 0.0016 0.0003
0.2000 1.0000 0.2000 0.0400 0.0080 0.0016
0.0400 0.2000 1.0000 0.2000 0.0400 0.0080
0.0080 0.0400 0.2000 1.0000 0.2000 0.0400
0.0016 0.0080 0.0400 0.2000 1.0000 0.2000
0.0003 0.0016 0.0080 0.0400 0.2000 1.0000
But when I set up a first-level model with standard parameters, including modeling of serial correlations, these are the resulting covariance components:
>> full(SPM.xVi.Vi{1}(1:6, 1:6))
ans =
1.0000 -0.0000 -0.0401 -0.0160 -0.0048 -0.0013
-0.0000 1.0000 -0.0000 -0.0401 -0.0160 -0.0048
-0.0401 -0.0000 1.0000 -0.0000 -0.0401 -0.0160
-0.0160 -0.0401 -0.0000 1.0000 -0.0000 -0.0401
-0.0048 -0.0160 -0.0401 -0.0000 1.0000 -0.0000
-0.0013 -0.0048 -0.0160 -0.0401 -0.0000 1.0000
>> full(SPM.xVi.Vi{2}(1:6, 1:6))
1.0000 0.4000 0.1201 0.0320 0.0080 0.0019
0.4000 1.0000 0.4000 0.1201 0.0320 0.0080
0.1201 0.4000 1.0000 0.4000 0.1201 0.0320
0.0320 0.1201 0.4000 1.0000 0.4000 0.1201
0.0080 0.0320 0.1201 0.4000 1.0000 0.4000
0.0019 0.0080 0.0320 0.1201 0.4000 1.0000
It turns out that the mean of these two matrices is actually the AR(1) matrix:
>> full(SPM.xVi.Vi{1}(1:6, 1:6) + SPM.xVi.Vi{2}(1:6, 1:6)) / 2
ans =
1.0000 0.2000 0.0400 0.0080 0.0016 0.0003
0.2000 1.0000 0.2000 0.0400 0.0080 0.0016
0.0400 0.2000 1.0000 0.2000 0.0400 0.0080
0.0080 0.0400 0.2000 1.0000 0.2000 0.0400
0.0016 0.0080 0.0400 0.2000 1.0000 0.2000
0.0003 0.0016 0.0080 0.0400 0.2000 1.0000
but the complement
>> full(SPM.xVi.Vi{2}(1:6, 1:6) - SPM.xVi.Vi{1}(1:6, 1:6)) / 2
ans =
0 0.2000 0.0801 0.0240 0.0064 0.0016
0.2000 0 0.2000 0.0801 0.0240 0.0064
0.0801 0.2000 0 0.2000 0.0801 0.0240
0.0240 0.0801 0.2000 0 0.2000 0.0801
0.0064 0.0240 0.0801 0.2000 0 0.2000
0.0016 0.0064 0.0240 0.0801 0.2000 0
is far from being the identity matrix. According to the comments in spm_Ce, this is the derivative of the AR(1) matrix with respect to a at a = 0.2:
% C{i} = AR(a) - a*dAR(a)/da;
% C{i + 1} = AR(a) + a*dAR(a)/da;
Could you shed light on this behavior?
Thank you,
Carsten
|