================== Welcome to 3dLMEr ==================
Program for Voxelwise Linear Mixed-Effects (LME) Analysis
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Version 1.1.1, Feb 18, 2025
Author: Gang Chen (gangchen@mail.nih.gov)
SSCC/NIMH, National Institutes of Health, Bethesda MD 20892, USA
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Introduction
------
### Overview of 3dLMEr
Linear Mixed-Effects (LME) analysis follows a traditional framework that
distinguishes between two types of effects:
- Fixed effects capture population-level components.
- Random effects account for lower-level variability, such as subjects, families,
or scanning sites.
3dLMEr is an advanced and more flexible successor to 3dLME. It enhances model
specification, particularly in handling random-effects components. While 3dLME was
built on the nlme R package, 3dLMEr leverages lme4, allowing for greater flexibility.
Additionally, statistical values for main effects and interactions are approximated
using Satterthwaite’s method.
### Key Differences Between 3dLMEr and 3dLME
1. Random-effects specification:
- In 3dLMEr, random effects are fully integrated into the model formula (via `-model ...`).
- The `-ranEff` option from 3dLME is no longer needed.
- Users must explicitly specify the model structure. See this blogpost for details:
How to Specify Individual-Level Random Effects in Hierarchical Modeling
https://discuss.afni.nimh.nih.gov/t/how-to-specify-individual-level-random-effects-in-hierarchical-modeling/6462
2. Simplified effect specification:
- Labels for simple and composite effects are now part of `-gltCode` and `-glfCode`,
eliminating the need for `-gltLabel`.
3. Output format for statistical values:
- Main effects, interactions, and composite effects (generated automatically by 3dLMEr)
are stored as chi-square statistics (with 2 degrees of freedom).
- Simple effects (specified by the user) are stored as Z-statistics.
- The fixed 2 degrees of freedom for chi-square statistics simplifies interpretation,
as the Satterthwaite method produces varying degrees of freedom.
### Citing 3dLMEr
If you use 3dLMEr in your analysis, cite:
- General LME approach:
Chen, G., Saad, Z.S., Britton, J.C., Pine, D.S., Cox, R.W. (2013).
Linear Mixed-Effects Modeling Approach to FMRI Group Analysis. *NeuroImage, 73*, 176-190.
[DOI: 10.1016/j.neuroimage.2013.01.047](http://dx.doi.org/10.1016/j.neuroimage.2013.01.047)
- Test-retest reliability using trial-level effect estimates (`-TRR` option):
Chen, G., Pine, D.S., Brotman, M.A., Smith, A.R., Cox, R.W., Haller, S.P. (2021).
Trial and error: A hierarchical modeling approach to test-retest reliability.
*NeuroImage, 245*, 118647.
[DOI: 10.1016/j.neuroimage.2021.118647](https://doi.org/10.1016/j.neuroimage.2021.118647)
### Input & Output Formats
Supported input formats:
- AFNI
- NIfTI
- Surface (`niml.dset`)
- 1D text files
To match the output format to the input, append an appropriate suffix to `-prefix`
(e.g., `.nii`, `.niml.dset`, or `.1D`).
### Model Specification & Considerations
Explanatory variables:
3dLMEr supports:
- Categorical variables (e.g., between- and within-subject factors)
- Quantitative variables (e.g., age, behavioral measures)
User responsibility:
- The burden of specifying lower-level effects is on the user.
- For clarifications, refer to this FAQ: [Mixed Models FAQ]
(https://bbolker.github.io/mixedmodels-misc/glmmFAQ.html).
Handling quantitative variables:
- Declare them explicitly using `-qVars`.
- Consider centering options:
- Global centering (across all subjects)
- Within-condition/group centering (depends on research context)
- More details on centering: https://afni.nimh.nih.gov/pub/dist/doc/htmldoc/statistics/center.html
### Example Scripts
Check out example scripts below that demonstrate different data structures. If one
matches your study, use it as a template. More examples will be added over time—
contributions are welcome!
### Installation Requirements
Before running 3dLMEr, install the following R packages:
```
install.packages("lmerTest")
install.packages("phia")
install.packages("snow")
```
Alternatively, use AFNI’s installer:
```
rPkgsInstall -pkgs "lmerTest,phia,snow"
```
### Running 3dLMEr
Once your script is ready, run it in the terminal:
```
nohup tcsh -x LME.txt &
```
or, to save the output log:
```
nohup tcsh -x LME.txt > diary.txt &
```
or, to display output live while saving it:
```
nohup tcsh -x LME.txt |& tee diary.txt &
```
Saving logs allows you to review output later if issues arise.
Example 1 --- Simplest case: LME analysis for one group of subjects each of
which has three effects associated with three emotions (pos, neg and neu),
and the effects of interest are the comparisons among the three emotions
at the population level (missing data allowed). This data structure is usually
considered as one-way repeated-measures (or within-subject) ANOVA if no
missing data occurred. The LME model is typically formulated with a random
intercept in this case. With the option -bounds, values beyond [-2, 2] will
be treated as outliers and considered as missing. If you want to set a range,
choose the bounds that make sense with your input data.
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is
to use 3dGLMM with a Student's t-distribution.
-------------------------------------------------------------------------
3dLMEr -prefix LME -jobs 12 \
-mask myMask+tlrc \
-model 'emotion+(1|Subj)' \
-SS_type 3 \
-bounds -2 2 \
-gltCode mean 'emotion : 0.333*pos +0.333*neg + 0.333*neu' \
-gltCode pos 'emotion : 1*pos' \
-gltCode neg 'emotion : 1*neg' \
-gltCode neu 'emotion : 1*neu' \
-gltCode pos-neg 'emotion : 1*pos -1*neg' \
-gltCode pos-neu 'emotion : 1*pos -1*neu' \
-gltCode neg-neu 'emotion : 1*neg -1*neu' \
-gltCode em-eff1 'emotion : 0.5*pos +0.5*neg -1*neu' \
-glfCode em-eff2 'emotion : 1*pos -1*neg & 1*pos -1*neu' \
-dataTable \
Subj emotion InputFile \
s1 pos s1_pos+tlrc \
s1 neg s1_neg+tlrc \
s1 neu s1_neu+tlrc \
s2 pos s2_pos+tlrc \
s2 neg s2_neg+tlrc \
s2 pos s2_neu+tlrc \
...
s20 pos s20_pos+tlrc \
s20 neg s20_neg+tlrc \
s20 neu s20_neu+tlrc \
...
**Note:** `3dLMEr` does not explicitly output the model intercept (overall mean).
However, you can extract it using the `-gltCode` option, as shown in the script above:
-gltCode mean 'emotion : 0.333*pos +0.333*neg +0.333*neu'
Example 2 --- LME analysis for one group of subjects each of which has
three effects associated with three emotions (pos, neg and neu), and the
effects of interest are the comparisons among the three emotions at the
population level. In addition, reaction time (RT) is available per emotion
from each subject. An LME model can be formulated to include both random
intercept and random slope. Be careful about the centering issue about any
quantitative variable: you have to decide which makes more sense - global
centering or within-condition (or within-group) centering?
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is
to use 3dGLMM with a Student's t-distribution.
-------------------------------------------------------------------------
3dLMEr -prefix LME -jobs 12 \
-mask myMask+tlrc \
-model 'emotion*RT+(RT|Subj)' \
-SS_type 3 \
-bounds -2 2 \
-qVars 'RT' \
-qVarCenters 0 \
-gltCode pos 'emotion : 1*pos' \
-gltCode neg 'emotion : 1*neg' \
-gltCode neu 'emotion : 1*neu' \
-gltCode pos-neg 'emotion : 1*pos -1*neg' \
-gltCode pos-neu 'emotion : 1*pos -1*neu' \
-gltCode neg-neu 'emotion : 1*neg -1*neu' \
-gltCode em-eff1 'emotion : 0.5*pos +0.5*neg -1*neu' \
-glfCode em-eff2 'emotion : 1*pos -1*neg & 1*pos -1*neu' \
-dataTable \
Subj emotion RT InputFile \
s1 pos 23 s1_pos+tlrc \
s1 neg 34 s1_neg+tlrc \
s1 neu 28 s1_neu+tlrc \
s2 pos 31 s2_pos+tlrc \
s2 neg 22 s2_neg+tlrc \
s2 pos 29 s2_neu+tlrc \
...
s20 pos 12 s20_pos+tlrc \
s20 neg 20 s20_neg+tlrc \
s20 neu 30 s20_neu+tlrc \
...
Example 3 --- LME analysis for one group of subjects each of which has three
effects associated with three emotions (pos, neg and neu), and the effects
of interest are the comparisons among the three emotions at the population
level. As the data were acquired across 12 scanning sites, we set up an LME
model with a crossed random-effects structure, one for cross-subjects and one
for cross-sites variability.
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is
to use 3dGLMM with a Student's t-distribution.
-------------------------------------------------------------------------
3dLMEr -prefix LME -jobs 12 \
-mask myMask+tlrc \
-model 'emotion+(1|Subj)+(1|site)' \
-SS_type 3 \
-bounds -2 2 \
-gltCode pos 'emotion : 1*pos' \
-gltCode neg 'emotion : 1*neg' \
-gltCode neu 'emotion : 1*neu' \
-gltCode pos-neg 'emotion : 1*pos -1*neg' \
-gltCode pos-neu 'emotion : 1*pos -1*neu' \
-gltCode neg-neu 'emotion : 1*neg -1*neu' \
-gltCode em-eff1 'emotion : 0.5*pos +0.5*neg -1*neu' \
-glfCode em-eff2 'emotion : 1*pos -1*neg & 1*pos -1*neu' \
-dataTable \
Subj emotion site InputFile \
s1 pos site1 s1_pos+tlrc \
s1 neg site1 s1_neg+tlrc \
s1 neu site2 s1_neu+tlrc \
s2 pos site1 s2_pos+tlrc \
s2 neg site2 s2_neg+tlrc \
s2 pos site3 s2_neu+tlrc \
...
s80 pos site12 s80_pos+tlrc \
s80 neg site12 s80_neg+tlrc \
s80 neu site10 s80_neu+tlrc \
...
Example 4 --- LME analysis with a between-subject factor (group: two groups of
subjects -- control, patient), two within-subject factros (emotion: 3 levels
-- pos, neg, neu; type: 2 levels -- face, word), one quantitative variable (age).
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is to
use 3dGLMM with a Student's t-distribution.
-------------------------------------------------------------------------
3dLMEr -prefix LME -jobs 12 \
-mask myMask+tlrc \
-model 'group*emotion*type+age+(1|Subj)+(1|Subj:emotion)+(1|Subj:type)' \
-SS_type 3 \
-bounds -2 2 \
-gltCode pat.pos 'gruop : 1*patient emotion : 1*pos' \
-gltCode pat.neg 'gruop : 1*patient emotion : 1*neg' \
-gltCode ctr.pos.age 'gruop : 1*control emotion : 1*pos age :' \
-dataTable \
Subj group emotion type age InputFile \
s1 control pos face 35 s1_pos+tlrc \
s1 control neg face 35 s1_neg+tlrc \
s1 control neu face 35 s1_neu+tlrc \
s2 control pos face 23 s2_pos+tlrc \
s2 control neg face 23 s2_neg+tlrc \
s2 control pos face 23 s2_neu+tlrc \
...
s80 patient pos word 28 s80_pos+tlrc \
s80 patient neg word 28 s80_neg+tlrc \
s80 patient neu word 28 s80_neu+tlrc \
...
Example 5 --- Test-retest reliability. LME model can be adopted for test-
retest reliability analysis if trial-level effect estimates (e.g., using
option -stim_times_IM in 3dDeconvolve/3dREMLfit) are available from each
subjects. The following script demonstrates a situation where each subject
performed same two tasks across two sessions. The goal is to obtain the
test-retest reliability at the whole-brain voxel level for the contrast
between the two tasks with the test-retest reliability for the average
effect between the two tasks as a byproduct.
WARNING: numerical failures may occur, especially for a contrast between
two conditions. The failures manifest with a large portion of 0, 1 and -1
values in the output. In that case, use the program TRR to conduct
region-level test-retest reliability analysis.
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is
to use 3dGLMM with a Student's t-distribution.
-------------------------------------------------------------------------
3dLMEr -prefix output -TRR -jobs 16
-qVars 'cond'
-bounds -2 2
-model '0+sess+cond:sess+(0+sess|Subj)+(0+cond:sess|Subj)'
-dataTable @data.tbl
With many trials per condition, it is recommended that the data table
is saved as a separate file in pure text of long format with condition
(variable 'cond' in the script above) through dummy coding of -0.5 and
0.5 with the option -qVars 'cond'. Code subject and session as factor
labels with labels. Below is an example of the data table. There is no
need to add backslash at the end of each line. If sub-brick selector
is used, do NOT use gzipped files (otherwise the file reading time would
be too long) and do NOT add quotes around the square brackets [] for the
sub-brick selector.
Subj sess cond InputFile
Subj1 s1 -0.5 Subj1s1c1_trial1.nii
Subj1 s1 -0.5 Subj1s1c1_trial2.nii
...
Subj1 s1 -0.5 Subj1s1c1_trial40.nii
Subj1 s1 0.5 Subj1s1c2_trial1.nii
Subj1 s1 0.5 Subj1s1c2_trial2.nii
...
Subj1 s1 0.5 Subj1s1c2_trial40.nii
Subj1 s2 -0.5 Subj1s2c1_trial1.nii
Subj1 s2 -0.5 Subj1s2c1_trial2.nii
...
Subj1 s2 -0.5 Subj1s2c1_trial40.nii
Subj1 s2 0.5 Subj1s2c2_trial1.nii
Subj1 s2 0.5 Subj1s2c2_trial2.nii
...
Subj1 s2 0.5 Subj1s2c2_trial40.nii
...
Options in alphabetical order:
------------------------------
-bounds lb ub: This option is for outlier removal. Two numbers are expected from
the user: the lower bound (lb) and the upper bound (ub). The input data will
be confined within [lb, ub]: any values in the input data that are beyond
the bounds will be removed and treated as missing. Make sure the first number
is less than the second. The default (the absence of this option) is no
outlier removal.
**NOTE**: Using the -bounds option to remove outliers should be approached
with caution due to its arbitrariness. A more principled alternative is
to use 3dGLMM with a Student's t-distribution.
-cio: Use AFNI's C io functions, which is the default. Alternatively, -Rio
can be used.
-dataTable TABLE: List the data structure with a header as the first line.
NOTE:
1) This option has to occur last in the script; that is, no other
options are allowed thereafter. Each line should end with a backslash
except for the last line.
2) The order of the columns should not matter except that the last
column has to be the one for input files, 'InputFile'. Unlike 3dLME, the
subject column (Subj in 3dLME) does not have to be the first column;
and it does not have to include a subject ID column under some situations
Each row should contain only one input file in the table of long format
(cf. wide format) as defined in R. Input files can be in AFNI, NIfTI or
surface format. AFNI files can be specified with sub-brick selector (square
brackets [] within quotes) specified with a number or label.
3) It is fine to have variables (or columns) in the table that are
not modeled in the analysis.
4) When the table is part of the script, a backslash is needed at the end
of each line (except for the last line) to indicate the continuation to the
next line. Alternatively, one can save the context of the table as a separate
file, e.g., calling it table.txt, and then in the script specify the data
with '-dataTable @table.txt'. However, when the table is provided as a
separate file, do NOT put any quotes around the square brackets for each
sub-brick, otherwise the program would not properly read the files, unlike the
situation when quotes are required if the table is included as part of the
script. Backslash is also not needed at the end of each line, but it would
not cause any problem if present. This option of separating the table from
the script is useful: (a) when there are many input files so that the program
complains with an 'Arg list too long' error; (b) when you want to try
different models with the same dataset.
-dbgArgs: This option will enable R to save the parameters in a
file called .3dLMEr.dbg.AFNI.args in the current directory
so that debugging can be performed.
-glfCode label CODING: Specify a general linear F-style (GLF) formulation
with the weights among factor levels in which two or more null
relationships (e.g., A-B=0 and B-C=0) are involved. The symbolic
coding has to be within (single or double) quotes. For example, the
coding '-glfCode AvBvc 'Condition : 1*A -1*B & 1*A -1*C Emotion : 1*pos''
examines the main effect of Condition at the positive Emotion with
the output labeled as AvBvC. Similarly the coding '-glfCode CondByEmo'
'Condition : 1*A -1*B & 1*A -1*C Emotion : 1*pos -1*neg' looks
for the interaction between the three levels of Condition and the
two levels of Emotion and the resulting sub-brick is labeled as
'CondByEmo'.
NOTE:
1) The weights for a variable do not have to add up to 0.
2) When a quantitative variable is present, other effects are
tested at the center value of the covariate unless the covariate
value is specified as, for example, 'Group : 1*Old Age : 2', where
the Old Group is tested at the Age of 2 above the center.
3) The absence of a categorical variable in a coding means the
levels of that factor are averaged (or collapsed) for the GLF.
4) The appearance of a categorical variable has to be followed
by the linear combination of its levels.
-gltCode label weights: Specify the label and weights of interest in a general
linear t-style (GLT) formulation in which only one null relationship is
involved (cf. -glfCode). The weights should be surrounded with quotes. For
example, the specification '-gltCode AvB 'Condition : 1*A -1*B' compares A
and B with a label 'AvB' for the output sub-bricks.
-help: this help message
-IF var_name: var_name is used to specify the column name that is designated for
input files of effect estimate. The default (when this option is not invoked
is 'InputFile', in which case the column header has to be exactly as 'InputFile'
This input file for effect estimates has to be the last column.
-jobs NJOBS: On a multi-processor machine, parallel computing will speed
up the program significantly.
Choose 1 for a single-processor computer.
-mask MASK: Process voxels inside this mask only.
Default is no masking.
-model FORMULA: Specify the model structure for all the variables. The
expression FORMULA with more than one variable has to be surrounded
within (single or double) quotes. Variable names in the formula
should be consistent with the ones used in the header of -dataTable.
In the LME context the simplest model is "1+(1|Subj)" in
which the random effect from each of the two subjects in a pair is
symmetrically incorporated in the model. Each random-effects factor is
specified within parentheses per formula convention in R. Any
effects of interest and confounding variables (quantitative or
categorical variables) can be added as fixed effects without parentheses.
-prefix PREFIX: Output file name. For AFNI format, provide prefix only,
with no view+suffix needed. Filename for NIfTI format should have
.nii attached (otherwise the output would be saved in AFNI format).
-qVarCenters VALUES: Specify centering values for quantitative variables
identified under -qVars. Multiple centers are separated by
commas (,) without any other characters such as spaces and should
be surrounded within (single or double) quotes. The order of the
values should match that of the quantitative variables in -qVars.
Default (absence of option -qVarCenters) means centering on the
average of the variable across ALL subjects regardless their
grouping. If within-group centering is desirable, center the
variable YOURSELF first before the values are fed into -dataTable.
-qVars variable_list: Identify quantitative variables (or covariates) with
this option. The list with more than one variable has to be
separated with comma (,) without any other characters such as
spaces and should be surrounded within (single or double) quotes.
For example, -qVars "Age,IQ"
WARNINGS:
1) Centering a quantitative variable through -qVarsCenters is
very critical when other fixed effects are of interest.
2) Between-subjects covariates are generally acceptable.
However EXTREME caution should be taken when the groups
differ substantially in the average value of the covariate.
-R2: Enabling this option will prompt the program to provide both
conditional and marginal coefficient of determination (R^2)
values associated with the adopted model. Marginal R^2 indicates
the proportion of variance explained by the fixed effects in the
model, while conditional R^2 represents the proportion of variance
explained by the entire model, encompassing both fixed and random
effects. Two sub-bricks labeled 'R2m' and 'R2c' will be provided
in the output.
-resid PREFIX: Output file name for the residuals. For AFNI format, provide
prefix only without view+suffix. Filename for NIfTI format should
have .nii attached, while file name for surface data is expected
to end with .niml.dset. The sub-brick labeled with the '(Intercept)',
if present, should be interpreted as the effect with each factor
at the reference level (alphabetically the lowest level) for each
factor and with each quantitative covariate at the center value.
-Rio: Use R's io functions. The alternative is -cio.
-show_allowed_options: list of allowed options
-SS_type NUMBER: Specify the type for sums of squares in the F-statistics.
Three options are: sequential (1), hierarchical (2), and marginal (3).
When this option is absent (default), marginal (3) is automatically set.
Some discussion regarding their differences can be found here:
https://sscc.nimh.nih.gov/sscc/gangc/SS.html
-TRR: This option will allow the analyst to perform test-retest reliability analysis
at the whole-brain voxel level. To be able to adopt this modeling approach,
trial-level effect estimates have to be provided from each subject (e.g.,
using option -stim_times_IM in 3dDeconvolve/3dREMLfit). Currently it works
with the situation with two conditions for a group of subjects that went
two sessions. The analytical goal to assess test-retest reliability across
the two sessions for the contrast between the two conditions. Check out
Example 4 for model specification. It is possible that numerical failures
may occur for a contrast between two conditions with values of 0, 1 or -1 in
the output. Use program TRR for ROI-level test-retest reliability analysis.
-vVarCenters VALUES: Specify centering values for voxel-wise covariates
identified under -vVars. Multiple centers are separated by
commas (,) within (single or double) quotes. The order of the
values should match that of the quantitative variables in -qVars.
Default (absence of option -vVarsCenters) means centering on the
average of the variable across ALL subjects regardless their
grouping. If within-group centering is desirable, center the
variable yourself first before the files are fed under -dataTable.
-vVars variable_list: Identify voxel-wise covariates with this option.
Currently one voxel-wise covariate is allowed only. By default
mean centering is performed voxel-wise across all subjects.
Alternatively centering can be specified through a global value
under -vVarsCenters. If the voxel-wise covariates have already
been centered, set the centers at 0 with -vVarsCenters.
