Anaticor: Task

Robin · November 13, 2018, 6:22pm

Just an observation.

Running task-analysis with anaticor can be manageable but also very very slow.

Using a typical pre-proc script (where the mask_epi_anat is applied to both 3dDeconvovle and REML) will vary in time depending on if you use band-pass or not.

Only using anaticor (also running 3dDeconvolve to get anaitcor free results):
real 30m52.777s

Not using anaticor but using bandpass (0.01 - 0.17):
real 67m49.350s

Using anaticor combined with bandpass:
real 436m21.949s

System:
SSD disks
Intel(R) Xeon(R) CPU E5-1650 v3 @ 3.50GHz
12 threads

The biggest time-consumer is the final loop of anaticor. I think it is this loop:
++ GLSQ loop:0123456789.0123456789.0123456789.0123456789.0123456789.
It only uses one thread.

Another interesting observation is that anaitcor makes expected results less signifciant.
In this example someone is stroking the right palm of the patient. Activating relevant motor areas on the left side of the brain. We do expect a similar but less significant activation on the opposite site but at the same threshold that is gone in anatcior. Anaticor also misses the activation in visual cortex. The q-scores are also different:
stats: q = 0.00087
reml: q = 0.0121

I’ll grant you that this is one single subject. And lowering the threshold for anaitcor is proabably still significant. The smoothenss might also be affected so that the 3dClustsim results might be different. But it’s to illustrate the difference between the two methods applied to exactly the same data. Any ideas?

Threshold: p < 0.0001 Cluster > 100
Top stats, bottom reml (no bandpass on either)

Robin · November 19, 2018, 11:03pm

Humble bump

rickr · November 26, 2018, 8:22pm

Hi Robin,

Band-passing in a task analysis will indeed be quite slow,
and not necessarily recommended. That can also give a
feel for how many regressors it is applying.

3dREMLfit does use multiple CPUs for the REML loop, but
just one for GLSQ (I had not been aware of that). Maybe
that is because the REML loop is much slower, so that
speed-up did not seem worthwhile.

Your numbers seem slow, but with many runs or high-res
data, they would be reasonable.

Consider going back to ponder why you are bandpassing.
While you would not omit it just for speed, it would
drastically reduce the processing time. There is also
a -regress_anaticor_fast option, which is indeed much
faster (maybe 50 times).

Is it REML that makes the (single subject, not group)
results less significant, or ANATICOR? For REML, it is
quite expected. For ANATICOR, it is not.

If ANATICOR alone changes the results much, then the
white matter mask or the anat/EPI registration might not
be good enough. If that is a concern, consider omitting
the option.

rick

Robin · November 27, 2018, 2:08pm

Hi Rick and thanks!!

We are using anaticor_fast. These numbers were from using the fast one, changed a long time ago :).

2). Wait a minute. What is the difference between REML and ANATICOR?
I have worked under the assumtion that everything with the REML extension (errts.REML and stats_REML) is the results you get WITH anaticor and the pure stats or errts are without. If we look at the afni.proc file:

regress-step: Produces the stats files via regression where the anaticor regressors are not used:


# run the regression analysis
3dDeconvolve -input pb03.$subj.r*.blur+tlrc.HEAD                              \
    -mask mask_epi_anat.$subj+tlrc                                            \
    -censor censor_${subj}_combined_2.1D                                      \
    -polort 4                                                                 \
    -local_times                                                              \
    -num_stimts 15                                                            \
    -stim_times_AM1 STIFILES GOHERE
      ...                 \                                                   \
    -stim_files FOR MOTION
      ....                                                     \                                                            \                                        \
    -fout -tout -x1D X.xmat.1D -xjpeg X.jpg                                   \
    -x1D_uncensored X.nocensor.xmat.1D                                        \
    -fitts fitts.$subj                                                        \
    -errts errts.${subj}                                                      \
    -bucket stats.$subj

The next step, (when using -regress_anaticor_fast) is to generate ANATICOR voxelwise regressors (anaticor regressors) using WM mask. This gives a REML-script.


# fast ANATICOR: generate local WMe time series averages
# create catenated volreg dataset
3dTcat -prefix rm.all_runs.volreg pb02.$subj.r*.volreg+tlrc.HEAD

# mask white matter before blurring
3dcalc -a rm.all_runs.volreg+tlrc -b mask_WMe_resam+tlrc                      \
       -expr "a*bool(b)" -datum float -prefix rm.all_runs.volreg.mask

# generate ANATICOR voxelwise regressors via blur
3dmerge -1blur_fwhm 30 -doall -prefix Local_WMe_rall                          \
    rm.all_runs.volreg.mask+tlrc

# -- execute the 3dREMLfit script, written by 3dDeconvolve --
# (include ANATICOR regressors via -dsort)
tcsh -x stats.REML_cmd -dsort Local_WMe_rall+tlrc -GOFORIT 8

Below is the REML-script thas ir run via tcsh -x stats.REML_cmd -dsort Local_WMe_rall+tlrc -GOFORIT 8


3dREMLfit -matrix X.xmat.1D -input pb03.sbj_02_A.r01.blur+tlrc.HEAD \
 -mask mask_epi_anat.sbj_02_A+tlrc -fout -tout \
 -Rbuck stats.sbj_02_A_REML -Rvar stats.sbj_02_A_REMLvar \
 -Rfitts fitts.sbj_02_A_REML -Rerrts errts.sbj_02_A_REML -verb $*

This gives me the output stats_sbj2_02_REML+tlrc.HEAD/BRIK.
So now I have:
stats_sbj2_02+tlrc.HEAD/BRIK (“normal” results, without anaticor regressors)
stats_sbj2_02_REML+tlrc.HEAD/BRIK (results given when also including the anaticor-regrressors)

So what do you mean when you say
"Is it REML that makes the (single subject, not group)
results less significant, or ANATICOR? For REML, it is
quite expected. For ANATICOR, it is not. "

The single subjects results are less significant in the stats_REML file than in the stats file. As I see it: ANATICOR vs No ANATICOR

Remember that this is task-based anaticor. We have had a similar discussion before and there we landed in that stats is just stats and stats_REML is the same X-matrix as the stats used but anaticor regressors added. For task.

Thanks!

rickr · November 27, 2018, 5:57pm

Hi Robin,

Thank you for clarifying your comparison.

ANATICOR and 3dREMLfit are 2 separate changes that affect
this single difference. ANATICOR is including local white
matter regressors (voxelwise). 3dREMLfit is applying an
ARMA(1,1) model in the context of a GLSQ one.

While ANATICOR should indeed be slower (in the context of
3dREMLfit), it should not be damaging to the stats (unless
there is some problem, such as poor selection of white matter
voxels).

But 3dREMLfit, as it accounts for temporal autocorrelation
in the noise, will strongly affect the single subject stats,
driving them downward. When you actually care about accurate
single subject t-stats, say, 3dREMLfit should be used.

But the beta maps should be reasonably similar between
3dREMLfit and 3dDeconvolve (plus any ANATICOR difference).
If you display the unthresholded beta maps between the 2
cases, how different are those?

At any rate, it is 3dREMLfit that will (appropriately)
drive down the t-stats, regardless of whether ANATICOR
is applied.

rick

Robin · November 27, 2018, 6:53pm

Thank!!

So, just to make sure I get this right.

I do get these datasets:

stats.sbj_A_+tlrc. = Normal stats file based on linear regression. Betas + t-stats where t-stats are not taking temporal auto correlation (important for singlesub)
errts.sbj_A+tlrc = The residuals (error-time-series), variance not accounted for in the model above

stats.sbj_A_REML+tlrc. = In this context it’s the same as stats but we have added ANATICOR voxe-wise regressors. The stats, however, should be different. “Lower” significance but “truer” since we accout for temporal autocorrelation.

errts.sbj_A_REML+tlrc = Should be similar to the error time-series from 3dDeconvolve (appart from ANATICOR differencies).

stats.sbj_A_REMLvar = Variance parameters (a,b, lam, stDev, -LogLik)

From my observations:

The beta-maps (stats and stats_REML betas) are very similar. If I threshold a bit on the betas the REML (anaticor-difference) dataset have a few clusters that the stats files does not have.
From before: The stats (t-scores) are not as significant in the stats_REML which is expected (since you so nicely explained it to me)
The errts and errts_REML graphs of the same region look very similar in the AFNI viewer.

So: You can simply see the stats and stats_REML as with and without anaticor (in this afni proc context) if you look at the betas, same for the errts and errts_REML. But when it comes to the t-stats they differ due to temporal autocorrelation being accounted for in REML. This means that it’s better to use the t-scores from REML when you use 3dMEMA since they are more correct.

This was very helpful! If I didn’t get it all wrong :).

This makes sense since we in another study found that ANATICOR (using betas from stats_REML in 3dMVM/3dttest++) gave “better” group maps. But in another, with large atrophy (alcoholics) REML (i.e. ANATICOR) killed some group differences. Which makes sense if the WM masks are not good due to atrophy.

Sounds about right?

If we don’t intend to use the t-scores (3dMEMA), but only the betas from stats_REML and errts_REML for smoothness estimates (for cluster corrections) is there a way to tell AFNI_proc not to run the time-consuming GLSQ part? The only info we give afni_proc is:


-regress_anaticor_fast \
  -regress_opts_reml \
        -GOFORIT 8 \
   -regress_reml_exec \

THANK YOU Rick!

rickr · November 28, 2018, 4:26pm

Hi Robin,

Your interpretation of the stats and errts datasets seems good.
And the comparisons of the results seems quite reasonable.

Yes, for 3dMEMA, it is appropriate to use the REML result
(though I really don’t think it would make much difference,
but still, REML is indeed considered to be more “correct”).

We will push more and more toward using REML in any case,
since there is little reason not to.

Yes, if you have the time, it might be good to verify the
WM masks in the case of ANATICOR. If they bleed into gray
matter (especially in areas that have BOLD responses of
interest), the results can get damaged, even worse than
using global signal regression. It would be much better
to not use ANATICOR than to use it with bad masks.

The GLSQ step cannot be separated from REML. It is the
step that actually generates the betas. The reason it is
so slow is because you are band passing (which I would be
inclined not to do in any case), which leads to very many
regressors.

rick

Robin · November 28, 2018, 5:35pm

That is perfect… We were trying to re-create some old results and tried bandpass (since we lost the original script). We use bandpass for resting state but not for task (any more).

How can I be sure I’m not using using global signal regression? It is across all voxel average signal right? As far as I can see we don’t use it.

We are, in a way, demeaning since we use polort A (usually 0, 1,3,4) that is including polort 0. But I guess that is not global signal regression. In a pretty normal afni-proc global signal regression is not used per default, right?

Thanks again Rick! No I’m comfortable using anaticor with task-data.

Best,
Robin

rickr · November 28, 2018, 5:53pm

No, it does not seem that you are trying to do GSR. It is
just that ANATICOR with a poor mask would be akin to it, but
likely even worse. The main point is that a WM mask going
into gray matter can be destructive, producing regressors of
no interest that contain the local gray matter signal.

Using polorts is quite different, it is just regressing slow
polynomial drifts. GSR means regressing the average time
series of brain voxels, which you are surely not trying to do.

rick