---
title: Contrasts in fmridesign
description: Guide to defining and using contrasts for hypothesis testing in fMRI
  designs.
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
params:
  family: lapis
  preset: adobe
  base_size: 13
  content_width: 80
vignette: |
  %\VignetteIndexEntry{Contrasts in fmridesign}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
css: albers.css
resource_files:
- albers.css
- albers.js
includes:
  in_header: |-
    <script src="albers.js"></script>
    <script>document.addEventListener('DOMContentLoaded',function(){document.body.classList.add('palette-lapis','preset-adobe');});</script>

---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE, comment = "#>", fig.align = "center", fig.retina = 2,
  out.width = "100%", fig.width = 7, fig.asp = 0.618, message = FALSE, warning = FALSE
)
set.seed(123)
old_opts <- options(pillar.sigfig = 7, width = 80)
old_theme <- ggplot2::theme_get()
suppressPackageStartupMessages({
  library(fmridesign)
  library(ggplot2)
})
ggplot2::theme_set(ggplot2::theme_minimal(base_size = params$base_size))
```

<!-- Removed per-page palette injection; default palette is red via CSS root variables. -->

```{r content-width, echo=FALSE, results='asis'}
cat(sprintf('<style>:root{--content:%sch}</style>', params$content_width))
```

## 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:

```{r basic_contrast}
# 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)
```

### Quick Validation

Validate all attached contrasts for a model in one call:

```{r validate_attached_contrasts}
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
```

## Advanced Contrast Specifications

### Pairwise Contrasts

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

```{r pairwise_contrasts}
# 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)
```

### 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.

```{r polynomial_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)
pc <- poly_contrasts[[1]]
round(head(pc$weights, 6), 3)

# If you want a specific component, select a column (e.g., linear = column 1)
linear_weights <- pc$weights[, 1, drop = FALSE]
head(linear_weights)
```

## Factorial Design Contrasts

### Main Effects and Interactions

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

```{r factorial_contrasts}
# 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))
```

## 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.

```{r parametric_contrasts}
# 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"]]
colnames(dm)[idx[["RT"]]]
```

## F-contrasts for Omnibus Tests

F-contrasts test multiple parameters simultaneously:

```{r f_contrasts}
# 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)
```

## Working with Contrast Weights

### Extracting and Manipulating Weights

```{r manipulate_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)
```

### 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()`:

```{r validate_contrasts}
# Validate the custom contrast against the design implied by the model
validate_contrasts(simple_model, weights = custom_contrast)
```

## Contrasts for Multi-Run Designs

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

```{r multirun_contrasts}
# 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)
```

## Contrast Specification Best Practices

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

### 2. Scaling and Normalization

```{r contrast_scaling}
# 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)
print(sum(scaled_contrast))  # Should be ~0
```

<!-- Section on multiple-comparison correction removed: not standard within-model in typical fMRI GLM workflows. Downstream group analysis packages handle control of false positives across space/contrasts. -->

<!-- Section on downstream statistical fitting removed to keep the vignette focused on
     building and validating contrasts. GLM fitting and inference vary by toolbox and
     are typically handled outside fmridesign. -->
## 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:

```{r custom_contrast_function}
# 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))
```

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

```{r rank_deficient}
# 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

```{r multicollinearity}
# 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")`

```{r cleanup, include=FALSE}
options(old_opts)
ggplot2::theme_set(old_theme)
```