Thursday, September 27, 2012

Duration Regressors with fMRI

This subject was brought to my attention by a colleague who wanted to know whether parametric modulation or duration modulation was a better way to account for RT effects. While it can depend on the question you are trying to answer, often duration modulation (referred to here as a "variable epoch model") works best. The following highlights the different approaches for modeling trials which involve a period of decision-making or an interval between presentation of a stimulus and the resulting response.


Over the past few years there has been a renewed interest in modeling duration in fMRI data. In particular, a methods paper by Grinband and colleagues (2008) compared the effects of modeling the duration of a trial - as measured by its reaction time (RT) - against models which used RT as a parametric modulator and against models which did not use RT at all. The argument against using RT-modulated regressors was that, at short time intervals (i.e., less than four seconds), using an impulse function was a good approximation to the resulting BOLD signal (cf. Henson, 2003).

Figure depicting relationship between stimulus onset, BOLD response, and underlying neural activity. Red lines: activity associated with punctate responses (e.g., light flashed at the subject). Blue lines: Activity associated with trials of unequal duration (e.g., decision-making).

However, for a few investigators, such assumptions were not good enough. To see whether different models of RT led to noticeable differences in BOLD signal, Grinband et al (2008) examined four types of modeling:
  1. Convolving the onset of a condition or response with the canonical HRF (constant impulse model);
  2. Separately modeling both the main effect of the condition as well as a mean-centered parametric modulator - in this case RT (variable impulse model);
  3. Binning each condition onset into a constant amount of time (e.g., 2 seconds) and convolving with the canonical HRF (constant epoch model); and
  4. Modeling each event as a boxcar function equal to the length of the subject's RT (variable epoch model).

Graphical summary of models from Grinband et al (2008). Top: Duration of cognitive process as indexed by reaction time. Constant Impulse:  Onset of each event treated as a punctate response. Constant Epoch: Onset of each event convolved with boxcar function of constant duration. Variable Impulse: Punctate response functions modulated by mean-centered parameter (here, RT). Variable Epoch: Each event modeled by boxcar of duration equal to that event's RT.

Each of these models was then compared using data from a decision-making task in which subjects determined whether a line was long or short. If this sounds uninteresting to you, you have obviously never done a psychology experiment before.

The authors found that the variable epoch model - in other words, convolving each event with a boxcar equal to the length of the subject's RT for that trial - captured more of the variability in the BOLD response, in addition to reducing false positives as compared to the other models. The variable epoch model also dramatically increased sexual drive and led to an unslakeable thirst for mindless violence. Therefore, these simulations suggest that for tasks requiring time - such as decision-making tasks - convolution with boxcar regressors is a more faithful representation of the underlying neuronal dynamics (cf. the drift-diffusion model of Ratcliff & McKoon, 2008). The following figures highlight the differences between the impulse and epoch models:

Comparison of impulse models and epoch models as depicted in Grinband et al (2008). A) For impulse models, the shape remains constant while the amplitude varies; for epoch models, increasing the duration of a trial leads to changes in both shape and amplitude. B) Under the impulse model, increasing the duration of a stimulus or cognitive process (as measured by RT) leads to a reduction in explained variance.

Figure from Grinband et al (2008) showing differential effects of stimulus intensity and stimulus duration. Left: Increasing stimulus intensity has no effect on the time to peak of the BOLD response. Right: Increasing stimulus duration (or the duration of the cognitive process) leads to a linear increase in the time for the BOLD response to peak.


One caveat: note well that both parametric modulation and convolution with boxcar functions will account for RT-related effects in your data; and although the Grinband simulations establish the supremacy of boxcar functions, there may be occasions that warrant parametric modulation. For example, one may be interested in the differences of RT modulation for certain trial types as compared to others; and the regressors generated by parametric modulation will allow the researcher to test them against each other directly.

27 comments:

  1. Great post. We did this recently to account for RT affects and it appeared to work well. It's worth noting that one of the issues with the parametric modulation approach is it's a pain to implement as you can't have separate modulators for each condition (since they then only account for RT effects within condition but not across). The RT as duration approach is simple to implement and allows one to keep the same design.

    ReplyDelete
    Replies
    1. Thanks for the information about parametric modulation across conditions; I was unaware of that fact!

      Delete
  2. I've just started to use the variable epoch model on a dataset where the RT varies both systematically across conditions, and to a much smaller degree within condition. The most interesting results so far is that my sexual appeal has increased considerably, and I can now tell the difference between various species of finch. Will keep you updated.

    ReplyDelete
    Replies
    1. Nice! I'm always glad to hear when that happens.

      Delete
  3. Great post. I found it very useful. Thank you!

    ReplyDelete
  4. Was there a reason (perhaps maybe due to a shortcoming in the analysis software used) why Grinband used a variable impulse function with a CONSTANT RT (i.e. averaged across trials), as opposed to using each trial's RT? This would be like the variable epoch they used but with a stick function (duration 0) instead of a boxcar.
    I read the paper several times in length but I could not find anything related to this...

    Finally, what about other possible models, i.e. using the trial-by-trial stimulus length as the duration of the cognitive process instead of RT (and from duration zero, as the stick function assumes) AND using RT as a parametric modulator?

    Thanks!

    ReplyDelete
  5. I am not sure I agree with this for all cases. I work on high level cognitive functions, and I think there are some processed which are modulated by RT and some which are not. This becomes very difficult if you have reaction time differences between conditions: you are explicitly modeling an amplitude difference in which the condition with longer RT has a larger modeled response. Thus when you run a simple A > B contrast you may get significance if the RT for B was larger than A, because B will then have a smaller beta value in a GLM model to equate to equal magnitude of BOLD signal changes in each condition.

    Like so many things in neuroimaging, the correct choice is dependent upon the circumstances!

    ReplyDelete
    Replies
    1. Hi Colin,

      I agree that each case is different, and that you should try modeling it different ways. However, if RT is an index for the duration of the cognitive process, and this is reflected by the BOLD response, then the variable duration model should best explain the observed timecourse. If you try to model a relatively small and brief BOLD response that isn't modulated by RT, then the beta estimate for that condition should be lower to reflect the difference.

      This would be more of a problem with parametric modulation by RT, where explained variance would need to be partitioned between the mean response of each condition and any modulation beyond that.

      Best,

      -Andy

      Delete
  6. Hi Andy,
    We are just starting to analyze data from a Mind Wandering Study and am seeking some advice. We used a target detection task with infrequent targets to engender mind wandering. The subjects were required to press a button on the appearance of the target. Our behavioural data shows that there is a relationship between the RT and how much the subject felt they had drifted off task. We would like to examine how the Bold activation changes pre- to post- target appearance, taking into account RT (as a measure of how 'far' the participant was off task). Have you any advice about how to model the RT in this case? I use FSL.

    David Watson

    ReplyDelete
    Replies
    1. Hi David,

      You could model the RT as a boxcar function, convolving a boxcar with how long it took the subject to make a response, time-locked to when the target appears.


      Best,

      -Andy

      Delete
  7. Thanks Andy
    If I wanted to do this in FSL with a three column custom EV would it be something like:
    onset duration (the RT) weighting

    with the weighting set to one for each event?

    David

    ReplyDelete
    Replies
    1. Hey David, that's correct! That should capture any neural activity correlated with the reaction time for that event.

      -Andy

      Delete
    2. Hi Andy, your suggestion works very well on my data set. I have another question which is also associated with using weighting.
      We have probed individuals just after target detection on three VAS scales related to attention. Each gives us an 'index'. Although these indexes are related to some extent they tend to focus on different aspects of their pre-target 'state'. Is there some way I can code these into an event structure? For example if I have an estimate of how much they were worrying about past or future events immediately prior to target detection. Would this just be entered as a value (like) RT with a weighting of 1? Could I do this for each of the indices separately? Does that make any sense?
      Thanks in advance. David

      Delete
    3. Hi David,

      It sounds as though you would want to model those indexes as parametric modulators. If you believe that they represent some psychological measurement that would affect their processing of the target, you could put all three in as modulators and see how they capture parts of the signal above and beyond that caused by the target.

      If you include more than one parametric modulator on a regressor, however, make sure that you mean-center it and that you don't orthogonalize it; there have been cases where orthogonalizing gives the most shared variance to the first modulator, and can distort your interpretation of the results: http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0126255

      -Andy

      Delete
    4. Hi Andy, thanks for your reply, that's exactly right. I wasn't completely sure how to set up the events for this and develop the appropriate contrasts. So:
      I create 1 unmodulated event file with weighting 1 for each event (simple unmodulated model).
      To look at possible modulation effects I create a modulated event file for each of my parametric modulators with weighting set by event (mean centred).
      I can then evaluate contrasts for all of these both 1 and -1.
      Hope this is correct!
      Finally, does it make any sense to look for differences between the unmodulated event and the parameter modulated events?

      Delete
    5. Hi David,

      You can include both the unmodulated event and the parametric modulators into the same model. Just specify a new parametric modulation for each condition you want to modulate, enter the values that modulate each event of that condition, and then estimate the model; you will get parameter estimates for both the regressor and the modulator. If you compare the beta maps of your main regressors between a model that has parametric modulators and one that doesn't, they are actually very similar; the only differences come from the parametric modulators soaking up additional variance, and thus affecting parameter estimates for the main regressors. You can set up contrasts the same way you would without the modulation, but you will have more regressors. Make sure that the weights line up with the regressors you're trying to contrast.

      You can take contrasts between the modulated event and the modulator, but I'm not sure how to interpret that. Ditto for taking contrasts between parametric modulators. Usually researchers are only interested in the effects of individual modulators.

      -Andy

      Delete
  8. This comment has been removed by a blog administrator.

    ReplyDelete
  9. Hello Andrew, thank you for the post! I have been trying to understand how does SPM12 create the PSTH plots for event-related responses (HRF). Could you please direct me to a tutorial that explains how does SPM extracts the time series to create it? Thank you! -Zhamilya

    ReplyDelete
    Replies
    1. Hi Zhamilya,

      I don't know the specifics about how the PSTH is generated, and I haven't been able to find any documentation on how it works. However, my best guess is that it does something like the following:

      1. Averages activity estimates at each TR for the deconvolved HRF for each occurrence of that condition;

      2. Calculates the SD for the activity at each TR;

      3. Plots both the mean activity and error bars for each TR, along with its best estimate of the shape of the HRF for that condition.

      -Andy

      Delete
  10. Hi Andrew, What is the way that I can simulate this using Matlab?

    ReplyDelete
  11. Hi Andrew, What is the way that I can simulate this using Matlab?

    ReplyDelete
    Replies
    1. Hi Dhinesh,

      I don't know how to simulate it using Matlab; I'll think about this.

      -Andy

      Delete
  12. Hi Andy,

    I'm new to neuroimaging and was hoping that you can help. In our study we have two groups with significantly different RTs and I am wondering if this is playing a role in the increased BOLD signal in the group with the longer RTs. I see that there are several options, but I'm not sure how to set up the GLM for any of them. I've tried entering the RTs as parametric modulators, but am still unsure how to set up the Contrast Manager to test my hypothesis. Your help would be appreciated

    -Greg

    ReplyDelete
    Replies
    1. Hi Greg,

      Once you enter the RTs as parametric modulators, you can weight the parametric modulator by one divided by the number of sessions (e.g., five sessions would equal a weight 1/5, or 0.2, for each parametric modulator regressor). This will produce a contrast image for each subject, which can then be submitted to a second level analysis.

      This approach calculates beta weights representing how much the variation of that modulator - in this case, RT - explains variability in the BOLD response beyond the unmodulated regressor. Voxels that systematically vary with the duration of RT will show a cleaner distinction between variability attributed to the unmodulated regressor and the parametric modulator when parametric modulators are included.

      Note that simply modeling RTs as a parametric modulator (or as a duration regressor) won't be able to determine whether the significance of the modulation is due primarily to shorter RTs being associated with lower amplitude of the BOLD response, or because longer RTs are associated with higher amplitude of the BOLD response, or both. If you want to make that claim, you could do one of the following:

      1. Separate your groups based on their having significantly different RTs, and compare them directly using a two-sample t-test. Look at the results when both groups having their RT modeled, and when they do not, in case modeling RT has a significant effect on one group as opposed to the other. If they are matched on all other characteristics - i.e., if RT is not confounded with anything else - any differences in activity can be attributed to differences in RT.

      2. Treat the entire sample as one group, and split your trials based on median RT into high-RT and low-RT trials. If you compare the high-RT against the low-RT trials and see a significant difference, one interpretation is that RT is driving the effect.

      See Grinband's 2011 Neuroimage paper for more details on their methods for accounting for RT differences.

      Hope this helps!

      -Andy

      Delete
    2. Hi Andy,

      Thank you for your help! I did as you suggested in the first paragraph. Although as I found out, with SPM12 you no longer have to divide by the number of sessions. Just putting in ones work.

      I got some interesting results, I think. When comparing two groups, the opposite activation pattern emerged when looking at the RT modulated and un modulated regressor. There were two levels of the independent variable (that were of interest) and we looked at the response to one compared to another (Stroop: Incongruent > Congruent). This is how I set up the contrasts that fed into the second level analyses.

      [-1 0 1 0]
      [0 -1 0 1]
      (the actual contrast contained zeros for the levels that were not of interest)

      If Group1 shows increased activation in the first contrast, but Group2 shows increased activation in the second (same brain regions), would a correct interpretation be that the results from the first analysis were primarily due to RT, but when RT is accounted for Group2 actually has increased activation?

      As a side note, I was surprised that the results from the first analysis were slightly different after redoing it with RT parametric modulators. Without the extra column the contrast just looked like this: [-1 1]. Does this indicate that I made a mistake somewhere? From what I read elsewhere parametric modulators do not do anything to analyses using the original columns; is this correct?

      Thank you also for pointing me to Grinband's article. I looked at it and Yeung's reply. If my interpretation of my results is correct, then there is still an effect even when RT is accounted for.

      Thanks again!

      Delete
  13. Hi Andrew,

    I want to regress out RT variance. I used the Posner cueing task. I made my GLM using 3 column format as follows.
    EV1 = valid target: onset time/ target duration / weight=1
    EV2 = invalid target: onset time/ target duration / weight=1
    EV3 = valid target RT: onset time(valid target)/ standardized RT / weight=1
    EV4 = invalid target RT: onset time(invalid target)/ standardized RT / weight=1

    To regress out RT variance, I made contrast like [-1, 1, 0, 0].
    I wonder if I did correctly.

    Thanks in advance.

    Sun-woo

    ReplyDelete
  14. Hi Andrew,

    I want to regress out RT variance. I used the Posner cueing task. I made my GLM using 3 column format in FSL as follows.
    EV1 = valid target: onset time/ target duration / weight=1
    EV2 = invalid target: onset time/ target duration / weight=1
    EV3 = valid target RT: onset time(valid target)/ standardized RT / weight=1
    EV4 = invalid target RT: onset time(invalid target)/ standardized RT / weight=1

    To regress out RT variance, I made contrast like [-1, 1, 0, 0].
    I wonder if I did correctly.

    Thanks in advance.

    Sun-woo

    ReplyDelete