HI Guillaume,
This bit of code seemed like a promising lead, until I realized that we are running this in SPM 5, not SPM 8. This code does not appear in SPM 5's version of the spm_spm.m file, so we are back to the drawing board....
Dan
On Jul 8, 2013, at 4:14 AM, Guillaume Flandin wrote:
> Hi Daniel,
>
> if sphericity is assumed and OLS estimates are used in SPM (SPM.xX.W is
> the identity matrix), then it is likely that the difference you observe
> is due to these lines in spm_spm.m:
>
> %-Modify ResMS (a form of shrinkage) to avoid problems of very low variance
> %----------------------------------------------------------------------
> try
> if ~strcmpi(spm_get_defaults('modality'),'fmri')
> ResMS = spm_read_vols(VResMS);
> ResMS = ResMS + 1e-3 * max(ResMS(isfinite(ResMS)));
> spm_write_vol(VResMS, ResMS);
> clear ResMS
> end
> end
>
> They refer to this article:
> http://www.ncbi.nlm.nih.gov/pubmed/22037420
> This is why I asked which modality are you running your estimation on as
> it will or not slightly rescale the residual sum of squares image.
>
> Best regards,
> Guillaume.
>
>
> On 06/07/13 01:20, Mathalon, Daniel wrote:
>> For non-sphericity assumptions, we assumed independent groups and equal variances. In this case, we do not find any indication that SPM is using ReML to estimate variances. It appears that the residual mean square values we calculate by applying the equations in the SPM code step by step do correspond to the residual mean square value provided by SPSS, but they do NOT correspond to the residual mean square value present for that voxel in the resMS image. So we think something else happens (e.g., some degrees of freedom adjustment) to alter the resMS image values beyond what we've been able to track in the code. Not sure what that is yet. More details will be sent by a member of my group further documenting the steps we've applied from the SPM script to calculate mean square residual value. Again, we seem to be missing a step that leads the resMS image to disagree with SPSS (and with our brute force calculation for a single voxel applying equations imbedded in the SPM
>
> code).
>>
>> Below are more details from Brian Roach from my lab.
>>
>> Dear SPMers,
>>
>> I checked some of our model parameters in the SPM.mat file, and I'm still not sure I understand the source of differences between the models. Everything that follows is based on SPM 5 code we're using to run this. Our SPM.xX.W (whitening matrix) is a sparse identity matrix so that when the spm_spm.m function whitens the data (lines 717-721):
>>
>> %-Whiten/Weight data and remove filter confounds
>> %--------------------------------------------------------------
>> fprintf('%s%30s',repmat(sprintf('\b'),1,30),'filtering')
>>
>> KWY = spm_filter(xX.K,W*Y);
>>
>> I think I should get values in KWY that are equal to my original data matrix Y (note: SPM.xX.K=1). Later, the betas and residual sum of squares are estimated (lines 732-739):
>>
>> %-General linear model: Weighted least squares estimation
>> %--------------------------------------------------------------
>> fprintf('%s%30s',repmat(sprintf('\b'),1,30),' estimation')
>>
>> beta = xX.pKX*KWY; %-Parameter estimates
>> res = spm_sp('r',xX.xKXs,KWY); %-Residuals
>> ResSS = sum(res.^2); %-Residual SSQ
>> clear KWY %-Clear to save memory
>>
>> In the belly of spm_sp.m, the calculation of the 'res' variable is given as (line 1403):
>>
>> Y = Y - sX.u(:,[1:r])*(sX.u(:,[1:r])'*Y); % warning('route1');
>>
>> If I repeat this calculation for my test voxel:
>>
>> spmR = KWY - (SPM.xX.xKXs.u*SPM.xX.xKXs.u'*KWY)
>>
>> I obtain residual values that are identical to my SPSS residuals or residuals returned by the regress.m function in matlab. This is in addition to obtaining identical beta values with all three methods. The only other thing I can see happening to this ResSS variable, is scaling (lines 821-824):
>>
>> %-Write ResSS into ResMS (variance) image scaled by tr(RV) above
>> %------------------------------------------------------------------
>> if length(Q), jj(Q) = CrResSS; end
>> VResMS = spm_write_plane(VResMS,jj,z);
>>
>> ...and lines 950-953:
>> %-Set VResMS scalefactor as 1/trRV (raw voxel data is ResSS)
>> %--------------------------------------------------------------------------
>> VResMS.pinfo(1) = 1/xX.trRV;
>> VResMS = spm_create_vol(VResMS);
>>
>> and based on that scaling variable (SPM.xX.trRV), I get a residual mean square value equal to that returned by SPSS or regress.m, but not equal to the actual value I see in the ResMS.img generated by running the spm job. What am I missing?
>>
>> thanks,
>> Brian
>> On Jul 5, 2013, at 5:02 AM, Guillaume Flandin wrote:
>>
>>> Dear Daniel,
>>>
>>> I have never used SPSS so I don't know what the options are but on the
>>> SPM side, which options did you specify for the non-sphericity
>>> assumptions (dependent/independent, equal/unequal variance)? Also, under
>>> which modality (fMRI/PET/EEG) were you running the model estimation?
>>>
>>> Best regards,
>>> Guillaume.
>>>
>>>
>>>
>>> On 04/07/13 20:12, Mathalon, Daniel wrote:
>>>> Dear SPMers,
>>>>
>>>> We are running a full factorial model in SPM 8 using a condition variable (representing two groups) and 8 additional covariates, including 1 Condition x Covariate interaction term.
>>>>
>>>> We've extracted the data from one voxel to compare results with SPSS linear regression or ANOVA programs, fully expecting them to agree (as a check that our script was working correctly).
>>>>
>>>> We find that the beta values estimated by SPM for that voxel for each variable in the model and the regression coefficients estimated by SPSS are in perfect agreement. The residual degrees of freedom in both SPM and SPSS also agree.
>>>>
>>>> However, we are finding that the t-test values generated in SPM and SPSS do not agree. For example, for the interaction term of interest, in SPSS t =1.459, but in SPM's t-map t =1.3414. We can't figure out why, but we suspect that the standard error of the beta must be estimated differently (since the betas and the degrees of freedom appear to be the same in SPSS and SPM).
>>>>
>>>> Does anyone have any ideas about why we are observing these discrepancies? Any insights or strategies for tracking this down would be greatly appreciated.
>>>>
>>>> Thanks,
>>>>
>>>> Dan
>>>>
>>>>
>>>> Daniel H. Mathalon, Ph.D., M.D.
>>>> Professor of Psychiatry
>>>> University of California, San Francisco
>>>>
>>>> Mail Address:
>>>> Psychiatry Service 116d
>>>> San Francisco VA Medical Center
>>>> 4150 Clement St.
>>>> San Francisco, CA 94121
>>>>
>>>> Office phone: (415) 221-4810, ext. 3860
>>>> Fax: (415) 750 6622
>>>> e-mail: [log in to unmask]
>>>>
>>>
>>> --
>>> Guillaume Flandin, PhD
>>> Wellcome Trust Centre for Neuroimaging
>>> University College London
>>> 12 Queen Square
>>> London WC1N 3BG
>>>
>>
>>
>>
>
> --
> Guillaume Flandin, PhD
> Wellcome Trust Centre for Neuroimaging
> University College London
> 12 Queen Square
> London WC1N 3BG
>
|