Spatial operations: masking, alignment, and mapping

Why mask?

In whole-brain analyses, some subjects may have signal dropout or incomplete coverage in certain brain regions. Those missing values show up as NA in the data arrays. If you run a group-level reduction on data that contains NAs, the missing values propagate into the result and you lose information at locations where most subjects do have data.

The mask() verb lets you decide, up front, which spatial locations have enough data to be worth analysing. It removes the rest before any downstream computation.

Setting up the example

Let’s inject some missing values into our test data to see masking in action. We create a fresh GDS with NAs already in place:

# ROI_3 is missing for subject 1; ROI_7 is missing for subject 2
assay(g, "beta")[c(3, 7), , 1]
#>          [,1]       [,2]      [,3]     [,4]
#> [1,]       NA -0.3888607 0.8280826 2.035104
#> [2,] 2.511522         NA 0.7427306 0.215541

Now ROI_3 and ROI_7 each have one subject with missing data.

Data availability across ROIs and subjects. Grey cells indicate missing values.

Intersection mask (the default)

The simplest policy is "intersection": keep only locations where every subject has valid (non-NA) data. This is the strictest rule—if even one subject is missing, that location is dropped.

g_masked <- mask(g) |> compute()
explain(g_masked)
#> GDS
#>   dims: 8 x 4 x 1 [sample x subject x contrast]
#>   space: parcels n=8
#>   subjects: 4 (sub-01, sub-02, sub-03, ...)
#>   contrasts: 1 (faces)
#>   assays: beta[location], var[variance], se[stdev], t[t]

ROI_3 and ROI_7 were dropped because they had at least one NA. The remaining 8 ROIs all have complete data.

Threshold mask

The intersection rule can be wasteful when you have many subjects—one missing subject out of 50 would discard the location entirely. The "threshold" rule is more permissive: it keeps a location if the fraction of non-missing values exceeds a threshold you choose.

# Keep ROIs where at least 75% of values are non-missing
g_thresh <- mask(g, MaskPolicy(rule = "threshold", threshold = 0.75)) |>
  compute()

explain(g_thresh)
#> GDS
#>   dims: 10 x 4 x 1 [sample x subject x contrast]
#>   space: parcels n=10
#>   subjects: 4 (sub-01, sub-02, sub-03, ...)
#>   contrasts: 1 (faces)
#>   assays: beta[location], var[variance], se[stdev], t[t]

With a 75% threshold, ROI_3 and ROI_7 both survive: each has 3 out of 4 subjects present (75%), which meets the cutoff.

Number of ROIs retained under different masking strategies.

Intersection masking is conservative—it discards any ROI with even one missing subject. Threshold masking recovers those ROIs when most subjects have data.

Other mask rules

Here is the full set of options:

A custom mask function receives the list of assay arrays and returns one TRUE/FALSE per sample. For example, to keep only ROIs with a strong average signal:

high_signal <- function(arrays) {
  rowMeans(arrays$beta, na.rm = TRUE) > 0.5
}

g_custom <- mask(g, MaskPolicy(rule = "custom", custom = high_signal)) |>
  compute()

dim(assay(g_custom, "beta"))
#> [1] 8 4 1

Masking in a pipeline

Because mask() is lazy, you can slot it into a larger pipeline. Here we mask first, then run a group analysis with FDR correction—all as a single plan:

result <- as_plan(g) |>
  mask(MaskPolicy(rule = "intersection")) |>
  reduce(method = "fixed") |>
  posthoc("fdr:bh") |>
  compute()

names(assays(result))
#>  [1] "beta"   "var"    "beta_g" "var_g"  "se_g"   "z_g"    "p_g"    "Q"     
#>  [9] "I2"     "n_eff"  "z"      "p"      "q"

The mask runs before the reduction, so the group model only sees complete cases.

Why align?

When you use dimensionality reduction (PCA, ICA) or functional alignment (hyperalignment) on fMRI data, each subject ends up in their own coordinate system. Component 1 for subject A is not the same as component 1 for subject B—the axes are rotated differently.

Before you can average across subjects, you need to apply a per-subject rotation matrix that maps each subject’s data into a shared “group” space. This is what align() does.

Creating test data in a latent space

Let’s create a GDS where each subject has 5 latent components, 2 experimental conditions, and a subject-specific coordinate system:

k <- 5       # number of latent components
n_sub <- 3   # number of subjects
n_con <- 2   # number of contrasts

betas_lat <- array(rnorm(k * n_sub * n_con), dim = c(k, n_sub, n_con))
vars_lat  <- array(runif(k * n_sub * n_con, 0.01, 0.1), dim = c(k, n_sub, n_con))

# Define the native and target spaces
native_space <- space_basis(k, basis_name = "native_pca")
group_space  <- space_basis(k, basis_name = "group_pca")

subs <- paste0("sub-", sprintf("%02d", 1:n_sub))

g_latent <- as_gds(
  list(beta = betas_lat, var = vars_lat),
  space     = native_space,
  subjects  = subs,
  contrasts = paste0("cond_", 1:n_con)
)

The space_basis() constructor says “these samples are latent components.” We have a native space (each subject’s own PCA axes) and a target group space (shared axes).

Providing per-subject rotation matrices

In a real analysis, rotation matrices come from Procrustes alignment or hyperalignment. Here we simulate them using random orthogonal matrices:

random_orthogonal <- function(k) {
  qr.Q(qr(matrix(rnorm(k * k), k, k)))
}

rotations <- setNames(
  lapply(1:n_sub, function(i) random_orthogonal(k)),
  subs
)

# Each rotation is a k-by-k orthogonal matrix
dim(rotations[[1]])
#> [1] 5 5

Each matrix is k x k (5 x 5 here). Multiplying a subject’s data by their rotation matrix maps them from native coordinates into group coordinates.

Bundling rotations into a MapFamily

fmrigds needs to know not just the matrices but also what spaces they connect and what properties they have. You bundle this information into a MapFamily object. For orthogonal rotations, use the OrthogonalFamily() helper:

alignment <- OrthogonalFamily(
  name                = "procrustes",
  from_space          = native_space,
  to_space            = group_space,
  matrices_by_subject = rotations
)

The orthogonal = TRUE trait tells fmrigds that these matrices satisfy R'R = I, which simplifies variance propagation—the variance of a rotated vector is just the rotated variance, with no cross-terms.

Running the alignment

Now apply the alignment and compute:

g_aligned <- align(g_latent, alignment) |> compute()

explain(g_aligned)
#> GDS
#>   dims: 5 x 3 x 2 [sample x subject x contrast]
#>   space: basis
#>   subjects: 3 (sub-01, sub-02, sub-03)
#>   contrasts: 2 (cond_1, cond_2)
#>   assays: beta[location], var[variance], se[stdev], t[t]
#>   map families: procrustes

Betas before (native space) and after (group space) alignment for subject 1, condition 1.

After alignment, the space descriptor has changed from native_pca to group_pca. Every subject’s data has been rotated through their respective matrix, and variances have been propagated through the rotation automatically. The component loadings change because the coordinate axes have been remapped.

Alignment followed by group analysis

The common workflow is to align first, then reduce across subjects. Because subjects are now in a common space, the group average is meaningful:

result <- as_plan(g_latent) |>
  align(alignment) |>
  reduce(method = "fixed") |>
  compute()

explain(result)
#> GDS
#>   dims: 5 x 1 x 2 [sample x subject x contrast]
#>   space: basis
#>   subjects: 1 (meta)
#>   contrasts: 2 (cond_1, cond_2)
#>   assays: beta[location], var[variance], beta_g[?], var_g[?], se_g[?], z_g[?], p_g[?], Q[?], I2[?], n_eff[n_eff], z[z], p[p]
#>   map families: procrustes

Other alignment types

fmrigds provides helpers for different kinds of per-subject maps:

All follow the same pattern: bundle per-subject operators into a family object, then pass it to align().

Why transform spaces?

Sometimes you want to change the spatial resolution of your data. Common scenarios:

• Voxels to parcels: aggregate thousands of voxels into a smaller set of brain regions.
• Fine parcellation to coarse: collapse a detailed atlas (e.g., 400 regions) into broader networks (e.g., 7 networks).
• Projection: map data from one coordinate system to another using a known linear transformation.

map_to() applies a matrix transformation to the sample axis. If the data has effect estimates and variances, the variances are propagated through the transformation automatically.

Setting up a mapping

The transformation matrix has dimensions target_samples x source_samples. Each row defines how one target location is computed from the source locations. For averaging, each row sums to 1:

# 10 ROIs that we want to average into 3 regions:
#   Region "frontal"   = average of ROIs 1-3
#   Region "parietal"  = average of ROIs 4-6
#   Region "occipital" = average of ROIs 7-10
g_clean <- as_gds(
  list(
    beta = array(rnorm(40, mean = 1.5), dim = c(10, 4, 1)),
    var  = array(runif(40, 0.01, 0.05), dim = c(10, 4, 1))
  ),
  space     = space_parcels(paste0("ROI_", 1:10)),
  subjects  = paste0("sub-", sprintf("%02d", 1:4)),
  contrasts = "faces"
)

# Build the averaging matrix (3 target regions x 10 source ROIs)
W <- matrix(0, nrow = 3, ncol = 10)
W[1, 1:3]  <- 1/3   # frontal = mean of ROI 1-3
W[2, 4:6]  <- 1/3   # parietal = mean of ROI 4-6
W[3, 7:10] <- 1/4   # occipital = mean of ROI 7-10

target <- space_parcels(c("frontal", "parietal", "occipital"))

The mapping matrix aggregates 10 fine-grained ROIs into 3 broader brain regions.

Each coloured cell indicates a non-zero weight. ROIs 1–3 map to “frontal”, 4–6 to “parietal”, and 7–10 to “occipital”.

Applying the transformation

Wrap the matrix in a map_linear() object that records the source and target spaces, then pass it to map_to():

mapping <- map_linear(
  from_space = space_parcels(paste0("ROI_", 1:10)),
  to_space   = target,
  operator   = W
)

g_regions <- map_to(g_clean, target_space = target, map = mapping) |>
  compute()

explain(g_regions)
#> GDS
#>   dims: 3 x 4 x 1 [sample x subject x contrast]
#>   space: parcels n=3
#>   subjects: 4 (sub-01, sub-02, sub-03, ...)
#>   contrasts: 1 (faces)
#>   assays: beta[location], var[variance], se[stdev], t[t]
assay(g_regions, "beta")
#> , , 1
#> 
#>          [,1]     [,2]      [,3]     [,4]
#> [1,] 1.478718 1.537494 1.3347244 1.814508
#> [2,] 1.059719 1.110867 1.3439407 1.649642
#> [3,] 1.206013 2.233578 0.8522503 1.461017

Effect estimates before (10 ROIs) and after (3 regions) spatial aggregation, for subject 1.

Each target region now holds the average of its source ROIs. The variances have also been propagated: under the default assumption that source variances are independent, Var(Wx) = W diag(v) W'.

How uncertainty propagation works

When you transform data through a linear map, the variances need to be transformed too. By default, map_to() assumes the source locations have independent variances (no spatial correlation). If your data has spatial correlations, you can provide a covariance function:

# Default: independent variances
map_to(g, target, mapping, uncertainty = UncertaintyRule("independent"))

# Correlated variances: provide a function that returns the covariance matrix
map_to(g, target, mapping, uncertainty = UncertaintyRule(
  "cov_provider",
  cov_provider = my_cov_function
))

Combining test statistics without effect estimates

Sometimes your GDS contains only test statistics (z-scores or p-values) without the underlying effect estimates and variances. In that case, linear averaging does not apply. Instead, map_to() can combine the statistics using established meta-analytic methods:

• Stouffer’s method: produces a weighted z-score for each target location.
• Fisher’s method: combines p-values using a chi-squared statistic.

# A GDS with only z-scores (no beta or var)
z_vals <- array(rnorm(40, mean = 2), dim = c(10, 4, 1))
g_zonly <- as_gds(
  list(z = z_vals),
  space     = space_parcels(paste0("ROI_", 1:10)),
  subjects  = paste0("sub-", sprintf("%02d", 1:4)),
  contrasts = "faces"
)

g_combined <- map_to(
  g_zonly,
  target_space = target,
  map = W,
  combine = "stouffer"
) |> compute()

names(assays(g_combined))
#> [1] "z" "p"
assay(g_combined, "z")
#> , , 1
#> 
#>          [,1]     [,2]     [,3]     [,4]
#> [1,] 4.345908 2.167603 4.749368 3.912684
#> [2,] 4.571568 3.726133 2.978678 4.240634
#> [3,] 4.361189 2.953622 2.991708 4.906879

The combined z-scores reflect the pooled evidence across the source locations in each target region.