Regression with Event-dependent Amplitudes
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. 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 auxiliary behavioral information (ABI) into the FMRI data analysis.
Two potential motives for using ABI in the FMRI data analysis can be distinguished. The first is simply to consider the ABI as 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 is to estimate quantitatively the effect of the ABI on the FMRI activation amplitude.
Binary ABI:
If the ABI were binary in nature ("on" and "off"), one method of carrying out
the analysis would be to split the tasks into two groups, and analyze these
separately. The statistical test for activation ignoring the ABI 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_file 1 regressor_on.1D -stim_label 1 'On' \ -stim_file 2 regressor_off.1D -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; presumably the files regressor_on.1D and regressor_off.1D were produced by program waver. If deconvolution were being performed, the GLT commands would be somewhat more complicated, depending on which lags the user wanted to keep in the test.
Continuous ABI:
More common is the case where the ABI 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 time series, and the second having
the amplitude for each event modulated by the event's ABI value
(or some function of the ABI value). To make these two regressors
be orthogonal, it is best to make the modulation be proportional to
the difference between each event's ABI value and the mean ABI value.
A standard way to generate the first regressor with waver is a command like
waver -GAM -peak 1.0 -TR 2.0 -tstim 8.3 17.5 21.9 ... > reg1.1DA new feature of waver (added 13 May 2005) lets the user scale the amplitude of each individual response; for example:
waver -GAM -peak 1.0 -TR 2.0 -tstim 8.3x1.7 17.5x-2.3 21.9x0.6 ... > reg2.1DHere, the scaling factor follow each time with the "x" character as a separator. (I thought of using the "*" character, but that would mean each time*amplitude combination would have to be enclosed in quotes, to prevent the shell from treating the "*" as a filename wildcard.) For the second regressor, the amplitude for the i-th event should be proportional to the i-th ABI value minus the mean ABI value.
N.B.: The 'xAMPLITUDE' feature was added to waver's -tstim option on 13 May 2005, and will not work in earlier editions of AFNI.
With this scheme, then a GLT can be used to test for the combined effect of both regressors (just as in the 'On&Off' case above). The default statistics for reg1.1D and for reg2.1D will provide the tests for the mean activation and the ABI-dependent activation, respectively.
Deconvolution with ABI-weighting:
If deconvolution is being computing using
-stim_maxlag>0, normally the -stim_file input
is a column of 0s and 1s. However, it is perfectly
legitimate to replace the 1s with any set of nonzero value. In this
application, you would again use two -stim_file options.
The first file would comprise 0s and 1s as normal. The second file
would replace the 1s with the de-meaned ABI amplitudes. You then have
to construct meaningful GLTs to test the results.
N.B.: I have not tried this, so
be sure to output the -fitts and -iresp files
to make sure you are getting reasonable results!
N.B.: At present, there is no
way to modulate the -stim_times inputs. Maybe someday.
Why Use 2 Regressors?:
One user asked the following question:
"Can't I just have the ABI weighted-regressor in model, or do you have to include the standard regressor in the full model to investigate the effect of ABI measurement values?" My answer follows.
The reasoning behind separating the regressors into 2 classes (mean activation and ABI-varying activation) is
- to allow for voxels where the ABI 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 ABI. However, if the second regressor has all positive ABI weights, then its weight will include some component of the mean activation. The net activation reported by the GLT summing the 2 regressors should be the same, since the two regressors would span the same linear space regardless of whether the mean is removed from the ABI scaling factor. But the interpretation of the individual activations or the difference between them is clouded when the second regressor isn't constructed as I described.
Suppose that you have 6 events, to be simple, and that the ABI values for these events are {1,2,1,2,1,2}, with mean ABI=1.5. Now suppose you have a voxel that IS active with the task, but whose activity is not dependent on the ABI at all. Say its activation level with each task is 6, so the activation vector is {6,6,6,6,6,6}. This vector is highly correlated with {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.
But if you use 2 regressors, they would be like {1.5,1.5,1.5,1.5,1.5,1.5} and {-0.5,+0.5,-0.5,+0.5,-0.5,+0.5}. The first is perfectly correlated with the "data vector" {6,6,6,6,6,6} and the second is not correlated with the data at all. So you would get an activation result saying "this voxel was activated by the task, but doesn't care about the ABI". You cannot easily make such a dissection without using 2 regressors.
Even if you don't care at all about such voxels, you must still include them. 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.
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 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.
Caveats and Issues:
One problem with the above idea is that one may not wish to assume that the
FMRI amplitude is any particular function of the ABI value. I don't
know at this time how to deal with this issue.