Project Voxel-Level Data using a HATSA Projector — project_voxels.hatsa

Projects voxel-level time series data into the common aligned space defined by a hatsa_projector object. This method uses Nyström extension.

Usage

# S3 method for class 'hatsa_projector'
project_voxels(
  object,
  voxel_timeseries_list,
  voxel_coords,
  parcel_coords,
  n_nearest_parcels = 10,
  kernel_sigma = "auto",
  W_vox_parc = NULL,
  data_type = c("timeseries", "coefficients"),
  ...
)

Arguments

object: A fitted hatsa_projector object.
voxel_timeseries_list: A list of voxel time-series matrices (T_i x V_v). It is assumed that the i-th element of this list corresponds to the i-th subject as stored in the hatsa_projector object (e.g., for U_original_list).
voxel_coords: A numeric matrix (V_v x 3) of voxel coordinates.
parcel_coords: A numeric matrix (V_p x 3) of parcel centroid coordinates corresponding to the parcellation used to fit the object.
n_nearest_parcels: Integer, number of nearest parcels for Nyström extension. Passed to compute_voxel_basis_nystrom.
kernel_sigma: Numeric or "auto", kernel bandwidth for Nyström extension. Passed to compute_voxel_basis_nystrom.
W_vox_parc: Optional pre-computed voxel-to-parcel affinity matrix (V_v x V_p). If provided, this matrix is reused for all subjects and the internal k-NN search and Gaussian kernel calculation are skipped. If NULL (default), the affinity matrix is computed once using n_nearest_parcels and kernel_sigma and then reused.
data_type: Character string, either "timeseries" (default) or "coefficients". If "timeseries", the projection involves scaling by `1/T_i` (number of time points) to estimate covariance with the basis. If "coefficients" (e.g., for projecting beta maps or statistical maps), this scaling is omitted.
...: Additional arguments passed to compute_voxel_basis_nystrom.

Value

A list of aligned voxel coefficient matrices. Each element `[[i]]` is a matrix of dimensions T_i x k, where T_i is the number of time points (or rows) in the input `voxel_timeseries_list[[i]]`, and k is the number of HATSA components. These represent the projection coefficients of the voxel data onto the aligned HATSA basis for each subject.

Details

It is assumed that voxel_coords and parcel_coords are in the same Cartesian coordinate system and units (e.g., millimeters in an MNI-aligned space). The function includes heuristic checks for gross inconsistencies in these coordinate systems (see ..validate_coordinate_inputs for details on checks performed), issuing messages if potential issues like substantially different ranges or scales are detected. The Nyström extension is formulated to be consistent with the unnormalized graph Laplacian used in the core HATSA algorithm for parcel-level decomposition. See compute_voxel_basis_nystrom for more details on parameters like row_normalize_W that control the Nyström calculation.

Examples

# This is a conceptual example. For it to run, you need a fitted hatsa_projector.
# First, let's set up parameters for a minimal run_hatsa_core call.
set.seed(456)
N_subj_fit <- 2
V_parc_fit <- 20 # Number of parcels in the fitted model
k_comp_fit <- 3   # Number of components in the fitted model
T_time_fit <- 40  # Number of time points for parcel data

# Generate mock parcel-level data for fitting HATSA
subject_parcel_data <- lapply(1:N_subj_fit, function(i) {
  matrix(stats::rnorm(T_time_fit * V_parc_fit), ncol = V_parc_fit)
})
anchor_idx_fit <- sample(1:V_parc_fit, min(V_parc_fit, 5))

# Fit a hatsa_projector object (requires Matrix and RSpectra)
fitted_hatsa_model <- NULL
if (requireNamespace("Matrix", quietly = TRUE) && 
    requireNamespace("RSpectra", quietly = TRUE) &&
    exists("run_hatsa_core")) {
  fitted_hatsa_model <- tryCatch({
    run_hatsa_core(
      subject_data_list = subject_parcel_data,
      anchor_indices    = anchor_idx_fit,
      spectral_rank_k = k_comp_fit,
      k_conn_pos        = min(5, V_parc_fit -1),
      k_conn_neg        = min(2, V_parc_fit -1),
      n_refine          = 1
    )
  }, error = function(e) NULL) # Return NULL on error for example
}
#> Proceeding to GPA with 2 subjects having valid anchor matrices (5 parcel + 0 task rows).

if (!is.null(fitted_hatsa_model)) {
  # Now, prepare data for project_voxels
  N_subj_proj <- N_subj_fit # Number of subjects for voxel projection
  V_vox_proj  <- 50         # Number of voxels
  T_time_vox  <- 35         # Number of time points for voxel data

  # Mock voxel time-series data
  voxel_ts_list <- lapply(1:N_subj_proj, function(i) {
    matrix(stats::rnorm(T_time_vox * V_vox_proj), ncol = V_vox_proj)
  })

  # Mock coordinates
  voxel_coords_map <- matrix(stats::rnorm(V_vox_proj * 3), ncol = 3)
  # Parcel coords must match the V_p used when fitting the model
  parcel_coords_map <- matrix(stats::rnorm(V_parc_fit * 3), ncol = 3) 

  # Project voxel data
  projected_vox_coeffs <- project_voxels(
    object = fitted_hatsa_model,
    voxel_timeseries_list = voxel_ts_list,
    voxel_coords = voxel_coords_map,
    parcel_coords = parcel_coords_map,
    n_nearest_parcels = 5, # Nystrom param
    kernel_sigma = "auto"    # Nystrom param
  )

  # print(str(projected_vox_coeffs, max.level=1))
  # if (length(projected_vox_coeffs) > 0) {
  #   print(dim(projected_vox_coeffs[[1]])) # Should be T_time_vox x k_comp_fit
  # }
} else {
  if (interactive()) message("Skipping project_voxels example: fitted_hatsa_model not created.")
}