Dear Takanori,
>I would like to ask you about a practical computation in spm_spm.
>I have been worried about whether the following is bug or not.
>My question is about the different calculation of UFp between spm99 and spm2.
>I got the different trRV used for UFp calculation in each version of spm.
>
>In spm 99, we can find the following expressions in spm_spm.m around l:443:
>%--------------------------------------------------------------
>[trRV trRVRV] = spm_SpUtil('trRV',xX.xKXs,V);
>[trMV trMVMV] = spm_SpUtil('trMV',spm_FcUtil('X1o',xCon(1),xX.xKXs),V);
>...
>tmp = (sum((h*beta).^2,1)/trMV) > UF*ResSS/trRV;
>%--------------------------------------------------------------
>Here, xX.xVi.Vi is i.i.d and
>KVi = spm_filter('apply',xX.K, xX.xVi.Vi); $B!! (B% l:439
>V = spm_filter('apply',xX.K,KVi');. % l:440
>
>On the other hand, spm_spm.m of spm 2 around l:445 says following:
>%--------------------------------------------------------------
>trRV = spm_SpUtil('trRV',xX.xKXs);
>X1o = spm_FcUtil('X1o', xCon(1),xX.xKXs);
>trMV = spm_SpUtil('trMV',X1o);
>....
>j = sum((Hsqr*beta).^2,1)/trMV > UF*ResSS/trRV;
>%--------------------------------------------------------------
>
>The main difference between them seems to be the use of spm_FcUtil.
>In spm99, trRV is calculated from filtered design matrix: xKXs and
"filtered" i.i.d matrix:KVK,
>and hence, trRV have a value of {total scans - rank(X) - rank(HPF)}
>if I do not use low pass filter.
>However, in spm2, trRV is calculated from only filtered design matrix.
>As a result, trRV have a value of {total scans - rank(X)}.
>I.e., the reduction of df due to high pass filtering is not taken into
account in spm2.
>
>For example, when data has 120 time points and the design matrix has 1
box-car regressor,
>1 constant term, and high pass filter with cut off =120sec containing 6th
order DCT basis functions,
>we get trRV = 112.000 for spm99, on the other hand, trRV = 118 in spm2.
>
>I am worried about this difference because UFp would affect the globally
pooled autocorrelation estimation.
>Which is correct? Or why are the expressions changed?
>I think that spm 99 procedure is rather correct
>because ResSS, the denominator of F-threshold inequality,
>is computed using the filtered Y and X.
Your analysis is absolutely correct. For computational expediency, the
'effects interest'
F-statistics threshold is computed under i.i.d. assumptions AFTER filtering
in SPM2.
Because the filter matrix is a residual forming matrix this assumption is,
strictly
speaking, not valid.
However, it does not really matter because, in SPM2 the threshold UFp is
only used
to estimate the serial correlations in a first pass, if V has not been
specified. In the
second pass, UFp is not computed and the appropriate degrees of freedom are
calculated correctly as in SPM99
i.e. (after the loop over voxels).
xX.V = spm_filter(xX.K,spm_filter(xX.K,W*V*W')'); % KWVW'K'
[trRV trRVRV] = spm_SpUtil('trRV',xX.xKXs,xX.V); % trRV (for X)
In principle this could affect the computation of V but the effect would be
equivalent to
using a slightly different UFp, which is arbitrary anyway.
I hope this explains things - Karl
|