Skip to contents

Introduction: Generating AFNI Models with afnireg

The afnireg package provides an R interface to AFNI’s 3dDeconvolve program for fMRI analysis. It allows you to specify complex experimental designs and HRF models in R using the fmridesign framework and execute them using AFNI’s optimized estimation routines.

This vignette demonstrates how to construct a full fMRI regression model within R and then translate it into the syntax and associated files required for AFNI’s 3dDeconvolve program. This allows you to use R’s flexible formula interface and HRF specification tools while utilizing AFNI for the final voxelwise model estimation.

The primary function for this translation is afni_lm().

A Simple Model Example

Let’s start with the simple single-run, single-factor design from the “Building Event Models” vignette.

TR <- 2
cond <- c("face", "scene", "tool", "object")
NSTIM <- length(cond) * 4

set.seed(123)
simple_design <- data.frame(
  stim = factor(sample(rep(cond, 4))),
  ISI = sample(10:20, NSTIM, replace = TRUE), # Using wider ISI from previous fixes
  run = rep(1, NSTIM),
  trial = factor(1:NSTIM)
)

simple_design$onset <- cumsum(c(0, simple_design$ISI[-NSTIM] + 2))
sframe <- sampling_frame(blocklens = 140, TR = TR)

# Define a contrast
con1 <- pair_contrast(~ stim == "face", ~stim == "scene", name = "face_scene")

# Create the event model
emodel <- event_model(onset ~ hrf(stim, contrasts = list(con1)), 
                     data = simple_design, 
                     block = ~ run, 
                     sampling_frame = sframe)

# Create the baseline model (B-spline drift)
bmodel <- baseline_model(basis = "bs", degree = 5, sframe = sframe)

# Combine into a full fmri_model
fmodel <- fmri_model(emodel, bmodel)
#> Error in eval(assertion, env): argument "dataset" is missing, with no default
# Define the dataset (using placeholder file names)
# Note: This requires actual neuroimaging files to run
dset <- fmri_dataset(scans = "scan01.nii.gz",
                     mask = "mask.nii.gz",
                     TR = TR,
                     run_length = 140,
                     event_table = simple_design,
                     base_path = ".") # Assuming files are in the working directory

Translating the Model with afni_lm

The afni_lm() function takes the fmri_model and fmri_dataset objects and generates an AFNI model specification.

# Note: This requires a valid fmri_dataset object (dset)
alm_spec <- afni_lm(fmodel, dset)
print(alm_spec)

This creates an afni_lm_spec object which contains:

  • The original fmri_model and fmri_dataset.
  • The generated 3dDeconvolve command string (alm_spec$cmd$cmd).
  • Parsed command line arguments (alm_spec$cmd$cmdlines).
  • Information about stimulus timing files, GLTs, etc.

Key Translations:

  • HRF Convolution: The package performs the HRF convolution specified in event_model (e.g., hrf(stim) uses HRF_SPMG1) and passes the resulting regressors to 3dDeconvolve via -stim_file arguments. Each condition level gets its own stimulus file (e.g., stim[face]_reg.1D).
  • Baseline Model: Regressors from the baseline_model (here, the B-spline basis columns) are included as nuisance regressors via the -ortvec argument.
  • Contrasts: R contrasts (like con1 defined with pair_contrast) are automatically converted into AFNI-style general linear tests (GLTs) using the -gltsym and -glt_label arguments. A separate text file is generated for each GLT (e.g., GLT_face_scene.txt).

Handling AFNI’s -polort

By default, afni_lm sets polort = -1 in the 3dDeconvolve command. This disables AFNI’s built-in polynomial baseline drift modeling. This is the desired behavior when you provide your own drift regressors via fmrireg’s baseline_model (like our B-spline model bmodel), as including both would lead to redundant and collinear regressors.

If you didn’t include drift terms in your baseline_model and wanted to use AFNI’s polynomial drift correction, you could specify polort in the afni_lm call:

# Example: No drift in bmodel, use AFNI polort
bmodel_no_drift <- baseline_model(basis = "constant", sframe = sframe) # Only intercept
fmodel_no_drift <- fmri_model(emodel, bmodel_no_drift)

# Specify polort=2 for quadratic drift
alm_spec_polort <- afni_lm(fmodel_no_drift, dset, polort = 2)
print(alm_spec_polort)

However, explicitly defining drift terms in R’s baseline_model is generally recommended for consistency.

Passing Custom 3dDeconvolve Options

You can pass additional options directly to the 3dDeconvolve command using the options argument in afni_lm. This accepts a list where element names correspond to 3dDeconvolve flags (without the leading -) and values are the arguments for those flags.

# Example: Request statistic and coefficient buckets, disable t-stat output
alm_spec_opts <- afni_lm(fmodel, dset, 
                         options = list(bucket = "stats_afni", 
                                        cbucket = "coefs_afni",
                                        tout = FALSE))
print(alm_spec_opts)

Using AFNI-Native HRF Basis Functions

A powerful feature is the ability to specify basis functions within event_model that directly map to AFNI’s built-in HRF models (like BLOCK, TENT, SPMG1, etc.). This avoids pre-convolution in R and lets 3dDeconvolve handle the HRF modeling using its optimized routines.

Use the afni_hrf() specifier within the event_model formula instead of hrf().

# Example using AFNI's BLOCK function for convolution
# Note: Using afni_hrf requires AFNI to be installed for 3dDeconvolve to run
emodel_afni_block <- event_model(
  onset ~ afni_hrf(stim, basis = "BLOCK(2,1)"), # BLOCK(duration, penalty)
  data = simple_design, 
  block = ~ run, 
  sampling_frame = sframe
)

fmodel_afni <- fmri_model(emodel_afni_block, bmodel)
alm_spec_afni_native <- afni_lm(fmodel_afni, dset)
print(alm_spec_afni_native)

# Check the command: It should now use -stim_times instead of -stim_file
# and specify 'BLOCK(2,1)' as the basis function.
#> Skipping AFNI-native HRF example as AFNI availability check is FALSE.

Refer to the afni_hrf() documentation and the 3dDeconvolve help for details on available AFNI basis functions and their parameters.

Handling Nuisance Regressors and Censoring

  • Nuisance Regressors: If your baseline_model includes nuisance regressors loaded from files (e.g., motion parameters), afni_lm will include them via the -ortvec argument, similar to how drift terms are handled.
  • Censoring: You can provide a censoring vector (or list of vectors per run) to the censor argument of afni_lm. This will generate a censor.1D file and add the -censor censor.1D option to the 3dDeconvolve command.
# Example: Censor the first 5 scans
censor_vec <- rep(1, 140) # Start with all scans included (1)
censor_vec[1:5] <- 0       # Censor first 5 scans (0)

alm_spec_censored <- afni_lm(fmodel, dset, censor = censor_vec)
print(alm_spec_censored)

Generating Files and Running the Model

Once you have the afni_lm_spec object, you can generate the necessary files and optionally execute the 3dDeconvolve command using the run() method.

# This command will:
# 1. Create the output directory 'glm_afni_output'.
# 2. Change into that directory.
# 3. Write all necessary files:
#    - Stimulus files (.1D) for each regressor if using standard hrf()
#    - Stimulus timing files (.1D) if using afni_hrf()
#    - Nuisance regressor files (.1D)
#    - GLT files (.txt)
#    - Censor file (censor.1D) if specified
#    - The 3ddeconvolve.sh script containing the full command
# 4. Execute the 3ddeconvolve.sh script.
# 5. Change back to the original directory.

run(alm_spec_opts, outdir = "glm_afni_output")

# To only generate the script and files without running:
# run(alm_spec_opts, outdir = "glm_afni_output", execute = FALSE)

This workflow allows you to define complex models in R using the afnireg package and then leverage AFNI’s 3dDeconvolve for efficient estimation.