Replying years later in case anyone else is interested... I had the same question, wanting to understand precisely how t-statistics are computed when using 3dREMLfit. To manually replicate a t-statistic you'll need to construct or estimate the correlation matrix R as described in the (handwritten...) 3dREMLfit notes. This matrix is neither saved in output nor explicitly stored in memory by the solver given its size.

As a test example to understand construction, given sigma, coef, a, and lambda from stats.REML and stats.REMLvar results datasets, a design matrix X, and a vector or matrix C (tested here as a vector, for a single contrast, for a single run, with no censored TRs) which extracts the contrast of interest, I am able to replicate t-statistics with reasonable precision with a Matlab function explicitly forming R:

function T = REML_tStat(a, lam, std, coef, X, C)
dim = size(X,1);
r = zeros(dim,1);
r(1) = 1; r(2) = lam;
for i = 3:dim
r(i) = lam*a.^(i-2);
end
R = toeplitz(r);
T = coef/sqrt(std^2*C*inv(X'*inv(R)*X)*C');
end

I'm not aware of a convenient/straightforward way to generate whole volumes with AFNI functions but do not believe that would be a worthwhile endeavor anyway.

If there are any glaring errors here I would appreciate those being pointed out by AFNI experts.