Aaron/Roberto/Stephen/Donald/Alexa/others -
I finally found some time to simulate some of the recent discussion on the SPM list about pooled vs partitioned errors in repeated-measures ANOVAs (rmANOVAs). Attached is a Zip file containing a matlab script ("check_pooled_error.m") and some JPEGs produced by that script, that show how, provided one can estimate the nonsphericity (correlation in the error/non-white error/non-IID error) accurately, then a model with a single pooled error term can be more powerful in testing true effects, yet still valid in rejecting null effects, than separate models that partition the error for each ANOVA effect. So these two issues - error partitioning and nonsphericity - are inter-related (as in email below).
Estimating nonsphericity is difficult with only one sample of data - one can get very imprecise estimates with few subjects relative to number of conditions. Partitioning the error reduces the degree of nonsphericity possible (indeed, completely removing any concern about nonsphericity when there is only one df in the contrast, ie numerator of the F-statistic for effect). This is probably why standard stats packages advise partitioning the error in rmANOVAs (and why the common corrections for nonsphericity, eg Greenhouse-Geisser, are often conservative).
In neuroimaging data however, we have multiple samples - one per voxel. If one assumes that the nonsphericity is common across those voxels (or a subset that show any sign of experimental effects; see Friston et al, 2002, Neuroimage), then one can pool the data across voxels to get a much more precise estimate of the error covariance. This nonsphericity estimate can then be inverted, to "prewhiten" the data, which recovers maximal efficiency, ie power to detect effects. (Note this sense of "pooling data across voxels" is different from "pooling error across conditions", though the issues are related, which is the main point of this email). Whether the correlation in the error is identical across (active) voxels is, of course, an assumption that might be questioned...
Anyway, the attached figures and matlab function should illustrate these points. At the top of the matlab script are instructions for running the 5 "demos" that produce the 5 graphs, together with their take-home messages. I cut and paste them here:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Demonstration 1: If the data are generated by a single IID error
% (albeit unlikely), then pooling is more efficient than partitioning error
%
% Run by:
%
% S.ftit = 'Single white error'; % Just title for plot
% [Pool,Part,QV] = check_pooled_error(S);
%
% ie using defaults below, ie a 2x2 ANOVA design with factors A and B
%
% Results:
%
% 1. When there is no true effect (B=0 for main effect of B; red lines),
% approx 5 percent of voxels survive p<.05 in both types of model - ie
% type I error controlled with both pooled and partitioned error
%
% 2. When there is an effect however (eg, big main effect of A, or smaller
% interaction), a greater percentage of voxels survive p<.05 for pooled
% than partitioned error model (ie more power, because greater df's with
% which to estimate error)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Demonstration 2: If the data are generated by partitioned but IID error
% terms, such that fixed correlation between errors induced, then pooling
% is still more efficient, with no bias (unless those error terms are
% non-IID, eg different (nonspherical) errors, as in demos below).
%
% Run by:
%
% S.GenPart = 1;
% S.ftit = 'Partitioned errors';
% [Pool,Part,QV] = check_pooled_error(S);
%
% Results:
%
% 1. Similar to above, noting that when there is no true effect
% (B=0 for main effect of B; red lines), pooled model still unbiased
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Demonstration 3: If there is greater correlation in the error (nonsphericity),
% owing to similar scores across conditions from the same subject (which is
% likely), then pooled error becomes biased, but partitioning error is less
% biased (indeed, UNbiased when only one numerator df per effect, as here)
%
% Run by:
%
% S.GenPart = 0;
% S.V = [1 0 0.5 0; 0 1 0 0; 0.5 0 1 0; 0 0 0 1];
% S.ftit = 'Correlated error';
% % Eg because only true effects of A and AxB, make these correlated
% [Pool,Part,AQV] = check_pooled_error(S);
%
% Results:
%
% 1. Now note that >5% of voxels have p<.05 under pooled error model (thin
% red line) - ie invalid - but still correct 5% for partitioned error model
% (thick red line) - ie still valid.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Demonstration 4: If you try to estimate error covariance for each voxel
% separately (using ReML), that estimate is poor (unless large number of
% subjects), resulting in "poor" p-values.
%
% Run by:
%
% S.GenPart = 0;
% S.PreWhiten = 1; % Separate ReML for each voxel
% S.V = [1 0 0.5 0; 0 1 0 0; 0.5 0 1 0; 0 0 0 1];
% S.ftit = 'Voxel-wise estimation of correlated error';
% [Pool,Part,AQV] = check_pooled_error(S);
%
% Results:
%
% 1. Warnings as ReML unable to estimate error covariance for some random
% data (some voxels), particularly when few subjects
%
% 2. Note that >5% of voxels have p<.05 under pooled error model (thin
% red line) - ie invalid.
%
% 3. Note that power to detect real effects also reduced for pooled versus
% partitioned error (green and blue lines lowered).
%
% This danger is probably why partitioning error is preferred in
% behavioural statistics when typically only one variate (one voxel) - any
% method (eg Greenhouse-Geisser) for correcting for correlated errors,
% based on one sample of data, will entail an inefficient estimate of that
% correlation (with few df's), so a poor correction (eg too conservative
% in case of G-G).
%
% (This is the same general problem that estimating covariance matrices
% efficiently requires large numbers of observations!)
%
% In imaging data however, we have many samples (voxels), so if we assume
% that error correlation across conditions is same across voxels, we can
% pool the samples to get a more precise estimate of that error
% correlation (ie effectively increase number of observations so as to
% increase accuracy of estimate of error covariance).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Demonstration 5: If you try to estimate the error covariance by pooling
% data across voxels, in a two-step procedure like that used in SPM, then
% you recover greater efficiency of pooled than partitioned model (when
% there is a true effect), together with valid results when there is no
% effect. **Note however that this is only true if the underlying error
% correlation IS the same across voxels (as in this code)...**
%
% Run by:
%
% S.GenPart = 0;
% S.PreWhiten = 2; % One ReML from data pooled across voxels
% S.V = [1 0 0.5 0; 0 1 0 0; 0.5 0 1 0; 0 0 0 1];
% S.ftit = 'Voxel-wide estimation of correlated error';
% [Pool,Part,AQV] = check_pooled_error(S);
%
% Results:
%
% 1. We recover the situation in Demonstration 1, where pooling error is
% more powerful when effects exist, but still valid when they do not.
%
% 2. [ To visualise nonsphericity: figure,imagesc(AQV),colormap('gray') ]
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Again, I am sure a proper statistician will explain more clearly / correctly (eg, which much more elegant code), but hope this code helps a bit, at least for those of us with a more intermediate level of understanding.
BW,R
-------------------------------------------------------
Dr Richard Henson
Assistant Director, Neuroimaging
MRC Cognition & Brain Sciences Unit
15 Chaucer Road
Cambridge, CB2 7EF, UK
Office: +44 (0)1223 355 294 x522
Mob: +44 (0)794 1377 345
Fax: +44 (0)1223 359 062
http://www.mrc-cbu.cam.ac.uk/people/rik.henson/personal
-------------------------------------------------------
From: Aaron Schultz [mailto:[log in to unmask]]
Sent: 10 January 2011 02:42
To: Rik Henson
Cc: [log in to unmask]
Subject: Re: [SPM] Design matrix for each or all subject(s)?
On Jan 9, 2011, at 3:03 PM, Rik Henson wrote:
Aaron -
A very clear exposition - am sure it will be very helpful to many people - so thanks. Just one minor thought:
I don't know for certain, but It does make sense to me that the partitioned and pooled variances would collapse onto one another in the special case where the Subject by within interactions are non-significant, but this will rarely be the case. Not only is the partitioned approach technically correct, but it will tend to yield higher t and F values since the partitioned error terms will almost always be less than the pooled error (after all the partitioned terms are just portions of the pooled error; in both cases the SS across all terms and errors must sum to the total SS of Y).
Though, as you say, the SS for each error partition must be less than the SS for the pooled error, the MS can be more (or less), since it is the SS divided by the df's for each error term - so surely an effect will not always have higher t/F values with a partitioned error?
But this raises an interesting point: if the first of two error partitions (say) has a higher MS than the other one, such that the pooled MS is less than that of the first error partition (and then presumably must be greater than that of the second), then by definition we have an inhomogeneous error, so should partition. In other words, in situations where the statistics for an effect using a pooled error are markedly better than using a partitioned error, then we should be worried about inhomogeneous error and use a partitioned error instead! (I guess the counter-argument is that we are only estimating the true error terms from a sample, and one might get unusually high or low MS for particular error terms simply due to poor estimation, eg too few dfs, bringing us back to the heart of the problem!)
In any case, my experience (with multifactorial but fully within-subject ANOVAs) is that the uncorrected statistics are fairly similar, but corrected statistics are often better with a pooled error, because RFT is less conservative (with small numbers of subjects). But perhaps my experiences are atypical (or I have misremembered).
I believe the partitioned approach to be the more general and more correct way to perform ANOVAs, and I would not trust any results where the full model residual error does not produce the same error term as the partitioned approach. My reasoning for this is that this is how ANOVA is presented in text books and how statistical packages perform and report ANOVA results.
Yes, that is my experience of standard stats packages too - and it is probably safer to partition (as in my previous email) - with an additional reason for this being that partitioning the error also entails simpler error covariance matrices, ie smaller nonsphericity problems. With mass-univariate approaches like SPM, however, we have the luxury of pooling across voxels to get more precise estimates of nonsphericity (assuming it is similar across voxels), which potentially allows more complex error covariance matrices (that arise with single, pooled error terms).
Feel free to forward to the list if you think this is interesting enough (others might be getting sick of the debate! (if I ever get enough time, I might try to resolve with the basic math / some simulations with small N!))
Best wishes
Rik
