Contrasts in fmridesign

Introduction to Contrasts

Contrasts specify linear combinations of GLM parameters (β) to test hypotheses about conditions or trends. In fmridesign, define contrasts inside hrf() calls using helper functions, then extract weights with contrast_weights() (and F-sets with Fcontrasts()).

Types of Contrasts

T-contrasts use a single contrast vector for directional questions (e.g., A > B). F-contrasts use a matrix to test sets of effects (e.g., any difference among levels), returning omnibus statistics.

Defining Contrasts in Event Models

Contrasts can be specified directly within the event_model() formula using the contrasts argument in hrf() terms:

# Create a simple two-condition experiment
set.seed(123)
sframe <- sampling_frame(blocklens = 200, TR = 2)

# Generate events
onsets <- sort(runif(40, 0, 350))
conditions <- rep(c("left", "right"), 20)

# Define contrasts within the model
emodel_with_contrast <- event_model(
  onset ~ hrf(hand, contrasts = pair_contrast(~ hand == "left", ~ hand == "right", name = "left_vs_right")),
  data = data.frame(
    onset = onsets,
    hand = factor(conditions),
    block = factor(rep(1, 40))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# Extract contrast weights
contrast_weights(emodel_with_contrast)
#> $`hand#left_vs_right`
#> contrast: left_vs_right 
#>  A:  ~hand == "left" 
#>  B:  ~hand == "right" 
#>  term:  hand 
#>  weights:  
#>            left_vs_right
#> hand.left              1
#> hand.right            -1
#>  conditions:  hand_hand.left hand_hand.right

Quick Validation

Validate all attached contrasts for a model in one call:

res <- validate_contrasts(emodel_with_contrast)
# 'full_rank' is only meaningful for F-contrasts; drop it when all are t-contrasts
if (!any(res$type == "F")) res$full_rank <- NULL
res
#>                 name type estimable sum_to_zero orthogonal_to_intercept
#> 2          hand#hand    t      TRUE        TRUE                    TRUE
#> 1 hand#left_vs_right    t      TRUE        TRUE                    TRUE
#>   nonzero_weights
#> 2               2
#> 1               2

Advanced Contrast Specifications

Pairwise Contrasts

For designs with multiple levels, you can specify all pairwise comparisons:

# Three-condition experiment
set.seed(456)
conditions_3 <- rep(c("easy", "medium", "hard"), each = 15)
onsets_3 <- sort(runif(45, 0, 350))

emodel_pairwise <- event_model(
  onset ~ hrf(difficulty,
              contrasts = pairwise_contrasts(c("easy","medium","hard"), facname = "difficulty", name_prefix = "pair")),
  data = data.frame(
    onset = onsets_3,
    difficulty = factor(conditions_3, levels = c("easy", "medium", "hard")),
    block = factor(rep(1, 45))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# View all contrasts
contrasts_list <- contrast_weights(emodel_pairwise)
names(contrasts_list)
#> [1] "difficulty#pair_easy_medium" "difficulty#pair_easy_hard"  
#> [3] "difficulty#pair_medium_hard"

Polynomial Contrasts

Test for linear, quadratic, or higher‑order trends across ordered conditions. Note: poly_contrast() with degree = k returns an F‑contrast matrix with k orthogonal columns (linear, quadratic, …, k‑th), so you usually only need a single call rather than separate one‑by‑one contrasts.

# Parametric design with 4 levels
set.seed(789)
intensity_levels <- rep(1:4, each = 12)
onsets_param <- sort(runif(48, 0, 350))

# Create a single polynomial contrast up to cubic (3 columns: linear/quadratic/cubic)
emodel_polynomial <- event_model(
  onset ~ hrf(intensity,
              contrasts = list(
                poly_trend = poly_contrast(~ intensity, name = "poly_trend", degree = 3)
              )),
  data = data.frame(
    onset = onsets_param,
    intensity = factor(intensity_levels),
    block = factor(rep(1, 48))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# Extract polynomial contrast weights (3 columns = linear, quadratic, cubic)
poly_contrasts <- contrast_weights(emodel_polynomial)
names(poly_contrasts)
#> [1] "intensity#poly_trend"
pc <- poly_contrasts[[1]]
round(head(pc$weights, 6), 3)
#>             poly_trend_1 poly_trend_2 poly_trend_3
#> intensity.1       -0.671          0.5       -0.224
#> intensity.2       -0.224         -0.5        0.671
#> intensity.3        0.224         -0.5       -0.671
#> intensity.4        0.671          0.5        0.224

# If you want a specific component, select a column (e.g., linear = column 1)
linear_weights <- pc$weights[, 1, drop = FALSE]
head(linear_weights)
#>             poly_trend_1
#> intensity.1   -0.6708204
#> intensity.2   -0.2236068
#> intensity.3    0.2236068
#> intensity.4    0.6708204

Factorial Design Contrasts

Main Effects and Interactions

For factorial designs, specify contrasts for main effects and interactions:

# 2x2 factorial design
set.seed(234)
n_trials <- 60
factor_A <- rep(c("A1", "A2"), each = 30)
factor_B <- rep(rep(c("B1", "B2"), each = 15), 2)
factorial_onsets <- sort(runif(n_trials, 0, 350))

emodel_factorial <- event_model(
  onset ~ hrf(A, B, contrasts = contrast_set(
    oneway_contrast(~ A, name = "main_A"),
    oneway_contrast(~ B, name = "main_B"),
    interaction_contrast(~ A * B, name = "A_by_B")
  )),
  data = data.frame(
    onset = factorial_onsets,
    A = factor(factor_A),
    B = factor(factor_B),
    block = factor(rep(1, n_trials))
  ),
  block = ~ block,
  sampling_frame = sframe
)

interaction_contrasts <- contrast_weights(emodel_factorial)
lapply(interaction_contrasts, function(x) round(x$weights, 3))
#> $`A_B#main_A`
#>           main_A_1
#> A.A1_B.B1       -1
#> A.A2_B.B1        1
#> A.A1_B.B2       -1
#> A.A2_B.B2        1
#> 
#> $`A_B#main_B`
#>           main_B_1
#> A.A1_B.B1       -1
#> A.A2_B.B1       -1
#> A.A1_B.B2        1
#> A.A2_B.B2        1
#> 
#> $`A_B#A_by_B`
#>       A_by_B_1
#> A1_B1        1
#> A2_B1       -1
#> A1_B2       -1
#> A2_B2        1

Parametric Modulators: Do You Need a Contrast?

For a simple parametric modulator (e.g., hrf(RT) with a single‑basis HRF), the effect is captured by a single regressor. In standard GLM software, the t‑test on that coefficient directly tests the parametric effect—no explicit contrast is required. Contrasts are useful if you need to combine multiple columns (e.g., multi‑basis HRF, or a set of RT‑by‑condition regressors) into a single hypothesis.

# Design with conditions and RT modulation
set.seed(567)
n_events <- 50
pm_conditions <- rep(c("congruent", "incongruent"), each = 25)
pm_onsets <- sort(runif(n_events, 0, 350))
pm_RT <- rnorm(n_events, mean = ifelse(pm_conditions == "congruent", 0.5, 0.7), sd = 0.1)

emodel_parametric <- event_model(
  onset ~ hrf(condition,
              contrasts = pair_contrast(~ condition == "incongruent", ~ condition == "congruent", name = "incongruent_gt_congruent")) +
          hrf(RT),
  data = data.frame(
    onset = pm_onsets,
    condition = factor(pm_conditions),
    RT = scale(pm_RT)[,1],
    block = factor(rep(1, n_events))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# Categorical contrast (condition) is attached; the RT parametric term is a single regressor.
# Identify its column(s) for reporting; the model t‑stat on this column tests the parametric effect.
dm <- design_matrix(emodel_parametric)
idx <- term_indices(dm)
idx[["RT"]]
#> [1] 3
colnames(dm)[idx[["RT"]]]
#> [1] "RT_RT"

F-contrasts for Omnibus Tests

F-contrasts test multiple parameters simultaneously:

# Four-condition design for omnibus test
set.seed(890)
conditions_4 <- rep(c("A", "B", "C", "D"), each = 12)
onsets_4 <- sort(runif(48, 0, 350))

# Using oneway_contrast for main effect
emodel_omnibus <- event_model(
  onset ~ hrf(condition, 
              contrasts = oneway_contrast(~ condition, name = "main_effect")),
  data = data.frame(
    onset = onsets_4,
    condition = factor(conditions_4),
    block = factor(rep(1, 48))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# Extract F-contrast
f_contrasts <- Fcontrasts(emodel_omnibus)
print(f_contrasts)
#> $`condition#condition`
#>                       c1 c2 c3
#> condition_condition.A  1  0  0
#> condition_condition.B  0  1  0
#> condition_condition.C  0  0  1
#> condition_condition.D -1 -1 -1
#> attr(,"term_indices")
#> [1] 1 2 3 4

Working with Contrast Weights

Extracting and Manipulating Weights

# Create a model
simple_model <- event_model(
  onset ~ hrf(stim),
  data = data.frame(
    onset = c(10, 30, 50, 70, 90),
    stim = factor(c("A", "B", "A", "B", "A")),
    block = factor(rep(1, 5))
  ),
  block = ~ block,
  sampling_frame = sampling_frame(60, TR = 2)
)

# Manually create contrast weights
design_mat <- design_matrix(simple_model)
n_cols <- ncol(design_mat)

# Create a custom contrast vector
custom_contrast <- rep(0, n_cols)
# Find columns for condition A and B
# Note: column names use dot notation for factor levels (e.g., "stim.A")
col_names <- colnames(design_mat)
# Match by suffix to handle term-tag prefixes like "stim_stim.A"
a_cols <- grep("stim\\.A$", col_names)
b_cols <- grep("stim\\.B$", col_names)

# A > B contrast
custom_contrast[a_cols] <- 1/length(a_cols)
custom_contrast[b_cols] <- -1/length(b_cols)

print(custom_contrast)
#> [1]  1 -1

Contrast Validation

Validate contrasts once the model (and design matrix) is constructed. You can validate built-in contrast specs or your own custom vectors using validate_contrasts():

# Validate the custom contrast against the design implied by the model
validate_contrasts(simple_model, weights = custom_contrast)
#>       name type estimable sum_to_zero orthogonal_to_intercept full_rank
#> 1 contrast    t      TRUE        TRUE                    TRUE        NA
#>   nonzero_weights
#> 1               2

Contrasts for Multi-Run Designs

When working with multiple runs, contrasts can be specified to test within-run or across-run effects:

# Two-run experiment with potential run differences
set.seed(345)
run1_onsets <- sort(runif(20, 0, 200))
run2_onsets <- sort(runif(20, 0, 200))
all_onsets <- c(run1_onsets, run2_onsets)
all_conditions <- rep(c("stim", "control"), 20)
all_blocks <- rep(1:2, each = 20)

emodel_multirun <- event_model(
  onset ~ hrf(condition, block,
              contrasts = list(
                overall   = pair_contrast(~ condition == "stim", ~ condition == "control", name = "overall"),
                run1_only = pair_contrast(~ condition == "stim", ~ condition == "control", name = "run1", where = ~ block == "1"),
                run2_only = pair_contrast(~ condition == "stim", ~ condition == "control", name = "run2", where = ~ block == "2")
              )),
  data = data.frame(
    onset = all_onsets,
    condition = factor(all_conditions),
    block = factor(all_blocks)
  ),
  block = ~ block,
  sampling_frame = sampling_frame(c(120, 120), TR = 2)
)

multirun_contrasts <- contrast_weights(emodel_multirun)
names(multirun_contrasts)
#> [1] "condition_block#overall" "condition_block#run1"   
#> [3] "condition_block#run2"

Contrast Specification Best Practices

Plan contrasts a priori, keep them simple and interpretable, and prefer orthogonal sets when possible.

2. Scaling and Normalization

# Properly scaled contrasts sum to zero
create_scaled_contrast <- function(n_positive, n_negative) {
  c(rep(1/n_positive, n_positive), rep(-1/n_negative, n_negative))
}

# Example: 3 vs 2 conditions
scaled_contrast <- create_scaled_contrast(3, 2)
print(scaled_contrast)
#> [1]  0.3333333  0.3333333  0.3333333 -0.5000000 -0.5000000
print(sum(scaled_contrast))  # Should be ~0
#> [1] -5.551115e-17

Advanced Topics

Custom Contrast Functions

You can generate complex patterns via helper generators that return contrast specifications. For example, a sliding-window set that compares adjacent, disjoint windows across an ordered factor:

# Five ordered levels
lvl <- paste0("L", 1:5)

# Build a small model with an ordered factor
set.seed(111)
emod_sliding <- event_model(
  onset ~ hrf(level, contrasts = sliding_window_contrasts(levels = lvl, facname = "level", window_size = 1)),
  data = data.frame(
    onset = sort(runif(50, 0, 350)),
    level = factor(sample(lvl, 50, replace = TRUE), levels = lvl, ordered = TRUE),
    block = factor(rep(1, 50))
  ),
  block = ~ block,
  sampling_frame = sframe
)

# Inspect the generated contrasts
names(contrast_weights(emod_sliding))
#> [1] "level#win_L1_vs_L2" "level#win_L2_vs_L3" "level#win_L3_vs_L4"
#> [4] "level#win_L4_vs_L5"

For more targeted patterns (e.g., specific basis functions or continuous terms), use column_contrast() with regex to match design-matrix columns.

Troubleshooting Common Issues

1. Rank-Deficient Contrasts

# This will cause issues - contrast is not estimable
bad_data <- data.frame(
  onset = c(10, 30),
  condition = factor(c("A", "A")),  # Only one level!
  block = factor(c(1, 1))
)

# This will warn about the issue
tryCatch({
  bad_model <- event_model(
    onset ~ hrf(condition, contrasts = pair_contrast(~ condition == "A", ~ condition == "B", name = "A_vs_B")),
    data = bad_data,
    block = ~ block,
    sampling_frame = sampling_frame(50, TR = 2)
  )
}, error = function(e) {
  print("Error: Cannot create contrast for single-level factor")
})

2. Multicollinearity in Contrasts

# Check for multicollinearity in your design matrix
cc <- check_collinearity(design_matrix(emodel_with_contrast), threshold = 0.9)
if (!cc$ok) cc$pairs

Summary

Define contrasts inline, extract weights cleanly, and validate before analysis. Use t-contrasts for directional tests and F-contrasts for omnibus effects. - Complex designs: handle factorial, parametric, and multi‑run contrasts - Validation tools: ensure contrasts are properly specified and estimable

Remember to: - Plan contrasts based on your hypotheses - Validate contrast properties before analysis - Document your contrast specifications for reproducibility

For more information on event and baseline models, see: - vignette("a_01_introduction") - vignette("a_03_baseline_model") - vignette("a_04_event_models")