Using 3dLMEr for fMRI prediction of symptom data

AFNI version info (afni -ver): AFNI_23.0.07


We are using pre-treatment brain activations to predict symptom change over time during treatment. Our 3dLMEr script is pasted below, with two subjects' (deidentifed) data shown as an example. The variable InputFile points to "blocks" of symptom data that we created for this purpose; for example, if subject 101 at timepoint 1 had a GAD-7 score of 5, the file GAD7/101-GAD7-Wk1+tlrc would have a value of 5 at every voxel. The variable VoxelActivity_T1 points to our pre-treatment fMRI data at the sub-brik of interest.

The 3dLMEr script ran successfully, and we are now wondering what ACF parameters to use in 3dClustSim. We ran 3dFWHMx on the residuals from 3dLMEr, then put this into 3dClustSim, and got a very low minimum cluster size. We then ran 3dFWHMx on each individual stats file and put the averaged parameters into 3dClustSim instead, which gave us a somewhat more reasonable minimum cluster size, but we are still looking at a lot more significant clusters than we expected. We have seen this across multiple different fMRI tasks, and it is making us wonder if we should be thinking differently about significance thresholding in this type of analysis. Since we are predicting "perfectly smooth" data, e.g., the block of 5's, does this mean we should choose ACF parameters in a different way? Or have we gone wrong somewhere else?

Thank you,

3dLMEr -prefix _MID_WB-GAD7_a_LME_Num \
-jobs 48 \
-model 'Time*Treatment*VoxelActivity_T1+enorm+GAD7_T1+(1|Subj)' \
-resid errts._MID_WB-GAD7_a_LME_Num \
-qVars 'enorm,Time,GAD7_T1' \
-qVarsCenters '0.1,5.99,11.97' \
-vVars 'VoxelActivity_T1' \
 -SS_type 3 \
 -mask MNI_mask_resampled+tlrc \
 -dataTable  \
 Subj enorm Time Treatment VoxelActivity_T1 GAD7_T1 InputFile \
 101   0.1298   1   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk1+tlrc \
 101   0.1298   2   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk2+tlrc \
 101   0.1298   3   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk3+tlrc \
 101   0.1298   4   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk4+tlrc \
 101   0.1298   5   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk5+tlrc \
 101   0.1298   6   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk6+tlrc \
 101   0.1298   7   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk7+tlrc \
 101   0.1298   8   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk8+tlrc \
 101   0.1298   9   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk9+tlrc \
 101   0.1298   10   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-Wk10+tlrc \
 101   0.1298   11   BA   ../../stats.101-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   10   GAD7/101-GAD7-PostWDLOCF+tlrc \
 102   0.0643   1   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk1+tlrc \
 102   0.0643   2   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk2+tlrc \
 102   0.0643   3   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk3+tlrc \
 102   0.0643   4   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk4+tlrc \
 102   0.0643   5   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk5+tlrc \
 102   0.0643   6   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk6+tlrc \
 102   0.0643   7   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk7+tlrc \
 102   0.0643   8   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk8+tlrc \
 102   0.0643   9   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk9+tlrc \
 102   0.0643   10   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-Wk10+tlrc \
 102   0.0643   11   BA   ../../stats.102-T1-_MID+tlrc'[a.rall_GLT#0_Coef]'   8   GAD7/102-GAD7-PostWDLOCF+tlrc \


Before discussing cluster-level inference, I'd like to address some modeling aspects in your 3dLMEr specifications.

  1. GAD-7 Score Relationship: Since the GAD-7 score is the response variable (last column of your data table), could you clarify its relationship with GAD7_T1?

  2. Within-Individual Variables: Are both Time and VoxelActivity_T1 considered within-individual variables? If they are, it might be worthwhile to modify the random effects structure from (1|Subj) to (Time + VoxelActivity_T1|Subj). This adjustment accounts for the slope-varying effects and potential correlations between these variables within each individual.

  3. Treatment as a Between-Individual Factor: Assuming that Treatment is a between-individual factor, do the variables enorm and GAD7_T1 vary across different treatment groups? If so, I recommend exploring centering techniques, as discussed in this resource.

  4. Complexities of Covariate Inclusion: The inclusion of covariates is often more intricate than commonly recognized. For details, refer to this preprint.

Gang Chen

Hi Gang,
Thank you for your reply! Responses below.

  1. Regarding GAD-7 scores: "GAD7_T1" represents pre-treatment GAD7 scores. Because we covary for pre-treatment GAD7 scores, the pre-treatment timepoint is not included in our "Time" variable. Time=1 corresponds to Week 1 of treatment, Time=2 corresponds to Week 2 of treatment, etc.

  2. Regarding VoxelActivity_T1: For a given subject, this is the same sub-brik in each row, representing pre-treatment voxel activity. To back up a bit, our effects of interest are the VoxelActivity_T1Time interaction, which indicates the voxels in which pre-treatment activity predicts symptom change, and the VoxelActivity_T1Time*Treatment interaction, which indicates the voxels in which pre-treatment activity predicts symptom change differentially across treatment-arms. I don't think VoxelActivity_T1 should be a considered a within-individual variable, unless I'm misunderstanding what 3dLME is doing?

  3. Just checked again, enorm (head motion) and GAD7_T1 (pre-treatment symptoms) do not vary significantly across treatment groups.

  4. Thank you for sending the preprint! Our rationale for covarying for pre-treatment symptom score (GAD7_T1) is that we wanted to identify how pre-treatment brain activity could predict trajectory of symptoms regardless of pre-treatment symptoms. We conceptualize the pre-treatment symptoms as a confounder: pre-treatment brain activity causes pre-treatment symptoms, and both pre-treatment brain activity and pre-treatment symptoms cause subsequent symptoms.



To clarify, it seems you intend to assume a linear relationship with respect to Time. In light of this, I recommend a slightly modified model:

-model 'Time*Treatment*VoxelActivity_T1+enorm+GAD7_T1+(1+Time|Subj)' \

Regarding the issue of stronger statistical evidence than expected, it might not be surprising because cluster-based inference can be somewhat arbitrary in common practice. A more focused approach would be defining regions of interest and conducting the analysis specifically within those regions. In doing so, all inferences would be confined to those relevant regions.

Gang Chen

Thank you Gang. You are correct that we are testing for a linear effect of Time.

Can you help us understand what (1+Time|Subj) does and why it is recommended? To my understanding this produces correlated random effects, equivalent to (1 | Subject) + (Time | Subject) (according to lmer text by D. Bates ). Does that make sense given that we expect data within-subjects to be correlated with the linear slope across time varying between subjects?

The notation 1+Time|Subj) indicates that 1) both the intercept and the linear trend vary across individuals, and 2) the individual-level intercept and the linear trend are correlated. These two assumptions seem reasonable in the current context. Let us know if the corresponding analysis results make sense.

Gang Chen