In principle, one can replace the data (FMRI amplitudes) with their ranks, and then do ANCOVA -- this is suggested in the paper
Analysis of Covariance Using the Rank Transformation,
W. J. Conover and Ronald L. Iman
Biometrics, Vol. 38, No. 3, Special Issue: Analysis of Covariance (Sep., 1982), pp. 715-724
How well this will work with FMRI data is an interesting question. I think you could try this using '3dTsort -rank' if you first created a single dataset comprising the desired betas for analysis in the various sub-bricks. Then use 3dLME to carry out the analysis. Maybe it would work.