I am wondering if a voxel-wise linear mixed-effects model via 3dLME is appropriate for my dissertation dataset.
I want to examine how response to uncertain threat might vary as a function of trauma symptom severity (PCL-5). I am using an uncertain picture task that examines valence (neutral, negative) and predictability (predictable, unpredictable), for a total of 4 within-subject conditions (predictable-neutral, unpredictable-negative, etc.). I also want to control for age and sex. I think that I could use 3dLME with valence (2) and predictability (2) as within-subjects factors, age and sex as between-subjects factors, and PCL-5 scores as a within-subject quantitative variable.
Does this seem appropriate? I read through the Chen 2013 paper but I still can’t quite figure out if this fits or not. The other option I considered was a hierarchical linear regression analysis.
Since you have a within-subject covariate (PCL-5), 3dLME is the way to go. One subtlety is whether PCL-5 differs between the levels of all the factors (including both between- and within-factors).
I am beginning to write the script for modeling Valence and Predictability as within-subjects factors, Age and Sex as between-subjects factors, and PCL-5 as a within-subjects quantitative variable (or would PCL-5 be a between-subjects factor if it is a measure of PTSD symptoms that is only measured once?). Below is a sketch of that script. The examples on https://afni.nimh.nih.gov/pub/dist/doc/program_help/3dLME.html don’t quite align with what I want to look at, so I am not sure what -model should look like, if I need to center my qVars, and if I need to spell out each glt. Any thoughts or advice?
3dLME -prefix unc_cl_test -jobs 1
-model “ValencePredictPCL5+Age+Sex”
-qVars “Age,Sex,PCL5”
-qVarCenters “”
-ranEff ‘~1’
-SS_type 3
-num_glt
-gltLabel 1
-dataTable
Subj Valence Predict Age Sex PCL5 InputFile
102_day1 neg p 49.3 m 14 /102_day1/unc_bl_results_peak/Pneg_block+tlrc
102_day1 neu p 49.3 m 14 /102_day1/unc_bl_results_peak/Pneu_block+tlrc
102_day1 neg u 49.3 m 14 /102_day1/unc_bl_results_peak/Uneg_block+tlrc
102_day1 neu u 49.3 m 14 /102_day1/unc_bl_results_peak/Uneg_block+tlrc
104_day1 neg p 24.8 f 35 /104_day1/unc_bl_results_peak/Pneg_block+tlrc
104_day1 neu p 24.8 f 35 /104_day1/unc_bl_results_peak/Pneu_block+tlrc
104_day1 neg u 24.8 f 35 /104_day1/unc_bl_results_peak/Uneg_block+tlrc
104_day1 neu u 24.8 f 35 /104_day1/unc_bl_results_peak/Uneg_block+tlrc
109_day1 neg p 21.7 f 46 /109_day1/unc_bl_results_peak/Pneg_block+tlrc
109_day1 neu p 21.7 f 46 /109_day1/unc_bl_results_peak/Pneu_block+tlrc
109_day1 neg u 21.7 f 46 /109_day1/unc_bl_results_peak/Uneg_block+tlrc
109_day1 neu u 21.7 f 46 /109_day1/unc_bl_results_peak/Uneg_block+tlrc
111_day1 neg p 25.6 m 49 /111_day1/unc_bl_results_peak/Pneg_block+tlrc
111_day1 neu p 25.6 m 49 /111_day1/unc_bl_results_peak/Pneu_block+tlrc
111_day1 neg u 25.6 m 49 /111_day1/unc_bl_results_peak/Uneg_block+tlrc
111_day1 neu u 25.6 m 49 /111_day1/unc_bl_results_peak/Uneg_block+tlrc
119_day1 neg p 44.6 f 26 /119_day1/unc_bl_results_peak/Pneg_block+tlrc
119_day1 neu p 44.6 f 26 /119_day1/unc_bl_results_peak/Pneu_block+tlrc
119_day1 neg u 44.6 f 26 /119_day1/unc_bl_results_peak/Uneg_block+tlrc
119_day1 neu u 44.6 f 26 /119_day1/unc_bl_results_peak/Uneg_block+tlrc
Depending on your research hypothesis, the above model might be fine. Or alternatively, you may try
-model “SexValencePredict*PCL5+Age”
-qVars “Age,Sex,PCL5”
Sex is not a quantitative variable, so remove from the above line.
-qVarCenters “”
Does the average PCL5 values varies between the two valences and between the two Predict levels? If so, you may center within each level separately, and then set
-qVarCenters “0,0”
If not, remove the option -qVarCenters and let the program center around the overall average of PCL5.
-ranEff ‘~1’
You want to account the random-slope effect of Age. So, change it to
I saw in a different post that 3dLME will run the main effects and interactions, but the user has to specify all of the follow-up tests. Does that include looking at direct comparisons of the main effects? For instance, would the following follow-ups via glt be appropriate if I wanted to know look at all of the main effects and interactions with and without trauma severity (PCL5; see script below)?
Also, I am a little confused about glf. The description in the help page of glfCode makes it seem like it is typically used for testing the main effects and interaction -but isn’t that what 3dLME is already doing?
Please let me know if you have any questions. Thanks again for your time.
would the following follow-ups via glt be appropriate if I wanted to know look at all of the main effects
and interactions with and without trauma severity (PCL5; see script below)?
The question does not make sense to me. Once you incorporate a quantitative variable (e.g., PCL5 in your case), you can make inferences by pretending the absence of its effect. Instead, you can only control the effect of the quantitative variable at a particular value. For example, the first 7 post-hoc tests in your script are conducted with PCL5 fixed at the center value. I’m not so sure whether the next 7 post-hoc tests are really what you wanted, but they give you the difference of the PCL5 effects between the two levels of a factor.
I am a little confused about glf. The description in the help page of glfCode makes it seem like it is
typically used for testing the main effects and interaction -but isn’t that what 3dLME is already doing?
They are meant to provide an F-test for a subset of a main effect or interaction. For example, the typical main effect for a factor with 4 levels tests whether there are differences among all the 4 levels, but you can use -glf specification to test whether there are differences across, for example, the first 3 levels.
Thank you very much for that information. I think I was able to figure out what I needed to include in my final script. I have been trying to run a test version of my script (just on 3 subjects and without the glfs) for awhile and finally got it to read the input files. However, now I am receiving this error: Model Test Failed. Here is the output. Any thoughts about how to proceed? I also put the script below as well.
Package nlme loaded successfully!
Package phia loaded successfully!
++++++++++++++++++++++++++++++++++++++++++++++++++++
***** Summary information of data structure *****
3 subjects : 102_day1 104_day1 109_day1
12 response values
2 levels for factor Valence : neg neu
2 levels for factor Predict : p u
12 centered values for numeric variable Age : 17.36667 17.36667 17.36667 17.36667 -7.133333 -7.133333 -7.133333 -7.133333 -10.23333 -10.23333 -10.23333 -10.23333
2 levels for factor Sex : 0 1
12 centered values for numeric variable PCL5 : -17.66667 -17.66667 -17.66667 -17.66667 3.333333 3.333333 3.333333 3.333333 14.33333 14.33333 14.33333 14.33333
14 post hoc tests
Contingency tables of subject distributions among the categorical variables:
Tabulation of subjects against all categorical variables
Subj vs Valence:
neg neu
102_day1 2 2
104_day1 2 2
109_day1 2 2
Subj vs Predict:
p u
102_day1 2 2
104_day1 2 2
109_day1 2 2
Subj vs Sex:
0 1
102_day1 4 0
104_day1 0 4
109_day1 0 4
***** End of data structure information *****
++++++++++++++++++++++++++++++++++++++++++++++++++++
Reading input files now...
Reading input files: Done!
If the program hangs here for more than, for example, half an hour,
kill the process because the model specification or the GLT coding
is likely inappropriate.
~~~~~~~~~~~~~~~~~~~ Model test failed ~~~~~~~~~~~~~~~~~~~
Possible reasons:
0) Make sure that R packages nlme and phia have been installed. See the 3dLME
help documentation for more details.
1) Inappropriate model specification with options -model, or -qVars.
2) In correct specifications in general linear test coding with -gltCode.
3) Mistakes in data table. Check the data structure shown above, and verify
whether there are any inconsistencies.
4) Inconsistent variable names which are case sensitive. For example, factor
named Group in model specification and then listed as group in the table hader
would cause grief for 3dLME.
#!/bin/tcsh
module load pkgsrc/2018Q3
module load pkgsrc/gcc-base
module load /sharedapps/LS/psych_imaging/modulefiles/afni/08.28.2017
3dLME -prefix unc_cl_test -jobs 6 \
-model "Valence*Predict*PCL5+Age+Sex" \
-qVars "Age,PCL5" \
-ranEff "~1+Age" \
-SS_type 3 \
-num_glt 14 \
-gltLabel 1 'neg-neu' -gltCode 1 'Valence : 1*neg -1*neu' \
-gltLabel 2 'u-p' -gltCode 2 'Predict : 1*u -1*p' \
-gltLabel 3 'Uneg-Uneu' -gltCode 3 'Valence : 1*neg -1*neu Predict : 1*u' \
-gltLabel 4 'Pneg-Pneu' -gltCode 4 'Valence : 1*neg -1*neu Predict : 1*p' \
-gltLabel 5 'Uneg-Pneg' -gltCode 5 'Predict : 1*u -1*p Valence : 1*neg' \
-gltLabel 6 'Uneu-Pneu' -gltCode 6 'Predict : 1*u -1*p Valence : 1*neu' \
-gltLabel 7 'Uneg-Pneu' -gltCode 7 'Predict : 1*u -1*p Valence : 1*neg -1*neu' \
-gltLabel 8 'neg-neu_PCL5' -gltCode 8 'Valence : 1*neg -1*neu PCL5 :' \
-gltLabel 9 'u-p_PCL5' -gltCode 9 'Predict : 1*u -1*p PCL5 :' \
-gltLabel 10 'Uneg-Uneu_PCL5' -gltCode 10 'Valence : 1*neg -1*neu Predict : 1*u PCL5 :' \
-gltLabel 11 'Pneg-Pneu_PCL5' -gltCode 11 'Valence : 1*neg -1*neu Predict : 1*p PCL5 :' \
-gltLabel 12 'Uneg-Pneg_PCL5' -gltCode 12 'Predict : 1*u -1*p Valence : 1*neg PCL5 :' \
-gltLabel 13 'Uneu-Pneu_PCL5' -gltCode 13 'Predict : 1*u -1*p Valence : 1*neu PCL5 :' \
-gltLabel 14 'Uneg-Pneu_PCL5' -gltCode 14 'Predict : 1*u -1*p Valence : 1*neg -1*neu PCL5 :' \
-dataTable \
Subj Valence Predict Age Sex PCL5 InputFile \
102_day1 neg p 49.3 0 14 102_day1/unc_trl_results/peak/Pneg_cl+tlrc \
102_day1 neu p 49.3 0 14 102_day1/unc_trl_results/peak/Pneu_cl+tlrc \
102_day1 neg u 49.3 0 14 102_day1/unc_trl_results/peak/Uneg_cl+tlrc \
102_day1 neu u 49.3 0 14 102_day1/unc_trl_results/peak/Uneu_cl+tlrc \
104_day1 neg p 24.8 1 35 104_day1/unc_trl_results/peak/Pneg_cl+tlrc \
104_day1 neu p 24.8 1 35 104_day1/unc_trl_results/peak/Pneu_cl+tlrc \
104_day1 neg u 24.8 1 35 104_day1/unc_trl_results/peak/Uneg_cl+tlrc \
104_day1 neu u 24.8 1 35 104_day1/unc_trl_results/peak/Uneu_cl+tlrc \
109_day1 neg p 21.7 1 46 109_day1/unc_trl_results/peak/Pneg_cl+tlrc \
109_day1 neu p 21.7 1 46 109_day1/unc_trl_results/peak/Pneu_cl+tlrc \
109_day1 neg u 21.7 1 46 109_day1/unc_trl_results/peak/Uneg_cl+tlrc \
109_day1 neu u 21.7 1 46 109_day1/unc_trl_results/peak/Uneu_cl+tlrc
I have been trying to run a test version of my script (just on 3 subjects and without the glfs) for awhile
The failure is most likely due to the fact you only fed the model with 3 subjects: not enough to estimate all the effects with such a meager amount of data.
I created a table text file “table_unc_cl.txt” that contains all of my usable subjects. I called upon it in my script but there must be something wrong with the table text file - I keep receiving the error: “The content under -dataTable is not rectangular ! 5574 14”. I tried to use file_tool to find an error in the table but that command is not found. I created the text file using excel and saved it as a Text (Tab Delimited). Could that be the problem? Is there a better way to export the text file so that it is readable? Below is my update script and the table (table_unc_cl.txt).
If there is missing data, how should I remedy that in the table?
Remove those rows with missing data from the table unless you could impute the missing values somehow.
The
National Institute of Mental Health (NIMH) is part of the National Institutes of
Health (NIH), a component of the U.S. Department of Health and Human
Services.