3dMSS with contrasts at individual level

[cltrtech]

AFNI version info: Precompiled binary linux_ubuntu_16_64: Feb 24 2026 (Version AFNI_26.0.09 'Pupienus Maximus')

Hello,

I'm testing the 3dMSS group-level modelling and had a few questions.

My first question concerns implementing and interpreting the predictions from contrasts at the individual level between groups. For example, in this post, there is a between-subjects factor pain modeled in 3dMSS and at the individual-level in 3dDeconvolve the stimulus is contrasted and modeled individually. Suppose for example's sake that there is a second within-subjects factor session at the individual level, modeling sessions of the task (the reasoning perhaps being to model habituation or whatever is reasonable for the example). So, there are multiple contrasts at the individual-level. My interpretation from the mentioned post (and this one) is that in the 3dMSS model one would run separate models to estimate each HRF for each session (session1 & session2) and to estimate HRF for each stimulus (taste, tasteless) using the individual-level data in the datatable. Additionally, to estimate contrasts at the group-level, one would run models on the individual-level contrasts as well (taste-tasteless_session2-session1 - how HRF for taste-tasteless changes from session1 to session2).

Am I correct in thinking of implementing such a complex model in this way? If not, where does the error arise? If so, a follow-up is would one have to use multiple prediction tables for the more complex contrasts or would using the HRF predictions from the simpler models suffice? Hopefully that makes sense.

My second question concerns \beta_0 and \beta_1 in the 3dMSS, specifically how they're detailed in Chen et al. (2023). The following quote caught my attention

One trivial case is the following: if only the first two splines of \beta_0(t) = 1 and \beta_1(t) = t are adopted, the smooth fitting model (3) degenerates into a linear regression: As ℎ′′( 𝑡 ) = 0 in this case, the penalty (4) becomes zero and no regularization is applied.

I haven't yet read any of the generalized additive model (GAM) literature, so the technical details are lost on me, but when the smooth fitting model degenerates how should the resulting HRF be treated? I noticed that in those cases the HRF is a linear function, although there are some between cases where it resembles a polynomial function as well. Would it not be classified as an HRF altogether? That is, would it be a simple linear or polynomial regression of the BOLD signal at each TR - degenerating into a fixed effects model with each TR as a separate predictor? Admittedly, I'm not committing to any dichotomy I'm implying here (prior two questions).

-Chris

[Gang]

Chris,

Regarding your first question, suppose you have four factors: (1) groups (G1 and G2), (2) task conditions (A1 and A2), (3) sessions (S1 and S2), and (4) hemodynamic response profiles.

A conventional modeling approach would treat this as a standard ANOVA-style factorial structure and extract post hoc estimates from a single full model. However, given the complexity of modeling smooth splines in 3dMSS, each research question must be addressed individually. For instance, to assess the interaction between session and task condition, you would first calculate the cross-session difference of the condition contrast for each participant (i.e., (S1A1 - S1A2) - (S2A1 - S2A2)) and then feed those contrast values into 3dMSS at the group level. Consequently, each specific research question will require its own model and corresponding prediction table.

Regarding your second question about using only the first two basis functions: this minimal setup is generally not practical for modeling hemodynamic responses. Because we know a typical hemodynamic response involves an initial slope, a peak, and a post-stimulus undershoot, a simple linear trajectory cannot capture this non-linear profile. As a general guideline, at least 5 smooth spline basis functions are needed to adequately reconstruct the response curve.

Gang Chen

[cltrtech]

I see, so for a complex factorial model one could target specific questions about the different levels and interactions using a variety of prediction tables.

So, my follow up question to that concerns designing the prediction tables. Suppose I've modeled the (S1A1-S1A2)-(S2A1-S2A2) contrast at the individual level and would like to examine how groups G1 and G2 differ in their HRFs from S1A1-S1A2. To examine how G1's contrast of (S1A1-S1A2) differs from G2's contrast, I imagine I would design 2 different prediction tables - one for each task condition during session S1. Then so on for every given question. Though, I'm lost at this step. In the prediction table below, am I to think that the pred column would only need to specify the specific levels identified in the individual-level data? Presumably it's case-sensitive to name the individual-level filenames.

pred	Groups	TR
G1.S1A1.TR0	0	0
G1.S1A1.TR0.25	0	0.25
G1.S1A1.TR0.5	0	0.5
....
G2.S1A1.TR0	0	0
G2.S1A1.TR0.25	0	0.25

If this is the case, the naming at the individual-level would have to be as meaningful as possible to identify the variable/values of interest. I don't imagine including predictors for contrasts levels in the datatable or 3dMSS function makes any sense (each datatable observation/row is a contrast so no specific level could be specified). Perhaps I'm overlooking a detail if this isn't correct.

-Chris

[Gang]

Hi Chris,

When using 3dMSS for group-level hemodynamic response analysis, the program serves two distinct primary purposes:

1. Statistical inference

This output provides statistical evaluations for group-level effects (e.g., presence of an overall group response or group comparisons across the response profile), expressed via \chi^2-statistics.

For example, with a three-factor design (group \times session \times condition), you can evaluate the three-way interaction by calculating the double-difference contrast per participant at each TR (i.e., (S1A1 - S1A2) - (S2A1 - S2A2) ) and feeding those discrete time-series responses into 3dMSS with a model specification like:

-mrr 's(TR,k=10) + s(TR,k=10,by=group) + s(TR,Subj,bs="fs",m=1)' \
-vt Subj 's(TR,Subj)' \
-qVars 'TR' \

(This assumes individual-level responses are estimated on the TR grid).

2. Estimating & visualizing group-level response profiles

Here, the goal is to reconstruct and visualize the full continuous response shape (e.g., in the AFNI GUI Graph window). This is specified using the -prediction option alongside a prediction table.

To obtain group-level curves for each experimental cell, you would run 4 separate models (one for each combination S1A1, S1A2, S2A1, and S2A2). For instance, using the individual-level S1A1 sub-bricks as input, you would append:

-prediction pred_S1A1.txt \

where pred_S1A1.txt has a format such as:

Label            Group   TR
G1.S1A1.TR0         -1    0
G1.S1A1.TR0.25      -1    0.25
G1.S1A1.TR0.5       -1    0.5
...             
G2.S1A1.TR0          1    0
G2.S1A1.TR0.25       1    0.25
G2.S1A1.TR0.5        1    0.5
...

Note: The format of the prediction table remains identical across all four models (S1A1, S1A2, S2A1, S2A2), except for updating the labels in the first column.

Let me know if this clarifies the workflow!

Gang Chen

[cltrtech]

Hello Gang,

Yes! It does clarify the workflow more and also thank you for including the addition of s(TR,Subj,bs="fs",m=1) in the model specification (for anyone reading and as confused as me when noticing the "fs" parameter value here's the documentation for its meaning in mgcv). I think my trouble was with distinguishing the two purposes. How you've described the purposes is how I originally thought about using the model but predictions from the double-difference contrast model didn't make sense unless conceptualized as HRF difference-in-differences (which are useful but not "proper" HRFs, I think). I have some more follow-up questions.

Can the full capabilities of the mgcv smoothing terms be utilized in 3dMSS? Seems like a silly question since you used it and the 3dMSS.R file seems to suggest as much, but it's worth asking, I think. If so, in the you mention two points regarding TRs, experimental design and HDR estimation in HDR detection sensitivity post. The second one catches my attention regarding flexible use of smoothing terms. In those cases that shorter TRs are challenging to obtain, could the use of different smooth terms also facilitate achieving "proper accuracy in estimating HDR profiles"?

My second follow up concerns the output of the 3dMSS specification you provided above. I notice that the results provide subject-level evidence similar to specifying a (1|subject) term in mixed models. So, could individual-level HRF profiles for subjects as they differ within groups be examined as well using this specification? Provided I design the proper prediction table?

-Chris

[Gang]

Chris,

Here are my thoughts on your follow-up questions:

Can the full capabilities of the mgcv smoothing terms be utilized in 3dMSS?

TBH, I haven't fully explored every smoothing specification available in the R package mgcv. So far, we have primarily focused on the most common smoothers applicable to hemodynamic response modeling. If you have specific smooth terms or basis combinations in mind that could be beneficial, please feel free to share them!

In those cases that shorter TRs are challenging to obtain, could the use of different smooth terms also facilitate achieving "proper accuracy in estimating HDR profiles"?

It might be most efficient if you could share a concrete example, such as your specific TR length and the estimated response duration you are working with. That way, we can discuss the most appropriate smoothing strategy for your temporal resolution.

I notice that the results provide subject-level evidence similar to specifying a (1|subject) term in mixed models. So, could individual-level HRF profiles for subjects as they differ within groups be examined as well using this specification? Provided I design the proper prediction table?

I assume that you're referring to the statistic output in the current 3dMSS implementation associated with the term s(TR,Subj,bs="fs",m=1) in the model specified in my previous message. The current statistical output evaluates the overall extent of variability in response profiles across individuals. This is equivalent to the statistical evidence regarding the standard deviation parameter \sigma associated with the distributional assumption \mathcal{N}(0, \sigma^2) in a conventional linear mixed-effects model, evaluating the random-intercept variance components in the (1|Subj) term, rather than estimating each individual's unique trajectory.

If your goal is to extract individual-specific response profiles, that is, obtaining each participant's specific deviation from the group mean (or each individual's predicted response) via partial pooling/shrinkage, it is conceptually possible, but it is not yet implemented in the current version of 3dMSS.

Gang Chen

[cltrtech]

Hello Gang,

I've not studied mgcv package/GAMMs yet, but once I do, and if I encounter any insights, I'll share them.

For the concrete example, the specific TR length is 2.38s and the estimated response duration I've been using is modeled by the TENT function 'TENT(0,16,9)', although that function was intended for a TR of 2s. So, I imagine it's slightly inappropriate for the 2.38s.

-Chris

[Gang]

Hi Chris,

Are the stimulus onset times randomized across participants? Also, are you using a relatively short trial duration (e.g., an event-related design)?

Assuming there is no instantaneous response at stimulus onset, you could specify the response model with a TR of 2.38 s as:

TENT(2.38, 16.66, 7)

To improve the sampling resolution at the group level, you might try a simple analytical trick: estimate the HDR on the standard TR grid for half of the participants, but artificially shift the real stimulus onset times by +0.5 TR for the other half. This effectively samples the response curve at half-TR offsets for the second subgroup. When combined in 3dMSS, the model fits across a grid with double the temporal resolution, sampling at TR boundaries for one half and mid-TR points for the other.

Gang Chen

[cltrtech]

Hello Gang,

The trials are not randomized across participants, but they do differ between versions of the task (two versions) and the ITIs are jittered. It is also an event-related design of 4.5s duration with the stimulus offsetting at 5s. I noticed the analytic trick in other posts. I'm going to try the TENT function above with the trick, unless the experiment design details justify a different TENT function setup.

I have one last question regarding HRF interpretation. Given what is known about a typical hemodynamic response, when estimated HRFs degenerate into linear functions wouldn't this be compatible with concluding that even if say G1 and G2 differed in their BOLD signals (say different slopes), there is not sufficient evidence to conclude that their hemodynamic responses are likely different? This is my reading of such cases so far, but I'm wondering if it's true to the underlying framework/theory.

-Chris

[Gang]

Chris,

Given what is known about a typical hemodynamic response, when estimated HRFs degenerate into linear functions wouldn't this be compatible with concluding that even if say G1 and G2 differed in their BOLD signals (say different slopes), there is not sufficient evidence to conclude that their hemodynamic responses are likely different?

Here is how I view this scenario--

You’ve raised a crucial question about whether we should let the statistical tail wag the scientific dog. It is theoretically possible for two groups to show a statistically "significant" difference in the form of linear trends across time. However, because we know from underlying physiology that a true BOLD response features a peak followed by an undershoot, a purely linear response profile indicates a fundamental model mis-specification or a non-biological artifact.

Even if the statistical evidence for a group difference in slopes is strong, such a finding lacks biological plausibility as a hemodynamic response. Rather than interpreting it as a meaningful difference in BOLD signaling, it would warrant caution, at most, something to shelve for further diagnostic investigation rather than present as a true hemodynamic effect.

Gang Chen

[cltrtech]

Hello Gang,

I think I'll use that version of "wag the dog" in the future if you don't mind and explore this diagnostic line of inquiry a bit later. Thank you for all the help!

-Chris