Whole-brain Multiple Regression in AFNI

Hi-

My understanding is that the functions of 3dRegAna have been subsumed under other programs.

Which is the best program now for doing whole-brain multiple regression in AFNI (say, examining whether a symptom or age, or the interaction between the two, predicts an evoked response in particular clusters of brain)? Is it 3dttest++, or 3dMVM?

Thanks,

Jim Waltz

3dttest++ is computationally more efficient. If you want an interaction between two explanatory variables, you would have to create it yourself as a separate explanatory variable. In case you’re interested in the intercept, be careful about the centering issue:

https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/STATISTICS/center.html

The choice between the 3dMVM and 3dttest++ programs depends on the statistical model you have for the signal inside your data. Rather than try to explain those points of the program here, I’ll point you to the group analysis presentations from which we teach. Follow the link below and get the afni24 and afni25 PDFs.

https://afni.nimh.nih.gov/pub/dist/edu/latest/afni_handouts/

Hi- Thanks for your responses so far.

I have some questions about the output of 3dttest++ with a covariate.

  • Say I have whole-brain single-subject betas from an amplitude-modulated regressor (say, trial-by-trial prediction error estimates.
  • Say my covariate is some score from a questionnaire, like childhood trauma, or Beck Depression Inventory (one score per subject).
  • Say I have no other between-subjects or within-subjects factors.

3dttest++ spits out a mean and a t-stat for both your dataset label (where the mean prediction-error-evoked beta for the group differs from zero) and your covariate (say, trauma score).

My question is: what do those means and t-stats REPRESENT? Do they represent the main effect of each factor, controlling for the other? Do they represent the main effect of each factor, NOT controlling for the other? Is it possible to define an interaction term, in this case? Is it possible to do this with 3dMVM? Or is it impossible, because it is not really multivariate?

Thanks,

Jim Waltz

Jim,

Is your model like the following?

Y = a + b*X + e

where Y is the beta that corresponds to average predictor error for each subject. Does the average prediction error vary substantially across subjects?

what do those means and t-stats REPRESENT? Do they represent the main effect of each factor, NOT controlling for the other?

Could you show the output of the following command?

3dinfo -verb outputFrom3dttest++

Is it possible to define an interaction term, in this case? Is it possible to do this with 3dMVM? Or is it impossible, because it is not really multivariate?

How many explanatory variables do you have in the model? It’s possible to include interaction terms in both 3dttest++ and 3dMVM as I mentioned in my previous response.

Hi Gang-

What I’m trying to do is (I think) super-simple. I am only interested in how a single covariate modifies the effect of Subj.

This is the output of the 3dttest++ script I ran.


Dataset File: Pessig_Ttest_PE_New_Patients_CovBPRS+tlrc
Identifier Code: XYZ_LtPNovy7stOZkRiDL-Winw Creation Date: Mon May 8 14:53:56 2017
Template Space: TLRC
Dataset Type: Func-Bucket (-fbuc)
Byte Order: LSB_FIRST [this CPU native = LSB_FIRST]
Storage Mode: BRIK
Storage Space: 21,959,824 (22 million [mega]) bytes
Geometry String: “MATRIX(1.5,0,0,-79.5,0,1.5,0,-79.5,0,0,1.5,-65):107,127,101”
Data Axes Tilt: Plumb
Data Axes Orientation:
first (x) = Right-to-Left
second (y) = Anterior-to-Posterior
third (z) = Inferior-to-Superior [-orient RAI]
R-to-L extent: -79.500 [R] -to- 79.500 [L] -step- 1.500 mm [107 voxels]
A-to-P extent: -79.500 [A] -to- 109.500 [P] -step- 1.500 mm [127 voxels]
I-to-S extent: -65.000 [I] -to- 85.000 [S] -step- 1.500 mm [101 voxels]
Number of values stored at each pixel = 4
– At sub-brick #0 ‘Patients_mean’ datum type is float: -0.160391 to 0.160074
– At sub-brick #1 ‘Patients_Tstat’ datum type is float: -5.70935 to 7.06039
statcode = fitt; statpar = 23
– At sub-brick #2 ‘Patients_bprs_tot’ datum type is float: -0.2191 to 0.257849
– At sub-brick #3 ‘Patients_bprs_tot_Tstat’ datum type is float: -5.669 to 5.80341
statcode = fitt; statpar = 23

----- HISTORY -----
[jwaltz@admins: Mon May 8 14:53:56 2017] 3dttest++ -prefix Pessig_Ttest_PE_New_Patients_CovBPRS -mask ./mask_group+tlrc -setA Patients 12380 ‘statsGN25_new_PEonly.12380+tlrc.[All_PE#1_Coef]’ 12390 ‘statsGN25_new_PEonly.12390+tlrc.[All_PE#1_Coef]’ 12425 ‘statsGN25_new_PEonly.12425+tlrc.[All_PE#1_Coef]’ 12463 ‘statsGN25_new_PEonly.12463+tlrc.[All_PE#1_Coef]’ 12468 ‘statsGN25_new_PEonly.12468+tlrc.[All_PE#1_Coef]’ 12481 ‘statsGN25_new_PEonly.12481+tlrc.[All_PE#1_Coef]’ 12505 ‘statsGN25_new_PEonly.12505+tlrc.[All_PE#1_Coef]’ 12609 ‘statsGN25_new_PEonly.12609+tlrc.[All_PE#1_Coef]’ 12614 ‘statsGN25_new_PEonly.12614+tlrc.[All_PE#1_Coef]’ 12616 ‘statsGN25_new_PEonly.12616+tlrc.[All_PE#1_Coef]’ 12645 ‘statsGN25_new_PEonly.12645+tlrc.[All_PE#1_Coef]’ 12679 ‘statsGN25_new_PEonly.12679+tlrc.[All_PE#1_Coef]’ 12787 ‘statsGN25_new_PEonly.12787+tlrc.[All_PE#1_Coef]’ 12790 ‘statsGN25_new_PEonly.12790+tlrc.[All_PE#1_Coef]’ 12833 ‘statsGN25_new_PEonly.12833+tlrc.[All_PE#1_Coef]’ 12834 ‘statsGN25_new_PEonly.12834+tlrc.[All_PE#1_Coef]’ 12863 ‘statsGN25_new_PEonly.12863+tlrc.[All_PE#1_Coef]’ 12864 ‘statsGN25_new_PEonly.12864+tlrc.[All_PE#1_Coef]’ 12874 ‘statsGN25_new_PEonly.12874+tlrc.[All_PE#1_Coef]’ 12882 ‘statsGN25_new_PEonly.12882+tlrc.[All_PE#1_Coef]’ 12906 ‘statsGN25_new_PEonly.12906+tlrc.[All_PE#1_Coef]’ 12953 ‘statsGN25_new_PEonly.12953+tlrc.[All_PE#1_Coef]’ 12962 ‘statsGN25_new_PEonly.12962+tlrc.[All_PE#1_Coef]’ 13063 ‘statsGN25_new_PEonly.13063+tlrc.[All_PE#1_Coef]’ 13125 ‘statsGN25_new_PEonly.13125+tlrc.[All_PE#1_Coef]’ -covariates /data/brutus_data12/jwaltz/pessig/scripts/group/BPRS_Totals.txt

I tried to do the same thing with 3dMVM, but I keep getting the following error:


Error in seq.default(wd + ii, len, wd) : wrong sign in ‘by’ argument
Calls: process.MVM.opts → data.frame → cbind → seq → seq.default
Execution halted
Test_Pessig_3dMVM_AllPE_BPRS.sh: line 13: 12380: command not found

This is that script:


#! bin/sh

cd /Users/jameswaltz/Analyses/AFNIana/PessigMRI/group.results.parametric

name=‘[All_PE#1_Coef]’

3dMVM -prefix Pessig_AllPE_BPRS_Param.3dMVM -mask ./mask_group+tlrc.
-bsVars “bprstot”
-qVars “bprstot”
-jobs 8
-dataTable
Subj bprstot InputFile \
12380 150 ./statsGN25_new_PEonly.12380+tlrc.${name}
12390 325 ./statsGN25_new_PEonly.12390+tlrc.${name}
12425 200 ./statsGN25_new_PEonly.12425+tlrc.${name}
12463 100 ./statsGN25_new_PEonly.12463+tlrc.${name}
12468 125 ./statsGN25_new_PEonly.12468+tlrc.${name}
12481 100 ./statsGN25_new_PEonly.12481+tlrc.${name}
12505 100 ./statsGN25_new_PEonly.12505+tlrc.${name}
12609 150 ./statsGN25_new_PEonly.12609+tlrc.${name}
12614 100 ./statsGN25_new_PEonly.12614+tlrc.${name}
12616 400 ./statsGN25_new_PEonly.12616+tlrc.${name}
12645 100 ./statsGN25_new_PEonly.12645+tlrc.${name}
12679 150 ./statsGN25_new_PEonly.12679+tlrc.${name}
12787 100 ./statsGN25_new_PEonly.12787+tlrc.${name}
12790 100 ./statsGN25_new_PEonly.12790+tlrc.${name}
12833 175 ./statsGN25_new_PEonly.12833+tlrc.${name}
12834 150 ./statsGN25_new_PEonly.12834+tlrc.${name}
12863 125 ./statsGN25_new_PEonly.12863+tlrc.${name}
12864 150 ./statsGN25_new_PEonly.12864+tlrc.${name}
12874 225 ./statsGN25_new_PEonly.12874+tlrc.${name}
12882 125 ./statsGN25_new_PEonly.12882+tlrc.${name}
12906 275 ./statsGN25_new_PEonly.12906+tlrc.${name}
12953 150 ./statsGN25_new_PEonly.12953+tlrc.${name}
12962 300 ./statsGN25_new_PEonly.12962+tlrc.${name}
13063 100 ./statsGN25_new_PEonly.13063+tlrc.${name}
13125 275 ./statsGN25_new_PEonly.13125+tlrc.${name}

Since your model is

Y = a + b*X + e

the sub-bricks #0 and #1 in your 3dttest++ output are for the intercept ‘a’ in the model above, which is the group average effect that is associated with the average value of X (‘bprstot’ in your case), and the sub-bricks #2 and #3 in your 3dttest++ output are for the effect of ‘bprstot’.

For your 3dMVM script, check if there is some extra space after the backslash on the following line:

Subj bprstot InputFile \

Hi-

Thank you. Your response answer both questions of mine. I got the 3dMVM script to run.

Last question: do the first 2 subbriks characterize the effect of subject, with the effect of BPRStot REGRESSED OUT?

Thanks,

Jim Waltz

do the first 2 subbriks characterize the effect of subject, with the effect of BPRStot REGRESSED OUT?

“Regressing out” is not an accurate description. Instead, as I mentioned in my previous post, the sub-bricks #0 and #1 in your 3dttest++ output are for the intercept ‘a’ in the model above, which is the group average effect that is associated with the average value of X (‘bprstot’ in your case).

Hi Gang-

So is it somewhat accurate to say that, if your regression model includes “brain” and BPRStot (which is a psychosis score), and your “brain” signals already reflect correspondence to a parametric regressor (reward prediction error, in this case), that the “BPRStot” coefficient that comes out of either 3dttest++ with a covariate or 3dMVM with a covariate reflects something like the MODULATORY EFFECT of BPRStot on the brain response to prediction errors (and thus a kind of, sort of interaction)?

Jim Waltz

Jim,

Word description tends to be elusive. Personally I prefer to use the model to demonstrate the specifics. If your model is

Y = a + b*X + e

where Y is the brain response, and X is “BPRStot”.

“brain” signals already reflect correspondence to a parametric regressor (reward prediction error, in this case)

In a modulation analysis, there are two betas in the output from, for example, 3dDeconvolve: which one are you referring to, the first or second?

The second (parametric) one.

Jim,

Read the following discussion:

https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/STATISTICS/center.html#centering-with-more-than-one-group-of-subjects

and see if it resolves the confusion.