Spring 2007 Changes to 3dDeconvolve

RW Cox — SSCC / NIMH / NIH / DHHS / USA / Earth

Click here to jump down to the description of Amplitude Modulated FMRI regression analysis.

Small changes: to the defaults, new options, etc.
- -nobout and -full_first are now the defaults. These changes mean that if you want the β weights for the baseline parameters in the output -bucket dataset, you have to specify -bout on the command line. If you want the full-model statistics to appear last in the dataset, you have to specify -nofull_first on the command line.
- Even if you do not give the -fout option on the command line (indicating you do not want F-statistics for various hypotheses to be calculated), the program will still compute the full model F-statistics. If you don't want that for some reason, you have to use the new -nofullf_atall option.
- If you do not give a -bucket option on the command line, then the program will act as if you had given -bucket Decon. (This is known as the "Ah need a bucket" change, with apologies to KFC.)
- The program now always outputs (to a file) the regression matrix X, even if you don't give a -x1D option. The default filename will be the same as the -bucket prefix, with the suffix .x1D added.
  - The matrix file format has been slightly altered to store column labels in XML-style comments in the header. (Previously, the matrix was just written out as an array of unlabeled numbers.) These labels will be useful in an upcoming regression matrix analysis program being planned by Ziad Saad. They are also useful in the new program 3dSynthesize (cf. infra).
- 3dDeconvolve used to fail with the -nodata option combined with -stim_times. This crash should be a thing of the past.
  - When using -nodata, the program needs to know the length of the (non-existent) imaging data (number of TRs) and it also needs to know the TR. The simplest and best way to specify these values is to put them immediately after the -nodata option; for example -nodata 300 2.5 to indicate 300 time points with TR=2.5 s.
  - If you don't do the above, then if you use -nlast, that value (+1) will be used as the number of TRs. If you don't give the TR in some way, then the default -nodata TR is 1.0 s. This TR is unimportant if you only use -stim_file, but is crucial if you use -stim_times with -nodata or with -input1D.
- New option -float (or -datum float) can be used to make all the output datasets be stored in floating point format. In the past, only scaled shorts were possible, and the limited (16-bit) precision of these sometimes caused problems. Shorts are still the default, but at some point in the future I may change the default to floats — if/when this happens, the option -short can be used if you like the more compact format.
- The program now reports when -stim_times time values are out of the time span of the dataset. These are not fatal errors, but can help notify you to potential problems of your timing files. (This problem is known as the PSFB syndrome — it's not as bad as the Mike Beauchamp syndrome, but try to avoid it.)
- The labels for the -bucket output dataset sub-bricks have been changed slightly to be more consistent and readable (e.g., Tstat instead of t-st to indicate a t-statistic).
- 3dDeconvolve now computes a recommended -polort value (1 degree for every 150 s of continuous imaging). If your input value is less than this, a non-fatal WARNING message is printed. If you use -polort A, then the program will automatically choose the polynomial degree to use for detrending (AKA high pass filtering).
- A new CSPLIN() model for -stim_times is now available. This function is a drop-in replacement for TENT(), with the same 3 arguments. The basis functions are cardinal cubic splines, rather than cardinal linear splines. CSPLIN() will produce smoother looking HRF curves, if plotted with -TR_times less than the dataset TR. (As always, if you are going to change your analysis methodology, run some data the old way and the new way, then compare the results to make sure you understand what is happening!)
Detection and processing of matrix errors; new -GOFORIT option
- 3dDeconvolve now makes several more checks for "bad things" in the regression matrix.
  - Besides checking the full matrix condition number, it also checks several sub-matrices: the signal sub-model, the baseline sub-model, the -polort sub-model, and the -stim_base sub-model.
  - Each check is printed out and labeled as to how good the program "thinks" it is. Potentially bad values are flagged with ** BEWARE **.
  - N.B.: 3dDeconvolve's condition number is not exactly the same as that computed by Matlab. 3dDeconvolve first scales the matrix columns to have L²-norm = 1, and then computes the condition number from the ratio of the extreme singular values of that matrix. This method prevents the pathology of saying that the matrix diag(1,10^–6) is ill-conditioned.
  - Other "bad things" that the program checks for include duplicate stimulus filenames, duplicate regression matrix columns, and all zero matrix columns.
- If "bad things" are detected in the matrix (each will be flagged in the text printout with a warning message containing the symbols '!!'), then 3dDeconvolve will not carry out the regression analysis. However, if you give the command line option -GOFORIT, then the program will proceed with the analysis. I strongly recommend that you understand the reason for the problem(s), and don't just blindly use -GOFORIT all the time.
- To help disentangle the ERROR and WARNING messages (if any) from the rest of the text output, they are now also output to a file named 3dDeconvolve.err.
New option for censoring time points
- The -CENSORTR option lets you specify on the command line time points to be removed from the analysis. It is followed by a list of strings; each string is of one of the following forms:
  - 37 ⇒ remove global time index #37
  - 2:37 ⇒ remove time index #37 in run #2
  - 37..47 ⇒ remove global time indexes #37-47
  - 37-47 ⇒ same as above
  - 2:37..47 ⇒ remove time indexes #37-47 in run #2
  - '*:0-2' ⇒ remove time indexes #0-2 in all runs
- Time indexes within each run start at 0.
- Run indexes start at 1 (just be to confusing, and also to be compatible with afni_proc.py).
- Multiple -CENSORTR options may be used, or multiple -CENSORTR strings can be given at once, separated by spaces or commas.
- N.B.: Under the above rules, 2:37,47 means index #37 in run #2 and then global time index 47; it does not mean index #37 in run #2 and then index #47 in run #2. To help catch this possible misuse, the program will print a warning message if you use some -CENSORTR strings with run numbers and some without run numbers.
New program 3dSynthesize for creating 3D+time datasets
- This program combines the β weights stored in the -cbucket output from 3dDeconvolve, and the regression matrix time series stored in the -x1D output, to produce model fit time series datasets. 3dDeconvolve itself has the -fitts option to produce the full model fit in each voxel. 3dSynthesize can be used to produce model fits from subsets of the full model.
- In the examples below, suppose that fred+orig is the output from -cbucket and that fred.x1D is the output from -x1D. Also suppose that there were two stimulus classes, given labels Face and House in 3dDeconvolve using -stim_label options.
- In general, if you want to "Double Plot" the resulting dataset on top of the original time series dataset (with the Dataset #N plugin), you'll need the baseline model component so that the 3dSynthesize output is on the same magnitude scale for graphing.
- For further details, see the -help output from 3dSynthesize: available here.
- [25 Jun 2007] Censoring
  3dDeconvolve and 3dSynthesize have been modified to work when the 3dDeconvolve run using a time point censoring option (i.e., '-censor' and/or '-CENSORTR'). The matrix files output by 3dDeconvolve (which files are now renamed to end in .xmat.1D) have information about which time points were censored. 3dSynthesize can use that information to generate sub-bricks to fill in those time points which are missing in the actual matrix. The options are:
  - -cenfill zero = fill censored time points with zeros [the new default]
  - -cenfill nbhr = fill censored time points with the average of their non-censored time neighbors
  - -cenfill none = don't put sub-bricks in for censored time points [what the program used to do]
  Another option is to use the new -x1D_uncensored filename option in 3dDeconvolve to output an uncensored version of the regression matrix, then use that matrix as the input the 3dSynthesize. Then the model fit that you choose will be computed at all the time points.
Amplitude Modulated FMRI regression analysis
—OR—
Analysis of event-related FMRI data when the amplitude of each event's BOLD response might depend on some externally observed data.
Conceptual Introduction:
When carrying out an FMRI experiment, the stimuli/tasks are grouped into classes. Within each class, the FMRI-measurable brain activity is presumed to be the same for each repetition of the task. This crude approximation is necessary since FMRI datasets are themselves crude, with low temporal resolution and a low contrast-to-noise ratio (i.e., the BOLD signal change is not very big). Therefore multiple measurements of the "same" response are needed to build up decent statistics.

In many experiments, with each individual stimulus/task a separate measurement of subject behavior is taken; for example, reaction time, galvanic skin response, emotional valence, pain level perception, et cetera. It is sometimes desirable to incorporate this Amplitude Modulation (AM) information into the FMRI data analysis.

Two potential motives for using AM in the FMRI data analysis can be distinguished. The first goal is simply to consider the AM as a nuisance variable that may affect the FMRI activation amplitude in some way, and one wishes to reduce this impact on the results as much as possible. The second goal is to estimate and map quantitatively the effect of the AM on the FMRI activation amplitude.

Binary AM:
If the AM were binary in nature ("on" and "off"), one method of carrying out the analysis would be to split the tasks into two classes, and analyze these stimulus classes separately (i.e., with two -stim_times options). The statistical test for activation ignoring the AM would then be a 2 DOF F-test, which could be carried out in 3dDeconvolve by using a 2 row GLT. The contrast between the two conditions ("on−off") could be carried out with a 1 row GLT. For example:
```
  3dDeconvolve ...                                                    \
      -stim_times 1 regressor_on.1D  'BLOCK(1,1)' -stim_label 1 'On'  \
      -stim_times 2 regressor_off.1D 'BLOCK(1,1)' -stim_label 2 'Off' \
      -gltsym 'SYM: On \ Off' -glt_label 1 'On+Off'                   \
      -gltsym 'SYM: On -Off'  -glt_label 2 'On-Off' ...
```
(A realistic 3dDeconvolve command line would, of course, have more options to specify the input and output filenames, etc.) The above example assumes that each case ("on" and "off") is being analyzed with simple (fixed-shape) regression — short 1-second blocks of activity.

Nothing more will be said here about binary AM, since it is just a standard application of 3dDeconvolve; the only (small) difference is that the stimulus class to which each individual stimulus is assigned is determined during the FMRI data acquisition itself, rather than determined by the investigator before the imaging session.

Continuous AM:
More complex is the case where the AM measurement values fall onto a continuous (or finely graded discrete) scale. One form of analysis is then to construct two regressors: the first being the standard constant-amplitude-for-all-events-in-the-same-class time series, and the second having the amplitude for each event modulated by that event's AM value (or some function of the AM value). To make these two regressors be orthogonal, it is best to make the modulation be proportional to the difference between each event's AM value and the mean AM value for that stimulus class.

The new -stim_times_AM2 option is designed to make this type of analysis easy. The 'AM' in the option suffix indicates that amplitude modulation for each time is expected in the input timing file. The '2' indicates that 2 regressors will be generated from 1 stimulus timing file.

The stimulus timing file for -stim_times_AM2 has a slightly different format than the stimulus timing file for the standard -stim_times option. Each stimulus time in the _AM2 file must have an amplitude "married" to it. For example:
```
 10*5 30*3 50*2 70*7 90*-3
```
indicates that the stimuli at times 10, 30, 50, 70, and 90 have amplitudes of 5, 3, 2, 7, and -3 (respectively). Note that if a stimulus time is given without an amplitude, the amplitude will be taken to be zero and 3dDeconvolve will print a warning message. (N.B.: the '*' separator can also be the 'x' character, if that is more convenient.)

The program 1dMarry can be used to "glue" two .1D formatted files together to produce a file appropriate for -stim_times_AM2. With the -divorce option, it can also split up a "married" file into 2 separate files — one with the times and one with the amplitudes. These features makes it relatively straightforward to run a standard 3dDeconvolve analysis with -stim_times and also the new -stim_times_AM2 type of analysis.

The same response models available with the standard -stim_times option are also usable with -stim_times_AM2. Two regression matrix columns will be generated for _AM2 for each one column specified by the response model (e.g., 'BLOCK(1,1)' generates 1 column normally, and 2 columns when used with _AM2). The first column will be created by giving equal weight (1) to each event in the stimulus timing file. The second column will have each event weighted by the difference between its individual amplitude and the mean of all amplitudes in the timing file. The significance of the output β weight for this second column (e.g., given by using the -tout option) can be used to map regions that are (linearly) sensitive to the amplitude information. The significance of the combined β weights for the two columns (e.g., given by using the -fout option) can be used to map regions that are sensitive the stimulus class as a whole.

It can be useful and enlightening to plot the columns of the regression matrix that correspond to the equal-weight and variable-weight model time series generated by -stim_times_AM2. For this purpose, program 1dplot can be applied to subsets of the .x1D file output by 3dDeconvolve.

It is possible to use the option -stim_times_AM1 if you want to just generate a single regression model where each event is simply scaled by its associated amplitude. There will be no separation of the model into the constant and varying components. I do not recommend this, for reasons given below, but the option is available. (If you can think of a good reason to use this option for analysis of FMRI time series, please let me know!)

Deconvolution with Amplitude Modulation:
It is also legal to use a deconvolution model (e.g., 'TENT()') with -stim_times_AM2. However, you must realize that the program will compute a separate HRF shape for the AM component of the response from the mean component. It is not possible to specify that the AM component has the same shape as the mean component, and just has a different response amplitude — that would be a nonlinear regression problem, and 3dDeconvolve isn't that flexible. Also, at present, the -iresp option will not output the HRF for the AM component of a -stim_times_AM2 deconvolution model. Nor have I actually tried using AM deconvolution myself on real data. If you are going to try to do this, you should (a) understand what you are doing, and (b) consult with someone here.
Why Use 2 Regressors?:
One user asked the following question: "Can't I just use the AM weighted-regressor in the model by itself? Why do you have to include the standard (mean amplitude) regressor in the full model to investigate the effect of event amplitude values?" In other words, why not use -stim_times_AM1?

The reasoning behind separating the regressor columns into 2 classes (mean activation and AM-varying activation) is
- to allow for voxels where the amplitude doesn't affect the result, and
- to allow for a cleaner interpretation; in voxels where both regressors have significant weight, you can use the coefficient (β weight) of the first regressor as the mean activation level, and the coefficient of the second as the dependence of the activation level on the amplitude.
A numerical example might help elucidate:

Suppose that you have 6 events, to be simple, and that the amplitudes for these events are {1, 2, 1, 2, 1, 2}, with mean=1.5. Now suppose you have a voxel that IS active with the task, but whose activity is not dependent on the amplitude at all. Say its activation level with each task is 6, so the "activity vector" (i.e., the BOLD response amplitude for each event) is {6, 6, 6, 6, 6, 6}. This vector is highly correlated with the AM vector {1, 2, 1, 2, 1, 2} (cc=0.9486), so you will get a positive activation result at this voxel when using a single regressor. You can't tell from the regression if this voxel is sensitive to the amplitude modulation or not.

But if you use 2 regressors, they would be proportional to {1, 1, 1, 1, 1, 1} and {-0.5, +0.5, -0.5, +0.5, -0.5, +0.5} (the differences of each event amplitude from the mean of 1.5). The first regression vector is perfectly correlated with the "activity vector" {6, 6, 6, 6, 6, 6} and the second regression vector is not correlated with the activity at all. So you would get an activation result saying "this voxel was activated by the task, but doesn't care about the amplitude". You cannot make such a dissection without using 2 regressors.

Even if you don't care at all about such non-AM-dependent voxels, you must still include them if you think this may be a significant effect in the data. You have to model the data as it presents itself. In a sense, the constant-activation model is like the baseline model (e.g., -polort stuff), in that it must be included in the fit since it does occur, but you are free to ignore it as you will. Interpreting the results is your problem.

What about Losing Degrees of Freedom?:
If you are concerned about losing degrees of freedom, since you will be adding regressors but not data, then I would run the analysis twice. Once with the mean regressors only, and then one with the mean and the variable regressors. Then decide if the maps from the mean regressors in the two cases differ markedly. My guess is that they will not, if you have a decent number of events in each case (30+). If they do not differ too much, then you are safe to use the double regressor (AM2) analysis. If they do differ a lot (e.g., you lose a lot of mean regressor activation when you set the F-statistic p-values the same), then you probably can't use the double regressor analysis. But it is easy enough to try.

You can open two AFNI controllers, and view the single and double regressor analyses side-by-side. You can set the threshold sliders to be locked together in p-value (using Edit Environment on variable AFNI_THRESH_LOCK). This should help you decide very quickly if the two results look the same or not — same, that is, from the viewpoint of interpreting the results. The maps will of course not be identical, since they will have been calculated with different models.

Caveats and Issues:
One problem with the above idea is that one may not wish to assume that the FMRI signal is any particular function of the event amplitude values. I don't know at this time how to deal with this issue in the context of linear regression. For example, the "linear in event amplitude" model could be extended to allow for a quadratic term (-stim_times_AM3?), but it is highly unclear that this would be useful. Some sort of combination of regression analysis with a mutual information measurement might be needed (to quantify if the BOLD response is "predictable" from the AM information), but I don't fully know how to formulate this idea mathematically.

15 March 2007 — rwcox@nih.gov