Dear AFNI gurus,
I’m trying to understand the permutation tests provided by 3dttest++, especially the ETAC method which is truly an ingenious solution to the cluster defining p trade-off problem. I have a question though, regarding the construction of the permuted, pure-noise dataset.
In the simplest permutation tests (comparing two-sample or one-sample mean), we permute or sign-flip the raw observations under the null hypothesis that there is no effect.
For more complex GLMs with covariates, one can permute the residuals of the reduced model (no effects of interest, the observations are explained only by the nuisance effects), as suggested in (Winkler et al., 2014).
However, the algorithm behind -ETAC and -Clustsim seems to permute the residuals of the full GLM model, according to the 2017 and 2019 paper. (I’m not sure whether I understood this part correctly.) This means, the permuted datasets will not be able to recreate the original observed data even in theory.
Would this somehow cause the reference distribution to contain fewer “extreme” samples and become narrower than it should be, and lead to an inflated p-value?
I did a simulation for the permutation test of one-sample mean using zero-mean Gaussian noise data. If I sign-flipped the raw data, the resulting p-values were uniformly distributed in [0,1] as expected. However, if I removed the mean of the sample (not zero due to random fluctuation) first, and permuted the residuals, the p distribution showed a peak below p=0.05. These results seem to support the above concern.
Could you make it more explicit about how the residuals are generated for the permutation procedure?
One rationale backing the current approach might be that we first need to regress out the (possibly existing) signals imposed by the driving input in order to generate data composed of pure noise. Is my understanding correct?
One rationale backing the current approach might be that we first need to regress out the (possibly existing) signals
imposed by the driving input in order to generate data composed of pure noise.
I’m not so familiar with the implementation, but I suspect that is likely the rationale. I wonder how to flip the sign with the observed data (rather than the residuals) when one or more covariates (e.g., age) are involved. Is your simulation process similar to what’s described in Dr. Cox’s following paper?
My guess is that the approach adopted in 3dttest++ is intended to address the issue of spatial neighborhood relatedness, which may cause some deviations in the case of a single voxel as you found out.
As the instigator/inventor/implementor of the ETAC approach, I have to say that I’ve never tested the accuracy of the tests for the false positive rate (FPR) for covariate estimates. It occurred to me to do so, and I even wrote scripts for that purpose. But didn’t carry those tests out.
What Anderson Winkler makes intuitive sense to me, but I never got that far in my implementation.
Personally speaking, I was very burned out when I got ETAC to the state it is in now. The whole clustering FPR subject is highly correlated in my mind with very unpleasant events. And now that I’m retired, I don’t plan to revisit it in any substantial fashion. I’m sorry if that leaves you hanging, but that’s the way it is.
The way the residuals are calculated for the case when there are covariates:
Before the analysis of a given voxel, the regression matrix for the covariate analysis is set up (for each of -setA and -setB) – call this matrix “X” – and pseudo-inverted – that is, matrix X has a column of all 1s (the “mean”) and then columns for all the covariates – 3dttest++.c function TT_matrix_setup.
In the analysis function regress_toz, the pseudo-inverse of X is applied to the values extracted from that voxel from each input dataset (in the corresponding -setA or -setB, of course) – this computes the (estimated) values of the covariates (“slopes”).
To get the residuals, the matrix X is multiplied back into the covariate estimates, giving the “fitted data vector”, which is then subtracted from the real data vector. Thus, the residuals are as “residual as they come” – they are what is left after ALL effects in the model are projected out.
I hope this helps, at least to reduce the confusion somewhat.
** Bob Cox