Thanks for your help, Rick!
Could you help us understand what's happening at Deconvolve versus the following steps (below) when dealing with local white matter? Is 3dDeconvolve first taking out a global white matter signal then dealing with local white matter in a second step? And if so, why would we want the white matter regressed that way? Is there a way regress local white matter only?


# ================================ regress =================================

# compute de-meaned motion parameters (for use in regression)
1d_tool.py -infile Rest_mt.1D -set_nruns 1 \
-demean -write motion_demean.1D

# compute motion parameter derivatives (just to have)
1d_tool.py -infile Rest_mt.1D -set_nruns 1 \
-derivative -demean -write motion_deriv.1D

# create bandpass regressors (instead of using 3dBandpass, say)
1dBport -nodata 296 2.0 -band 0.01 0.1 -invert -nozero > bandpass_rall.1D

# create ROI regressor: FSWe
# (get each ROI average time series and remove resulting mean)
foreach run ( $runs )
3dmaskave -quiet -mask follow_ROI_FSWe+orig \
pb00.$subj.r$run.tcat+orig \
| 1d_tool.py -infile - -demean -write rm.ROI.FSWe.r$run.1D
end
# and catenate the demeaned ROI averages across runs
cat rm.ROI.FSWe.r*.1D > ROI.FSWe_rall.1D

# ------------------------------
# create ROI PC ort sets: FSvent

# create a time series dataset to run 3dpc on...

# detrend, so principal components are not affected
foreach run ( $runs )
3dDetrend -polort 4 -prefix rm.det_pcin_r$run \
pb00.$subj.r$run.tcat+orig
end

# catenate runs
3dTcat -prefix rm.det_pcin_rall rm.det_pcin_r*+orig.HEAD

# make ROI PCs : FSvent
3dpc -mask follow_ROI_FSvent+orig -pcsave 3 -prefix roi_pc_01_FSvent \
rm.det_pcin_rall+orig

# ------------------------------
# run the regression analysis
3dDeconvolve -input pb00.$subj.r*.tcat+orig.HEAD \
-ortvec bandpass_rall.1D bandpass \
-ortvec ROI.FSWe_rall.1D ROI.FSWe \
-ortvec roi_pc_01_FSvent_vec.1D ROI.PC.FSvent \
-polort 4 -float \
-num_stimts 6 \
-stim_file 1 motion_demean.1D'[0]' -stim_base 1 -stim_label 1 roll \
-stim_file 2 motion_demean.1D'[1]' -stim_base 2 -stim_label 2 pitch \
-stim_file 3 motion_demean.1D'[2]' -stim_base 3 -stim_label 3 yaw \
-stim_file 4 motion_demean.1D'[3]' -stim_base 4 -stim_label 4 dS \
-stim_file 5 motion_demean.1D'[4]' -stim_base 5 -stim_label 5 dL \
-stim_file 6 motion_demean.1D'[5]' -stim_base 6 -stim_label 6 dP \
-fout -tout -x1D X.xmat.1D -xjpeg X.jpg \
-fitts fitts.$subj \
-errts errts.${subj} \
-x1D_stop \
-bucket stats.$subj

# -- use 3dTproject to project out regression matrix --
3dTproject -polort 0 -input pb00.$subj.r*.tcat+orig.HEAD \
-ort X.xmat.1D -prefix errts.${subj}.tproject



# if 3dDeconvolve fails, terminate the script
if ( $status != 0 ) then
echo '---------------------------------------'
echo '** 3dDeconvolve error, failing...'
echo ' (consider the file 3dDeconvolve.err)'
exit
endif


# display any large pairwise correlations from the X-matrix
1d_tool.py -show_cormat_warnings -infile X.xmat.1D |& tee out.cormat_warn.txt

# create an all_runs dataset to match the fitts, errts, etc.
3dTcat -prefix all_runs.$subj pb00.$subj.r*.tcat+orig.HEAD

# --------------------------------------------------
# generate ANATICOR result: errts.$subj.anaticor+orig

# --------------------------------------------------
# ANATICOR: generate local FSWe time series averages
# create catenated volreg dataset
3dTcat -prefix rm.all_runs.volreg pb00.$subj.r*.tcat+orig.HEAD
3dLocalstat -stat mean -nbhd 'SPHERE(10)' -prefix Local_FSWe_rall \
-mask follow_ROI_FSWe+orig -use_nonmask \
rm.all_runs.volreg+orig

# -- use 3dTproject to project out regression matrix --
3dTproject -polort 0 -input pb00.$subj.r*.tcat+orig.HEAD \
-dsort Local_FSWe_rall+orig \
-ort X.xmat.1D -prefix errts.$subj.anaticor

# --------------------------------------------------
# create a temporal signal to noise ratio dataset
# signal: if 'scale' block, mean should be 100
# noise : compute standard deviation of errts
3dTstat -mean -prefix rm.signal.all all_runs.$subj+orig
3dTstat -stdev -prefix rm.noise.all errts.$subj.anaticor+orig
3dcalc -a rm.signal.all+orig \
-b rm.noise.all+orig \
-expr 'a/b' -prefix TSNR.$subj

# --------------------------------------------------------
# compute sum of non-baseline regressors from the X-matrix
# (use 1d_tool.py to get list of regressor colums)
set reg_cols = `1d_tool.py -infile X.xmat.1D -show_indices_interest`
3dTstat -sum -prefix sum_ideal.1D X.xmat.1D"[$reg_cols]"

# also, create a stimulus-only X-matrix, for easy review
1dcat X.xmat.1D"[$reg_cols]" > X.stim.xmat.1D