1 Overview

NeuroVault, is an open repository for sharing statistical maps, parcellations, and atlases of the human brain. Scientists can upload the brain images that they use for their studies to the online repository, to allow for interactive visualization of images while also contributing data entries for meta-analysis.

The images in Neurovault are highly heterogeneous and susceptible to human reporting errors in manually enter image metadata. This document automatically curates the Neurovault repository to produce a final dataset that is amenable to voxel-level meta-analysis. Here, we develop and run a preprocessing pipeline to ensure the Neurovault repository entries correspond to group-level statistical images and that images are registered to a common template (MNI 2mm). Details of preprocessing are described below. As a final step, we convert the images to effect sizes (the robust effect size index; RESI) and perform statistical analysis to remove outliers.

NeuroVault is living growing database of images. This preprocessing pipeline is automatically run on the first of each month to build a comprehensive preprocessed dataset that is freely available for download. The curated dataset is available on Zenodo. In our preprint, we demonstrate the utility of the curated dataset through unsupervised and supervised meta-analyses. This document and the dataset were last updated on 2026-09-30.

2 Download the Neurovault data

2.1 Metadata download and exclusionary criteria

The Neurovault metadata include manually entered information at the study- and image-level describing the type of data that were uploaded, including:

analysis_level: indicates whether the image is a group-level or individual-level statistical map.
not_mni: indicates whether the image is in MNI space.
is_thresholded: indicates whether the image is thresholded or unthresholded.
map_type: indicates the statistical type of the image (e.g., T, Z, F).
modality: indicates the imaging modality (e.g., structural, functional).
perc_bad_voxels: Proportion of voxels with 0 intensity or NaN values.
perc_voxels_outside: Proportion of voxels outside the brain.
number_of_subjects: Sample size of the study
DOI: The DOI for the publication associated with the repository entry.
paper_url: The URL for the publication associated with the repository entry.

We download the Neurovault study-level and image-level metadata, then apply the following quality exclusionary criteria to the metadata to ensure uniformity of the group-level statistical images.

Retain only group-level images, which represent group-level statistical maps rather than individual-level contrasts (analysis_level == "group").
Exclude images that are not in MNI space (not_mni == FALSE).
Keep only unthresholded maps (is_thresholded == FALSE), which ensures all statistical values computed in the image are included. This avoids biased caused by only reporting significant results
Restrict to statistical image types T and Z images, which are most commonly used in group-level inference. This step removes less frequently reported F images, and non-statistical images.
Exclude those that have sample size of 0 or NaN.

See Table 1 and Figure 1 for a summary of the intersection of these exclusionary criteria.

## This chunk would take some time, as it creates the file neurovaultMetadata, containing all the study-level images

# Check if data was already downloaded this month
rds_file_path <- get_monthly_rds_path("neurovaultMetadata")
file_exists_this_month <- file.exists(rds_file_path)

if (!file_exists_this_month) {
  # Downloads metadata, saves it internally, and returns it so `nv` is populated
  download_and_save_metadata_fn <- function(rds_dir) {
    library(httr)
    library(jsonlite)
    
    # Helper function to get date suffix
    get_date_suffix <- function() {
      format(Sys.Date(), "%Y-%m")
    }
    
    # Helper function to get monthly RDS path
    get_monthly_rds_path <- function(base_name, dir = rds_dir) {
      date_suffix <- get_date_suffix()
      current_file <- file.path(dir, paste0(base_name, "_", date_suffix, ".rds"))
      
      pattern <- paste0("^", base_name, "_[0-9]{4}-[0-9]{2}\\.rds$")
      all_files <- list.files(dir, pattern = pattern, full.names = TRUE)
      
      old_files <- all_files[all_files != current_file]
      if (length(old_files) > 0) {
        invisible(lapply(old_files, file.remove))
      }
      
      return(current_file)
    }
    
    Sys.setenv(CURL_CA_BUNDLE = "/etc/ssl/certs/ca-certificates.crt")
    httr::GET("https://api.github.com")
    
    url <- 'https://neurovault.org/api/collections/?format=json'
    
    nv = list()
    i = 0
    while (!is.null(url)) {
      i = i + 1
      cat(sprintf("[%s] Downloading collections batch %d\n", Sys.time(), i))
      # httr::config, not plotly::config (plotly is loaded earlier and masks it)
      response <- GET(url, httr::config(cainfo="/etc/ssl/certs/ca-certificates.crt"))
      data <- content(response, "text")
      json_data <- fromJSON(data, flatten = TRUE)
      nv[[i]] = json_data$results
      
      url <- json_data[['next']]
    }
    nv = do.call(rbind, nv)
    # subset to studies that have more than 0 images
    nv = nv[nv$number_of_images>0,]
    
    # Save directly to versioned RDS file
    if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
    rds_path <- get_monthly_rds_path("neurovaultMetadata")
    saveRDS(nv, file = rds_path)
    
    cat(sprintf("[%s] Neurovault metadata saved to %s\n", Sys.time(), rds_path))
    return(nv)
  }
  
  logfiles = c('logfiles/neurovault_collection_download.log', 'logfiles/neurovault_collection_download.err')
  unlink(logfiles)
  nv <- run_with_logging(download_and_save_metadata_fn(rds_dir), logfiles[1], logfiles[2])
  
  message(sprintf("[%s] Download complete. Output saved to %s", 
          Sys.time(), rds_file_path))
} else {
  # Load already downloaded metadata
  nv <- load_monthly_rds("neurovaultMetadata")
  message("Loaded previously downloaded metadata from ", rds_file_path)
}

  nv <- load_monthly_rds("neurovaultMetadata")
nv$DOIexists = factor(ifelse(is.na(nv$DOI), "No", "Yes"), levels=c('Yes', 'No'))
nv$paperURLexists = factor(ifelse(is.na(nv$paper_url), "No", "Yes"), levels=c('Yes', 'No'))
label(nv$paperURLexists) = "Paper URL exists"
label(nv$number_of_images) = "Number of images"
label(nv$DOIexists) = "DOI exists"
table1(~number_of_images + DOIexists + paperURLexists, data = nv, caption = "Study-level metadata.")

## The following two chunks would take some time, as it extracts and downloads all the image-level data and creates the first image tabular dataset that we work with. More than 2 hours

## May 2026: 10743 collections

api_cores = 4
# Check if image list was already downloaded this month
rds_file_path <- get_monthly_rds_path("imageList")
file_exists_this_month <- file.exists(rds_file_path)

if (!file_exists_this_month) {
  # Downloads image list, saves it internally, and returns it so `image_list` is populated
  download_and_save_image_list_fn <- function(collection_ids, n_cores, rds_dir) {
    library(httr)
    library(jsonlite)
    library(parallel)
    
    # Ensure collection_ids is a proper character vector (fix for list coercion error)
    if (is.list(collection_ids)) {
      collection_ids <- unlist(collection_ids)
    }
    collection_ids <- as.character(collection_ids)
    
    cat(sprintf("[%s] [BG] Starting download with %d collection IDs\n", Sys.time(), length(collection_ids)))
    
    # Helper function to get monthly RDS path
    get_monthly_rds_path <- function(base_name, dir = rds_dir) {
      get_date_suffix <- function() {
        format(Sys.Date(), "%Y-%m")
      }
      date_suffix <- get_date_suffix()
      current_file <- file.path(dir, paste0(base_name, "_", date_suffix, ".rds"))
      
      pattern <- paste0("^", base_name, "_[0-9]{4}-[0-9]{2}\\.rds$")
      all_files <- list.files(dir, pattern = pattern, full.names = TRUE)
      
      old_files <- all_files[all_files != current_file]
      if (length(old_files) > 0) {
        invisible(lapply(old_files, file.remove))
      }
      
      return(current_file)
    }
    
    # Download function with pagination support and error handling
    download_collection_safe <- function(collection_id) {
      tryCatch({
        collection_id <- as.character(collection_id)
        cat(sprintf("[%s] [BG] Downloading collection %s\n", Sys.time(), collection_id))
        
        url <- paste0("https://neurovault.org/api/collections/", collection_id, "/images/?format=json")
        all_results <- list()
        page <- 0
        
        # Paginate through all results
        while (!is.null(url)) {
          page <- page + 1
          response <- GET(url, timeout(30))
          
          if (status_code(response) != 200) {
            cat(sprintf("[%s] [BG] Collection %s page %s: HTTP %s\n", Sys.time(), collection_id, page, status_code(response)))
            break
          }
          
          data <- content(response, "text")
          json_data <- fromJSON(data, flatten = TRUE)
          
          # Extract results from this page
          if (!is.null(json_data$results) && is.data.frame(json_data$results) && nrow(json_data$results) > 0) {
            all_results[[page]] <- json_data$results
          }
          
          # Check for next page
          url <- json_data[['next']]
        }
        
        # Combine all pages
        if (length(all_results) == 0) {
          cat(sprintf("[%s] [BG] Collection %s: No results\n", Sys.time(), collection_id))
          return(data.frame())
        }
        
        result <- do.call(rbind, all_results)
        
        # Validate result is a dataframe
        if (!is.data.frame(result)) {
          cat(sprintf("[%s] [BG] Collection %s: invalid result type (%s)\n", Sys.time(), collection_id, class(result)[1]))
          return(data.frame())
        }
        
        cat(sprintf("[%s] [BG] Collection %s: OK (%s pages, %s total rows, %s cols)\n", Sys.time(), collection_id, page, nrow(result), ncol(result)))
        return(result)
        
      }, error = function(e) {
        cat(sprintf("[%s] [BG] Collection %s: ERROR - %s\n", Sys.time(), collection_id, e$message))
        return(data.frame())
      })
    }
    
    # Download all collections
    cat(sprintf("Starting mclapply", Sys.time(), n_cores))
    image_list <- mclapply(collection_ids, download_collection_safe, mc.cores = n_cores)
    
    # Validate we got a list
    if (!is.list(image_list)) {
      stop(sprintf("ERROR: mclapply did not return a list, got %s", class(image_list)))
    }
    cat(sprintf("mclapply returned results\n", Sys.time(), length(image_list)))
    
    # Keep only data frames with >0 rows
    cat(sprintf("Filtering dataframes with >0 rows\n", Sys.time()))
    is_valid <- sapply(image_list, function(x) is.data.frame(x) && nrow(x) > 0)
    valid_count <- sum(is_valid)
    
    if (valid_count == 0) {
      stop("ERROR: No valid dataframes downloaded. Check collection_ids and API.")
    }
    
    image_list <- image_list[is_valid]
    
    # Save list directly to versioned RDS file
    if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
    rds_path <- get_monthly_rds_path("imageList")
    
    cat(sprintf("Saving image list (%d collections) to %s\n", length(image_list), rds_path))
    saveRDS(image_list, file = rds_path, compress = "gzip")
    
    cat("Image list saved successfully\n")
    return(image_list)
  }
  
  
  if (is.null(nv) || nrow(nv) == 0) {
    stop("ERROR: neurovaultMetadata is empty. Run getStudyLevelData chunk first.")
  }
  
  # Extract and validate collection IDs
  collection_ids <- nv$id
  if (is.list(collection_ids)) {
    collection_ids <- unlist(collection_ids)
  }
  collection_ids <- as.character(collection_ids)
  
  cat(sprintf("[%s] Extracted collection IDs\n", Sys.time(), length(collection_ids)))
  
  # Create logfiles directory if needed
  logdir <- "logfiles"
  if (!dir.exists(logdir)) dir.create(logdir, showWarnings = FALSE)
  logfiles = c('logfiles/neurovault_download.log', 'logfiles/neurovault_download.err')
  unlink(logfiles)
  
  image_list <- run_with_logging(
    download_and_save_image_list_fn(collection_ids, api_cores, rds_dir),
    logfiles[1], logfiles[2]
  )
  
  message(sprintf("[%s] Downloaded %d collections. Output: %s", 
            Sys.time(), length(collection_ids), logfiles[1]))
} else {
  # Load already downloaded list
  image_list <- load_monthly_rds("imageList")
  message("Loaded previously downloaded image list from ", rds_file_path)
}

# Check if merged image-level data already exists this month
rds_file_path <- get_monthly_rds_path("fullimageMetadata")
file_exists_this_month <- file.exists(rds_file_path)

#res = unlist(by(images$collection_id, images$collection_id, length)); res = data.frame(collection_id = names(res), numberinImages=res) 
#test = merge(res, nv[, c('id', 'number_of_images')], by.x='collection_id', by.y='id')

if (!file_exists_this_month) {
  
  if (is.null(image_list) || length(image_list) == 0) {
    stop("ERROR: imageList is empty. Run getImageLevelData_section1.1_download chunk first and wait for it to complete.")
  }
  
  cat(sprintf("[%s] Loaded image list with dataframes\n", Sys.time(), length(image_list)))
  
  # Find common column names across all dataframes (columns in ALL collections)
  cat(sprintf("[%s] Analyzing column overlap across collections\n", Sys.time(), length(image_list)))
  all_colnames <- lapply(image_list, colnames)
  
  # Columns that appear in ALL collections
  # SNV: This excludes important columns
  # common_cols <- Reduce(intersect, all_colnames)
  # cat(sprintf("[%s] Found %d columns in ALL collections\n", Sys.time(), length(common_cols)))
  
  # Find columns that appear in at least 80% of collections
  col_freq <- table(unlist(all_colnames))
  common_cols <- names(col_freq)[col_freq >= 0.8 * length(image_list)]
  cat(sprintf("[%s] Fallback: columns in ≥80%% of collections\n", Sys.time(), length(common_cols)))
  
  # Subset each dataframe to common columns only
 cat(sprintf("[%s] Subsetting to %d common columns\n", Sys.time(), length(common_cols)))
  image_list <- lapply(image_list, function(df) {
    # Get columns that exist in this dataframe
    cols_exist <- intersect(common_cols, colnames(df))
    
    # Get missing columns
    cols_missing <- setdiff(common_cols, colnames(df))
    
    # Start with existing columns
    result <- df[, cols_exist, drop = FALSE]
    
    # Add missing columns filled with NA
    if (length(cols_missing) > 0) {
      result[, cols_missing] <- NA
    }
    
    # Reorder to match common_cols order
    result[, common_cols, drop = FALSE]
  })
  
  # Combine dataframes using rbind (much faster than merge)
  cat(sprintf("[%s] Combining %d dataframes via rbind...\n", Sys.time(), length(image_list)))
  start_time <- Sys.time()
  images <- do.call(rbind, image_list)
  elapsed <- difftime(Sys.time(), start_time, units = "secs")
  cat(sprintf("[%s] rbind complete in %.1f sec: %d rows × %d columns\n", 
              Sys.time(), elapsed, nrow(images), ncol(images)))
  
  # Sort all columns by missingness
  cat(sprintf("[%s] Computing column completeness for %d columns\n", Sys.time(), ncol(images)))
  col_miss <- colMeans(is.na(images))
  col_miss = sort(col_miss)
  col_keep <- names(col_miss)[col_miss < 0.8]
  images <- images[, col_keep]
  
  
  # Save to versioned RDS file
  if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
  rds_file_path <- get_monthly_rds_path("fullimageMetadata")
  
  cat(sprintf("[%s] Saving merged data to %s\n", Sys.time(), rds_file_path))
  saveRDS(images, file = rds_file_path, compress = "gzip")
  
  cat(sprintf("[%s] ✓ Full image metadata saved successfully\n", Sys.time()))
} else {
  # Load already merged metadata
  images <- load_monthly_rds("fullimageMetadata")
  message("Loaded previously merged image metadata from ", rds_file_path)
}

images <- images %>% mutate(
  image_type      = factor(image_type,      levels = names(sort(table(image_type),      decreasing = TRUE))),
  map_type        = factor(map_type,        levels = names(sort(table(map_type),        decreasing = TRUE))),
  analysis_level  = factor(analysis_level,  levels = names(sort(table(analysis_level),  decreasing = TRUE))),
  is_thresholded  = factor(is_thresholded,  levels = names(sort(table(is_thresholded),  decreasing = TRUE))),
  not_mni         = factor(not_mni,         levels = names(sort(table(not_mni),         decreasing = TRUE)))
)

table1Form = ~ number_of_subjects + modality + image_type + map_type + analysis_level + is_thresholded + not_mni + perc_bad_voxels + perc_voxels_outside

table1(table1Form,
       data = images,
       caption = "Table 1: Neurovault image-level metadata summary.")

Table 1: Neurovault image-level metadata summary.
	Overall (N=624521)
number_of_subjects
Mean (SD)	194 (61600)
Median [Min, Max]	1.00 [-275, 32200000]
Missing	350692 (56.2%)
modality
	33055 (5.3%)
Diffusion MRI	1056 (0.2%)
EEG	49 (0.0%)
fMRI-BOLD	512049 (82.0%)
fMRI-CBF	276 (0.0%)
fMRI-CBV	26 (0.0%)
MEG	86 (0.0%)
Other	71843 (11.5%)
PET FDG	143 (0.0%)
PET other	280 (0.0%)
Structural MRI	4459 (0.7%)
Missing	1199 (0.2%)
image_type
statistic_map	623322 (99.8%)
NIDM results statistic map	1181 (0.2%)
atlas	18 (0.0%)
map_type
univariate-beta map	148231 (23.7%)
variance	138086 (22.1%)
Z map	135497 (21.7%)
T map	83459 (13.4%)
other	61525 (9.9%)
ROI/mask	46084 (7.4%)
P map (given null hypothesis)	7674 (1.2%)
parcellation	1260 (0.2%)
F map	1072 (0.2%)
1-P map ("inverted" probability)	555 (0.1%)
multivariate-beta map	507 (0.1%)
anatomical	497 (0.1%)
Chi squared map	56 (0.0%)
Missing	18 (0.0%)
analysis_level
single-subject	433734 (69.5%)
group	61835 (9.9%)
meta-analysis	55678 (8.9%)
other	2333 (0.4%)
	27 (0.0%)
Missing	70914 (11.4%)
is_thresholded
FALSE	521245 (83.5%)
TRUE	95822 (15.3%)
Missing	7454 (1.2%)
not_mni
FALSE	621704 (99.5%)
TRUE	2799 (0.4%)
Missing	18 (0.0%)
perc_bad_voxels
Mean (SD)	76.9 (15.1)
Median [Min, Max]	77.5 [0, 100]
Missing	18 (0.0%)
perc_voxels_outside
Mean (SD)	20.8 (32.1)
Median [Min, Max]	8.52 [0, 100]
Missing	974 (0.2%)

Table 1 summarizes the full NeuroVault repository prior to any exclusions (N = 624521). Several patterns are worth noting. The repository is dominated by single-subject images and meta-analytic maps; group-level images — the target of this pipeline — comprise only 11.2% of the total. The most common map types are non-statistical maps such as univariate-beta map and variance; T and Z maps together account for a minority of all images. A large proportion of images have no recorded sample size (56.2% missing), reflecting the manually curated nature of the repository. The high median perc_bad_voxels reflects expected zero-valued voxels outside the brain mask rather than data corruption.

# SNV: Warning mapping missing to FALSE in upset plot is suppressed.
library(ComplexUpset)
## UpSet table
images$binarized_group_level <- ifelse(images$analysis_level == "group", TRUE, FALSE)
images$binarized_map_type <- ifelse(images$map_type == "T map" | images$map_type == "Z map", TRUE, FALSE)
images$unthresholded <- ifelse(images$is_thresholded == FALSE, TRUE, FALSE)
images$mni <- ifelse(images$not_mni == FALSE, TRUE, FALSE)
images$nonzero_sample <- ifelse(is.na(images$number_of_subjects) | images$number_of_subjects<=0, FALSE, TRUE)

images_clean <- images %>%
  transmute(
    "Group-level analysis" = as.logical(binarized_group_level),
    "T/Z map"  = as.logical(binarized_map_type),
    "Unthresholded" = as.logical(unthresholded),
    "MNI space" = as.logical(mni),
    "Nonzero sample" = as.logical(nonzero_sample)
  )

upset(
  images_clean,
  intersect = c(
    "Group-level analysis",
    "T/Z map",
    "Unthresholded",
    "MNI space",
    'Nonzero sample'
  ),
  name='Group',
  min_size = 50,
  base_annotations = list(
    "Intersection size" = intersection_size(
      text = list(
        angle = 90,
        vjust = 0.5,
        hjust = -0.02
      )
    )
  )
)

Figure 1: Intersection of NeuroVault image-level exclusionary criteria.

# SNV: clean up larger unused objects
ggsave("poster_upset.png", width=10, height=8, dpi=150)
rm(images_clean, image_list)

Figure 1 shows the an UpSet plot showing the intersection of the five exclusionary criteria applied to the raw NeuroVault image metadata. Each column in the upper bar chart represents a unique combination of criteria, with the bar height indicating how many images fall into that combination. The dot matrix below each bar shows which criteria are active — filled dots connected by a vertical line indicate criteria that are simultaneously satisfied. The horizontal bars on the left show the total number of images satisfying each criterion individually (set size), regardless of the others. The final retained images satisfy all five criteria simultaneously. See Table 2 below for summary statistics after exclusion.

Note:

The analysis_level variable has other types of data, and group-level images make up less than 10% of data.
Including only T and Z maps removes approximately 70% of the data.
0 images have a reported sample size greater than 100,000 (maximum: 32,000). Our approach here is by spot checking the entry with 32,222,222, and it turns out to be a sample size of 32. We have manually fixed it in our dataset. However, for everything else and future additional studies to our pipeline, a study with more than 0 images can reflect metadata entry errors. Thus, implausibly large values will be flagged and removed here (except for the one that we spot checked.

# SNV: I deleted reloading the dataset, here
images_before <- nrow(images)
# subset to group-level analysis
images = images[ which(images$analysis_level=='group'), ]
# subset to unthresholded images
images = images[which(images$is_thresholded==FALSE),]
# Including only T and Z map types
images = images[images$map_type %in% c("T map", "Z map"), ]
# Exclude all not mni
images = images[which(images$not_mni == FALSE), ]
# Exclude sample size 0
images <- images[which(!is.na(images$number_of_subjects) & images$number_of_subjects>0),]

table1(table1Form, data = images, caption = "Table 2: Summary statistics after metadata exclusions.")

Table 2: Summary statistics after metadata exclusions.
	Overall (N=33199)
number_of_subjects
Mean (SD)	66.4 (201)
Median [Min, Max]	28.0 [1.00, 30000]
modality
Diffusion MRI	65 (0.2%)
fMRI-BOLD	17803 (53.6%)
fMRI-CBF	14 (0.0%)
fMRI-CBV	2 (0.0%)
MEG	39 (0.1%)
Other	14672 (44.2%)
PET FDG	45 (0.1%)
PET other	43 (0.1%)
Structural MRI	516 (1.6%)
image_type
statistic_map	33199 (100%)
NIDM results statistic map	0 (0%)
atlas	0 (0%)
map_type
univariate-beta map	0 (0%)
variance	0 (0%)
Z map	2597 (7.8%)
T map	30602 (92.2%)
other	0 (0%)
ROI/mask	0 (0%)
P map (given null hypothesis)	0 (0%)
parcellation	0 (0%)
F map	0 (0%)
1-P map ("inverted" probability)	0 (0%)
multivariate-beta map	0 (0%)
anatomical	0 (0%)
Chi squared map	0 (0%)
analysis_level
single-subject	0 (0%)
group	33199 (100%)
meta-analysis	0 (0%)
other	0 (0%)
	0 (0%)
is_thresholded
FALSE	33199 (100%)
TRUE	0 (0%)
not_mni
FALSE	33199 (100%)
TRUE	0 (0%)
perc_bad_voxels
Mean (SD)	75.2 (11.0)
Median [Min, Max]	78.9 [0, 86.1]
perc_voxels_outside
Mean (SD)	10.3 (19.6)
Median [Min, Max]	0.820 [0, 100]

After applying these exclusions, the dataset is reduced from 624,521 to 33,199 images (Table 2). The retained images are exclusively group-level, unthresholded T and Z maps in MNI space with a valid nonzero sample size. As expected, variables used as exclusionary criteria — analysis_level, is_thresholded, not_mni, and map_type — show no remaining variation, confirming the filters were applied correctly. The dataset is dominated by fMRI-BOLD images (53.6%), with a median sample size of 28 subjects. The mean sample size (66) is notably larger than the median, suggesting a small number of studies with disproportionately large sample sizes that will be addressed in the next section.

2.2 Extracting titles and abstract from Study-Level data

The collection-level title (name) and abstract (description) fields live in the study-level data (nv), not the image-level data (images, whose own name/description columns are per-image map labels) – so we join them in here by collection_id. These are NeuroVault’s own title/description for the collection (often terse or blank), which we supplement below with the associated paper’s actual title/abstract fetched from PubMed by DOI.

2.3 Excluding meta-analytic studies

NeuroVault also hosts collections that are themselves meta-analyses (e.g. coordinate-based/ALE meta-analyses) rather than primary studies, which we do not want in a dataset intended for voxel-level meta-analysis. We flag these with a simple, case-insensitive search for “meta-analysis” and a case-sensitive search for “ALE” (as a whole word, to avoid matching substrings of other words) in the NeuroVault collection title/description or the associated paper’s PubMed title/abstract.

2.4 Imaging data download

Here, we download Nifti images from each NeuroVault collection after applying the metadata exclusions. After preparing local file paths, it checks for missing files and downloads them in parallel using 36 cores.

3 Image derived exclusionary criteria

After downloading the images, we compute image dimensions (imgdim), voxel dimensions (voxdim), number of nonzero voxels (nvox), image intensity range (range), and orientation from the Nifti images. The image and voxel dimensions are used to compute the physical dimensions of each image in millimeters. We screen for duplicate images based on both the qualitative and quantitative aspects, including the file’s basename imgdim, voxeldim, nvox, range, and bounding_box. It would be rare for two different images to have exact identical values across all these metrics, especially bounding_box and range. After identifying duplicates, we check whether the image dimensions in mm are within 2.5% of the MNI template image. Finally, we check the image range; if it is less than 0.01, we exclude the image (Figure 2). This threshold was chosen because for images where the maximum and minimum intensity values differ by less than 0.01, they contain effectively no signal variation, and are therefore unlikely to represent valid statistical maps. We do not exclude based on large ranges because the test statistic under the alternative converges to infinity as the sample size gets large. Further effect size based exclusionary criteria are applied below. Tables 4 and 5 show the dataset before and after these exclusionary criteria, respectively.

3.1 Duplicates handling

# Check if image dimensions/voxels metadata was already computed this month
rds_file_path <- get_monthly_rds_path("imageMetadataWithDimensions")
file_exists_this_month <- file.exists(rds_file_path)

if (!file_exists_this_month) {
  # Computes image dimensions, saves them internally, and returns the updated data frame
  compute_and_save_image_dimensions_fn <- function(rds_dir, images, n_cores, c3d_bin) {
    library(parallel)
    
    # Helper function to get monthly RDS path
    get_monthly_rds_path <- function(base_name, dir = rds_dir) {
      get_date_suffix <- function() {
        format(Sys.Date(), "%Y-%m")
      }
      date_suffix <- get_date_suffix()
      current_file <- file.path(dir, paste0(base_name, "_", date_suffix, ".rds"))
      
      pattern <- paste0("^", base_name, "_[0-9]{4}-[0-9]{2}\\.rds$")
      all_files <- list.files(dir, pattern = pattern, full.names = TRUE)
      
      old_files <- all_files[all_files != current_file]
      if (length(old_files) > 0) {
        invisible(lapply(old_files, file.remove))
      }
      
      return(current_file)
    }
    
    # Helper function to extract value by keyword from a string
    extract_info_field <- function(lines, keyword) {
      val <- sapply(lines, function(x) {
        first_line <- x[1]
        splt <- strsplit(first_line, split = ";")[[1]]
        match <- grep(keyword, splt, value = TRUE)
        if (length(match) > 0) {
          return(trimws(sub(paste0("^.*", keyword, "\\s*=\\s*"), "", match[1])))
        } else {
          return(NA)
        }
      })
      return(val)
    }
    
    # Get image dimensions and voxel info
    cat(sprintf("[%s] [BG] Computing image dimensions with c3d -info...\n", Sys.time()))
    cmd = paste(c3d_bin, images$localfile, '-info')
    voxdim = mclapply(cmd, system, intern=TRUE, mc.cores=n_cores)
    images$voxdim = unlist(lapply(voxdim, function(x){splt = strsplit(x[1], split = ';')[[1]]; splt[grep('vox', splt)] }))
    images$imgdim = unlist(lapply(voxdim, function(x){splt = strsplit(x[1], split = ';')[[1]]; splt[grep('dim', splt)] }))
    images$voxdim <- trimws(images$voxdim)
    images$imgdim <- trimws(images$imgdim)
    cat(sprintf("[%s] [BG] Image dimensions computed\n", Sys.time()))
    
    # Count non-zero voxels
    cat(sprintf("[%s] [BG] Counting non-zero voxels with c3d -threshold...\n", Sys.time()))
    cmd_vox <- paste(c3d_bin, images$localfile, '-threshold 0 0 0 1 -voxel-sum')
    nvox <- unlist(mclapply(cmd_vox, system, intern = TRUE, mc.cores = n_cores))
    images$nvox <- as.numeric(gsub("[^0-9.]", "", nvox))
    cat(sprintf("[%s] [BG] Non-zero voxel count complete\n", Sys.time()))
    
    # Orientation, range, and bounding box are parsed from the same `-info`
    # output already fetched above (`voxdim`) -- avoids a second c3d -info pass
    images$orient <- extract_info_field(voxdim, "orient")
    images$range <- extract_info_field(voxdim, "range")
    images$bounding_box <- extract_info_field(voxdim, "bb")
    images$voxdim <- trimws(gsub("vox = ", "", images$voxdim))
    images$imgdim <- trimws(gsub("Image #1: dim = ", "", images$imgdim))
    
    # Save to versioned RDS file
    if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
    rds_path <- get_monthly_rds_path("imageMetadataWithDimensions")
    
    cat(sprintf("[%s] Saving enhanced metadata to %s\n", Sys.time(), rds_path))
    saveRDS(images, file = rds_path, compress = "gzip")
    
    cat(sprintf("[%s] ✓ Image dimension metadata saved successfully\n", Sys.time()))
    return(images)
  }
  
  # Create logfiles directory if needed
  logdir <- "logfiles"
  if (!dir.exists(logdir)) dir.create(logdir, showWarnings = FALSE)
  logfiles = c('logfiles/c3d_dimensions.log', 'logfiles/c3d_dimensions.err')
  unlink(logfiles)
  
  images <- run_with_logging(
    compute_and_save_image_dimensions_fn(rds_dir, images, n_cores, c3d_bin),
    logfiles[1], logfiles[2]
  )
  
  message(sprintf("[%s] Computed image dimensions. Output: %s", 
            Sys.time(), logfiles[1]))
} else {
  # Load already computed metadata
  images <- load_monthly_rds("imageMetadataWithDimensions")
  message("Loaded previously computed image dimensions from ", rds_file_path)
}

Table 3: Duplicate images identified and removed based on identical intensity range and filename.

images$localfile_base <- basename(as.character(images$localfile))
cand_keys <- c("range", "localfile_base")

keyDF <- as.data.frame(lapply(images[, cand_keys], function(x) {
  if (is.list(x)) sapply(x, paste, collapse = ",") else as.character(x)
}))

images$duplicated <- duplicated(keyDF)
n_dup <- sum(images$duplicated)

reactable(
  images[images$duplicated == TRUE, c('collection', 'id', 'range', 'localfile_base')],
  striped = TRUE,
  highlight = TRUE,
  bordered = TRUE,
  searchable = TRUE,
  defaultPageSize = 10,
  columns = list(
    collection     = colDef(name = "Collection", minWidth = 150),
    id             = colDef(name = "ID",         minWidth = 80),
    range          = colDef(name = "Range",      minWidth = 150),
    localfile_base = colDef(name = "Filename",   minWidth = 200)
  )
)

images <- images[!images$duplicated, ]
if (exists("keyDF")) rm(keyDF)

We screen for duplicate images based on identical intensity range and filename, as two images sharing both values are almost certainly the same statistical map uploaded to multiple collections. A total of 2066 duplicate images were identified and removed, leaving 31068 images for downstream analysis. The duplicates are displayed in the table above. The table above shows an interactive table of the objects being dropped here. Note: Only the second occurrence of each duplicate pair is shown here; the first occurrence is retained in the dataset and will not appear in this table, so it is normal for no two rows in the table to look identical to each other.

3.2 Proportionality checking

Using a proportionality tolerance of 0.025, we then check whether the physical dimensions of each image are approximately proportional to the MNI template (182x218x182 mm) as shown in Figure 2.

# compute dim_mm by multiplying imgdim with voxdim after converting to numeric values
imgdims = lapply(gsub('\\]', ')', gsub('\\[', 'c(', as.character(images$imgdim))), function(x) eval(parse(text=x)))
voxdims = lapply(gsub('\\]', ')', gsub('\\[', 'c(', as.character(images$voxdim))), function(x) eval(parse(text=x)))
dims_mat = t(mapply(function(x,y){ x*y}, imgdims, voxdims))
colnames(dims_mat) <- c("dim_x", "dim_y", "dim_z")

# Split dim_mm into x, y, z components
# SNV: dim_mm was not defined and only existed because it was loading an old dataset version.
images$dim_mm = apply(dims_mat, 1, function(x) paste(x, collapse = "*"))


# Compute ratios
ratio_xy <- dims_mat[,1] / dims_mat[,2]
ratio_xz <- dims_mat[,1] / dims_mat[,3]
ratio_yz <- dims_mat[,2] / dims_mat[,3]


# Define reference ratios (from MNI: 182x218x182)
ref_xy <- 182 / 218       # ≈ 0.8349
ref_xz <- 1               # since x = z
ref_yz <- 218 / 182       # ≈ 1.1978

images[, paste('log10', c('ratio_xy', 'ratio_xz', 'ratio_yz'), sep = '_')] <- cbind(ratio_xy/ref_xy, ratio_xz/ref_xz, ratio_yz/ref_yz)

images[, c('ratio_x', 'ratio_y', 'ratio_z') ] <- sweep(dims_mat, 2, c(182, 218, 182), FUN = "/")


# Reference thresholds by dimension
ref_thresholds <- list(
  ratio_x = c(0.75, 1.25),
  ratio_y = c(0.75, 1.25),
  ratio_z = c(0.6, 1.2)
)

images$is_proportional = apply(images[, c('ratio_x', 'ratio_y', 'ratio_z')], 1, function(r) {
  all(r >= ref_thresholds$ratio_x[1] & r <= ref_thresholds$ratio_x[2]) &
    all(r >= ref_thresholds$ratio_y[1] & r <= ref_thresholds$ratio_y[2]) &
    all(r >= ref_thresholds$ratio_z[1] & r <= ref_thresholds$ratio_z[2])
})

# Logical mask for proportional dimensions
# Set epsilon
epsilon <- 0.025
# images$is_proportional <- abs(ratio_xy - ref_xy) < epsilon &
#                    abs(ratio_xz - ref_xz) < epsilon &
#                    abs(ratio_yz - ref_yz) < epsilon

# SNV: Changed to not repeat code you had above and use the same (more compact) syntax you used above
images <- images %>% mutate(
  dim_mm      = factor(dim_mm,      levels = names(sort(table(dim_mm),      decreasing = TRUE)))
)

prop_near_mni <- sum(images$is_proportional)
mni_total <- length(images$is_proportional)


# PLOT

library(GGally)
library(ggplot2)
library(scales)
library(rlang)

# Create ratio dataframe
ratio_df <- images[, c('ratio_x', 'ratio_y', 'ratio_z')]
ratio_df <- na.omit(ratio_df)

# Diagonal function (histograms)
diag_fun_ratio <- function(data, mapping, ...) {
  p <- GGally::ggally_barDiag(data, mapping, ...)
  
  var <- as_label(mapping$x)
  thrs <- ref_thresholds[[var]]
  
  if (!is.null(thrs)) {
    for (thr in thrs) {
      p <- p + geom_vline(xintercept = thr, linetype = "dashed", color = "red", alpha = 0.7)
    }
  }
  p
}

# Lower function (scatter plots)
lower_fun_ratio <- function(data, mapping, ...) {
  p <- GGally::ggally_points(data, mapping, ...)
  
  xvar <- as_label(mapping$x)
  yvar <- as_label(mapping$y)
  
  x_thrs <- ref_thresholds[[xvar]]
  y_thrs <- ref_thresholds[[yvar]]
  
  if (!is.null(x_thrs)) {
    for (thr in x_thrs) {
      p <- p + geom_vline(xintercept = thr, linetype = "dashed", color = "red", alpha = 0.7)
    }
  }
  
  if (!is.null(y_thrs)) {
    for (thr in y_thrs) {
      p <- p + geom_hline(yintercept = thr, linetype = "dashed", color = "red", alpha = 0.7)
    }
  }
  
  p
}

# Upper function (correlations)
upper_fun_ratio <- function(data, mapping, ...) {
  GGally::ggally_cor(data, mapping, ...) +
    theme_void()
}

# Create the plot
GGally::ggpairs(
  ratio_df,
  diag = list(continuous = diag_fun_ratio),
  lower = list(continuous = lower_fun_ratio),
  upper = list(continuous = upper_fun_ratio),
  title = "Figure 2: Image dimension ratios (relative to MNI 182×218×182 template)"
)

We find that 30665 of 31068 images are approximately proportional to the MNI template.

3.3 Range Exclusion

We convert the image dimensions to mm by multiplying the number of voxels by the pixel dimension. Using a proportionality tolerance of \(\varepsilon = 0.025\), we then check whether the physical dimensions of each image are approximately proportional to the MNI template (182×218×182 mm). We find that 30665 of 31068 images are approximately proportional to the MNI template.

Images with a low intensity range (\(|\max - \min|\)) are unlikely to be valid statistical images, so we exclude all images that have a range smaller than 0.01. A visual illustration is available in Figure 3 below. Everything below the red dashed line are removed due to small range.

library(stringr)
library(ggplot2)
library(scales)

# SNV: removed loading and saving dataset in this chunk

# parse [lo, hi]
hi <- str_match(
  images[["range"]],
  "\\[\\s*([+-]?(?:\\d*\\.?\\d+(?:[eE][+-]?\\d+)?))\\s*,\\s*([+-]?(?:\\d*\\.?\\d+(?:[eE][+-]?\\d+)?))\\s*\\]"
)
lo <- as.numeric(hi[, 2])
hi <- as.numeric(hi[, 3])

# attach & compute range
images$range_low  <- lo
images$range_high <- hi
images$range_diff <- images$range_high - images$range_low
rm(lo, hi)

# histogram of log10(range_diff) (positive values only)
xpos <- images$range_diff
# SNV: commented this because it was already satisfied.
# xpos <- xpos[is.finite(xpos) & xpos > 0]

# xpos already defined above as positive, finite range_diff values
df_sc <- data.frame(
  idx = seq_along(xpos),
  l10 = log10(xpos)
)

tick_min <- floor(min(df_sc$l10, na.rm = TRUE))
tick_max <- ceiling(max(df_sc$l10, na.rm = TRUE))
ticks    <- tick_min:tick_max

# SNV: changed hard coding of horizontal line
# SNV: Change this to -2, since it is thresholding on the log scale
thr <- -2
ggplot(df_sc, aes(idx, l10)) +
  geom_point(alpha = 0.5, size = 0.6) +
  scale_y_continuous(breaks = ticks, labels = scales::math_format(10^.x)) +
  geom_hline(yintercept = ticks, linetype = "dotted", linewidth = 0.3, color = "grey70") +
  labs(
    title = "Figure 3: Scatterplot of image range (log10 scale) by row index",
    x = "Image index",
    y = "Image numerical range (log10 scale)",
    subtitle = sprintf("n = %d ", length(xpos))
  ) +
  geom_hline(yintercept = thr,linetype = "dashed", color = "red") + 
  theme_minimal(base_size = 12)

rm(df_sc)

# drop rows with log10(range_diff) < 0.1 (guard against non-positive / NA)
images$range_exclusion = log10(images$range_diff) < thr

Table 4 below summarizes the dataset after image-level exclusions, showing the distribution of our dataset after the three additional criteria derived from the downloaded images: duplicated (identical range and filename across collections), is_proportional (physical dimensions within 2.5% of the MNI 182×218×182 template), and range_exclusion (intensity range below 0.01). The dim_mm column shows the distribution of image dimensions in millimeters across the dataset.

The remaining images are proportionally sized relative to the MNI template, contain meaningful intensity variation, and have no detected duplicates. The dim_mm distribution in Table 4 reflects only the image dimensions present in the retained dataset, with the most common dimension being 182218182.

# SNV: perc_of_voxels_outside was not included in table1Form, so I did not include it here.
table2Form = update.formula(table1Form, ~ duplicated + is_proportional + range_exclusion + dim_mm + .)
images <- images[which(images$is_proportional & !images$duplicated & !images$range_exclusion), ]
dimtab = table(images$dim_mm); dimtab = dimtab[dimtab > 0]
images$dim_mm = factor(images$dim_mm, levels = names(sort(dimtab,      decreasing = TRUE)))

table1(table2Form, data = images, caption = "Table 4: After image-level exclusionary criteria with prior metadata.")

Table 4: After image-level exclusionary criteria with prior metadata.
	Overall (N=30307)
duplicated
Yes	0 (0%)
No	30307 (100%)
is_proportional
Yes	30307 (100%)
No	0 (0%)
range_exclusion
Yes	0 (0%)
No	30307 (100%)
dim_mm
182218182	18131 (59.8%)
195231196	2329 (7.7%)
194230194	2064 (6.8%)
195231195	1174 (3.9%)
192228192	1124 (3.7%)
158190158	840 (2.8%)
195232.5195	650 (2.1%)
159189138	536 (1.8%)
181.5217.5181.5	469 (1.5%)
195231198	460 (1.5%)
196.749230.755194.4	391 (1.3%)
159189156	299 (1.0%)
183219183	298 (1.0%)
194230194.4	288 (1.0%)
150186160	108 (0.4%)
184.25220184	77 (0.3%)
157.5192.5157.5	76 (0.3%)
162192150	76 (0.3%)
180216180	63 (0.2%)
194.4230.4194.4	59 (0.2%)
158190156	56 (0.2%)
158190138	45 (0.1%)
194230195	45 (0.1%)
196.218231.894196.42	33 (0.1%)
170210145	30 (0.1%)
197233189	28 (0.1%)
193.2229.2193.2	26 (0.1%)
181217181	24 (0.1%)
180217.5180	20 (0.1%)
188220184	19 (0.1%)
158190168	18 (0.1%)
182217182	18 (0.1%)
157.5190.5157.5	17 (0.1%)
169.5205.5169.5	17 (0.1%)
148184156	16 (0.1%)
202242212	16 (0.1%)
196232197.6	15 (0.0%)
157.5190157.5	14 (0.0%)
162.4191.4150.8	14 (0.0%)
194.247230.346195	14 (0.0%)
195.8231195.8	14 (0.0%)
157.5190137.5	13 (0.0%)
182.5220182.5	13 (0.0%)
159195156	12 (0.0%)
194230195.8	12 (0.0%)
192.5229.25192.5	11 (0.0%)
196.875231.25198	11 (0.0%)
153.9189160.38	10 (0.0%)
154190164	10 (0.0%)
164192140	10 (0.0%)
196.552231.384195.25	10 (0.0%)
152184152	9 (0.0%)
161190.75150.5	9 (0.0%)
195.816231.636198.75	9 (0.0%)
156186147	8 (0.0%)
158190136.5	8 (0.0%)
194.4230.4196.56	8 (0.0%)
153189147	7 (0.0%)
156.6189158.4	6 (0.0%)
158190174	6 (0.0%)
195231197.12	6 (0.0%)
157.5189136.5	5 (0.0%)
159195159	5 (0.0%)
196.472232.412194.4	5 (0.0%)
198.75232.5197.6	5 (0.0%)
198231195	5 (0.0%)
157.5195171	4 (0.0%)
182218183	4 (0.0%)
193.5229.5195	4 (0.0%)
195231196.56	4 (0.0%)
195231198.8	4 (0.0%)
139.01492173.76865151.65264	3 (0.0%)
183249183	3 (0.0%)
192.5227.5192.5	3 (0.0%)
196.735231.77197.1	3 (0.0%)
144180146	2 (0.0%)
157189156	2 (0.0%)
159192156	2 (0.0%)
160192140	2 (0.0%)
164192160	2 (0.0%)
169205169	2 (0.0%)
171207171	2 (0.0%)
178.5214.5178.5	2 (0.0%)
182.5217.5182.5	2 (0.0%)
192.5230192.5	2 (0.0%)
193.952229.216193.6	2 (0.0%)
195.468230.373195.5	2 (0.0%)
195.966233.784197.2	2 (0.0%)
195232.5198	2 (0.0%)
200235191.4	2 (0.0%)
202232192	2 (0.0%)
150.5185158	1 (0.0%)
153.6189.6158.4	1 (0.0%)
154192158	1 (0.0%)
159192138	1 (0.0%)
159192159	1 (0.0%)
161.25191.25136.8	1 (0.0%)
161.25191.25140	1 (0.0%)
161191151	1 (0.0%)
180219184	1 (0.0%)
183219174	1 (0.0%)
183219183.6	1 (0.0%)
184.8220.8182.4	1 (0.0%)
188224184	1 (0.0%)
190.9227.7190.9	1 (0.0%)
192.5230.3125192.5	1 (0.0%)
193.5229.5193.5	1 (0.0%)
193.6228.8193.6	1 (0.0%)
193229193	1 (0.0%)
198234190	1 (0.0%)
number_of_subjects
Mean (SD)	67.5 (207)
Median [Min, Max]	30.0 [1.00, 30000]
modality
Diffusion MRI	65 (0.2%)
fMRI-BOLD	15140 (50.0%)
fMRI-CBF	14 (0.0%)
fMRI-CBV	1 (0.0%)
MEG	39 (0.1%)
Other	14662 (48.4%)
PET FDG	26 (0.1%)
PET other	42 (0.1%)
Structural MRI	318 (1.0%)
image_type
statistic_map	30307 (100%)
NIDM results statistic map	0 (0%)
atlas	0 (0%)
map_type
univariate-beta map	0 (0%)
variance	0 (0%)
Z map	2505 (8.3%)
T map	27802 (91.7%)
other	0 (0%)
ROI/mask	0 (0%)
P map (given null hypothesis)	0 (0%)
parcellation	0 (0%)
F map	0 (0%)
1-P map ("inverted" probability)	0 (0%)
multivariate-beta map	0 (0%)
anatomical	0 (0%)
Chi squared map	0 (0%)
analysis_level
single-subject	0 (0%)
group	30307 (100%)
meta-analysis	0 (0%)
other	0 (0%)
	0 (0%)
is_thresholded
FALSE	30307 (100%)
TRUE	0 (0%)
not_mni
FALSE	30307 (100%)
TRUE	0 (0%)
perc_bad_voxels
Mean (SD)	75.9 (9.49)
Median [Min, Max]	78.9 [0, 86.1]
perc_voxels_outside
Mean (SD)	9.33 (19.3)
Median [Min, Max]	0 [0, 100]
is_meta_analysis
Yes	0 (0%)
No	30307 (100%)

3.4 Image value processing

Aside from the above three approaches, there are other attributes of certain images that stops registration.

There could be NaN values in the images that causes FLIRT to fail. We want to remove them as well as some values that are very close to zero due to prior registration artifact.
There are images with active voxels everywhere, or with very extremely low intensity voxels. We set a fixed \(\varepsilon\) (0.001) here, to filter out voxels with values around zero by \(\varepsilon\) in all images.

This value was chosen by eroding the non-brain mask, defined below, and looking at a histogram of the 90th quantile of the values outside the brain.

# Define the image processing function
process_images_c3d_fn <- function(images_data, n_cores, eps_val) {
  library(parallel)
  
  # Create output filenames
  images_data$localfileadj <- sub("\\.nii(\\.gz)?$", "_eps.nii.gz", 
                                   images_data$localfile, ignore.case = TRUE)
  
  # Only (re)process images whose output doesn't already exist, so re-renders
  # don't redo this step for files already processed in a previous run.
  needs_processing <- !file.exists(images_data$localfileadj) | is.na(file.size(images_data$localfileadj)) | file.size(images_data$localfileadj) == 0
  cat(sprintf("[%s] Processing %d of %d images (skipping %d already processed)...\n",
              Sys.time(), sum(needs_processing), nrow(images_data), sum(!needs_processing)))
  
  localfile_todo    <- images_data$localfile[needs_processing]
  localfileadj_todo <- images_data$localfileadj[needs_processing]
  
  if (length(localfile_todo) > 0) {
    # Combined into a single c3d invocation (NaN-replace + threshold + multiply)
    # instead of two, so each file is only read from and written to disk once.
    cat(sprintf("[%s] Converting NaN to zeros and applying threshold filter (eps = %g)...\n", Sys.time(), eps_val))
    cmd <- sprintf("%s %s -replace NaN 0.0 -dup -threshold -%g %g 0 1 -multiply -o %s",
                   shQuote("/usr/local/c3d-1.1.0-Linux-x86_64/bin/c3d"),
                   shQuote(localfile_todo), eps_val, eps_val,
                   shQuote(localfileadj_todo))
    invisible(mclapply(cmd, function(c) system(paste(c, ">/dev/null")), 
                       mc.cores = n_cores))
    cat(sprintf("[%s] NaN conversion and threshold filtering complete\n", Sys.time()))
  }
  
  # Return updated filenames
  return(images_data$localfileadj)
}

# Create logfiles directory if needed
logdir <- "logfiles"
if (!dir.exists(logdir)) dir.create(logdir, showWarnings = FALSE)
logfiles = c('logfiles/image_processing_c3d.log', 'logfiles/image_processing_c3d.err')
unlink(logfiles)

# Run and wait for completion
eps <- 1e-3
cat(sprintf("[%s] Starting image processing (%d images)...\n", Sys.time(), nrow(images)))

## [2026-09-30 11:11:22.194307] Starting image processing (30307 images)...

localfileadj_result <- run_with_logging(
  process_images_c3d_fn(images, n_cores, eps),
  logfiles[1], logfiles[2]
)

# Update images dataframe with processed filenames
images$localfileadj <- localfileadj_result
cat(sprintf("[%s] ✓ Image processing completed. Output: %s\n", Sys.time(), logfiles[1]))

## [2026-09-30 11:11:22.737393] ✓ Image processing completed. Output: logfiles/image_processing_c3d.log

4 Image realignment and Registration

Despite the images being reported as MNI space. After visual inspection, we noticed they were not properly aligned with the MNI template. Examples of some misaligned images from the image dimension categories shown in Table 4 are shown in Figure 3.

Figure 4: Examples of image misalignment.

To resolve image misalignment, we evaluated registration of varying complexity:

Reslicing – no registration is performed; the image is simply resampled into the MNI 2mm grid.
Rigid registration (DOF = 6) – allows only rotation and translation, preserving the original shape and scale of the brain.
Affine registration (DOF = 12) – adds scaling and shearing to the rigid transform.

For both rigid and linear registrations, we use mutual information as the cost function (-cost mutualinfo using FLIRT) to accommodate differences in contrast and intensity between the template and statistical images. Results are shown for the seven most common image dimensions covering the majority of the dataset.

4.1 Registration snapshots

4.1.1 158x190x158

4.1.2 182x218x182

4.1.3 192x228x192

4.1.4 194x230x194

4.1.5 195x231x195

4.1.6 195x231x196

4.1.7 195x232.5x195

Based on these results, we chose rigid registration, as a balance between flexibility and accuracy across the image dimension categories. With the registered images, we further extracted information from the data:

voxdim_new: indicates the new voxel dimension of the images, now all images should be 2mm isotropic.
imgdim_new: indicates the new image dimension, now all images should be 91*109*91.
range_new: indicates the new range of voxel intensity.

4.2 Run registration

library(RNifti)

# Add resliced output filenames to dataframe (fast operation, no need for background)
cat(sprintf("[%s] Adding resliced output filenames\n", Sys.time()))

## [2026-09-30 11:11:22.868706] Adding resliced output filenames

images$resliced <- gsub('.nii.gz', '_rigid_cost_step1.nii.gz', images$localfile)

# Define the registration function
register_images_fn <- function(images_data, template_path, n_cores, outdir) {
  library(parallel)
  
  ## FSL Setup - Robust
  if (!nzchar(Sys.which("flirt"))) {
    if (!nzchar(Sys.getenv("FSLDIR"))) {
      fsl_candidates <- c("/usr/local/fsl", "/usr/share/fsl", "/opt/fsl")
      existing <- file.exists(fsl_candidates)
      if (any(existing)) {
        Sys.setenv(FSLDIR = fsl_candidates[which(existing)[1]])
      } else {
        stop("FSL not found. Please install FSL and set FSLDIR environment variable")
      }
    }
    Sys.setenv(FSLOUTPUTTYPE = "NIFTI_GZ")
    fslbin <- file.path(Sys.getenv("FSLDIR"), "bin")
    Sys.setenv(PATH = paste(fslbin, Sys.getenv("PATH"), sep = ":"))
  }
  
  ## Reference image
  ref <- file.path(Sys.getenv("FSLDIR"), "data", "standard", "MNI152_T1_2mm_brain.nii.gz")
  
  ## Helpers
  validate_flirt <- function() {
    flirt <- Sys.which("flirt")
    if (!nzchar(flirt)) stop("FLIRT not found in PATH")
    if (system2(flirt, args = "-version", stdout = FALSE, stderr = FALSE) != 0) {
      stop("FLIRT exists but failed to run. Please check FSL installation")
    }
    flirt
  }
  
  flirt <- validate_flirt()
  
  reslice_image <- function(input, output, ref) {
    if (!file.exists(output)) {
      message(sprintf("[%s] Reslicing: %s", Sys.time(), basename(input)))
      cmd <- sprintf("%s -in '%s' -ref '%s' -out '%s' -dof 6 -cost mutualinfo -interp trilinear",
                     flirt, input, ref, output)
      status <- system(cmd, ignore.stdout = TRUE, ignore.stderr = TRUE)
      if (status != 0 || !file.exists(output) || file.size(output) == 0) {
        warning("FLIRT failed for: ", input)
        return(NA_character_)
      }
    }
    return(output)
  }
  
  ## Reslice all images
  cat(sprintf("[%s] Starting reslicing with %d cores...\n", Sys.time(), n_cores))
  results <- mcmapply(
    FUN = reslice_image,
    input = images_data$localfileadj,
    output = images_data$resliced,
    MoreArgs = list(ref = ref),
    mc.cores = n_cores
  )
  
  cat(sprintf("[%s] Image registration completed\n", Sys.time()))
}

# Create logfiles directory if needed
logdir <- "logfiles"
if (!dir.exists(logdir)) dir.create(logdir, showWarnings = FALSE)
logfiles = c('logfiles/flirt_registration.log', 'logfiles/flirt_registration.err')
unlink(logfiles)

run_with_logging(
  register_images_fn(images, template, n_cores, outdir),
  logfiles[1], logfiles[2]
)

message(sprintf("[%s] Image registration completed. Output: %s", 
                Sys.time(), logfiles[1]))

library(stringr)

# Check if post-registration extraction and coverage data already exists this month
rds_file_path <- get_monthly_rds_path("imageExtractionAndCoverage")
file_exists_this_month <- file.exists(rds_file_path)

if (!file_exists_this_month) {
  # ============================================================================
  # PART 1: FEATURE EXTRACTION FROM REGISTERED IMAGES
  # ============================================================================
  
  # Count non-zero voxels with error handling
  cmd_vox <- paste(c3d_bin, images$resliced, '-threshold 0 0 0 1 -voxel-sum')
  nvox_result <- tryCatch({
    unlist(mclapply(cmd_vox, function(cmd) {
      tryCatch(system(cmd, intern = TRUE), error = function(e) "")
    }, mc.cores = n_cores))
  }, error = function(e) {
    cat(sprintf("[%s] Warning: c3d voxel count failed\n", Sys.time()))
    rep("", length(cmd_vox))
  })
  images$nvox_new <- as.numeric(gsub("[^0-9.]", "", nvox_result))
  images$nvox_new[is.na(images$nvox_new)] <- NA
  
  # Extract new info with c3d -info
  cat(sprintf("[%s] Extracting image info (voxdim, imgdim, orientation, range)...\n", Sys.time()))
  cmd_info <- paste(c3d_bin, images$resliced, "-info")
  info_result <- tryCatch({
    mclapply(cmd_info, function(cmd) {
      tryCatch(system(cmd, intern = TRUE), error = function(e) "")
    }, mc.cores = n_cores)
  }, error = function(e) {
    cat(sprintf("[%s] Warning: c3d -info failed\n", Sys.time()))
    rep(list(""), length(cmd_info))
  })
  
  # Enhanced extractField function with error handling
  extractField <- function(x, field_regex) {
    tryCatch({
      if (length(x) == 0 || all(x == "")) return(NA_character_)
      line <- x[1]
      if (is.na(line) || line == "") return(NA_character_)
      match <- regmatches(line, regexpr(field_regex, line))
      if (length(match) == 0) return(NA_character_)
      return(match)
    }, error = function(e) NA_character_)
  }
  
  # Extract fields with safe mapping
  cat(sprintf("[%s] Parsing extracted fields...\n", Sys.time()))
  images$voxdim_new <- sapply(info_result, extractField, "vox = \\[.*?\\]")
  images$imgdim_new <- sapply(info_result, extractField, "dim = \\[.*?\\]")
  images$orient_new <- sapply(info_result, extractField, "orient = [A-Z]+")
  images$range_new  <- sapply(info_result, extractField, "range = \\[.*?\\]")
  
  # Count how many extractions were successful
  na_count <- sum(is.na(images$voxdim_new))
  cat(sprintf("[%s] Extraction complete: %d successful, %d with missing/failed output\n", 
              Sys.time(), nrow(images) - na_count, na_count))
  
  # ============================================================================
  # PART 2: COMPUTE COVERAGE STATISTICS
  # ============================================================================
  
  # Get files for coverage computation
  files_all <- images$resliced
  keep <- !is.na(files_all) & file.exists(files_all)
  files <- files_all[keep]
  idx   <- which(keep)
  
  # Function to compute mask overlap statistics for all masks in one pass per
  # file, instead of re-reading every image from disk once per mask (4x fewer reads)
  gm_idx      <- which(readNifti(gm_mask) == 1)
  wm_idx      <- which(readNifti(wm_mask) == 1)
  outside_idx <- which(readNifti(nonbrain_mask) == 1)
  eroded_idx  <- which(readNifti(nonbrain_eroded_mask) == 1)
  
  coverage_stats <- mclapply(files, function(file){
    img <- readNifti(file)
    c(
      nnz_gm      = sum(img[gm_idx] != 0),
      nnz_wm      = sum(img[wm_idx] != 0),
      nnz_outside = sum(img[outside_idx] != 0),
      q90_outside = unname(quantile(abs(img[eroded_idx]), probs = 0.9))
    )
  }, mc.cores = n_cores)
  coverage_mat <- do.call(rbind, coverage_stats)
  
  # Write back to images (order-preserving; NAs for missing files)
  images$nvox_outside[idx] <- coverage_mat[, "nnz_outside"]
  images$nvox_gm[idx]      <- coverage_mat[, "nnz_gm"]
  images$nvox_wm[idx]      <- coverage_mat[, "nnz_wm"]
  images$q90_outside[idx]  <- coverage_mat[, "q90_outside"]
  
  # Save to versioned RDS file
  if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
  rds_file_path <- get_monthly_rds_path("imageExtractionAndCoverage")
  
  cat(sprintf("[%s] Saving extracted features and coverage data to %s\n", Sys.time(), rds_file_path))
  saveRDS(images, file = rds_file_path, compress = "gzip")
  
  cat(sprintf("[%s] ✓ Extraction and coverage data saved successfully\n", Sys.time()))
} else {
  # Load previously computed data
  images <- load_monthly_rds("imageExtractionAndCoverage")
  message("Loaded previously computed extraction and coverage data from ", rds_file_path)
}

## Loaded previously computed extraction and coverage data from /media/disk2/ree/metadataset/rds/imageExtractionAndCoverage_2026-09.rds

# Plot histogram of q90_outside for reference
#hist(log10(images$q90_outside[images$q90_outside>0]), xlab='log10(90th quantile)', main='90th image quantile')

4.3 Exclusionary criteria based on template overlap

After registration, we compute the proportion of nonzero voxels overlapping with three spatial masks derived from the MNI template: gray matter, white matter, and the region outside both (i.e., non-brain voxels). These coverage metrics serve as a proxy for registration quality. A well-registered brain image should have high overlap with gray and white matter and minimal signal outside the brain mask.

Figure 5 shows the distribution of these three coverage proportions across all registered images. Based on these histograms, we selected thresholds that separate the bulk of the distribution from a low-coverage tail, and visually evaluated excluded images to confirm the thresholds were not discarding well-registered data.

# thresholds
GM_MIN  <- 100000
WM_MIN  <- 20000
OUT_MAX <- 100000

GM_total <- sum(readNifti(gm_mask))
WM_total <- sum(readNifti(wm_mask))

out_total <- 91*109*91 - GM_total - WM_total
images$gm_perc <- images$nvox_gm / GM_total
images$wm_perc <- images$nvox_wm / WM_total
images$outside_perc <- images$nvox_outside / out_total


library(GGally)
library(ggplot2)
library(scales)
library(rlang)

# --- proportions already computed above ---
coverage_df <- images[, c("gm_perc", "wm_perc", "outside_perc")]
coverage_df <- na.omit(coverage_df)

colnames(coverage_df) <- c("Gray Matter", "White Matter", "Outside Mask")

# convert to percent for plotting
coverage_pct <- coverage_df * 100

# convert voxel-count thresholds to percent thresholds
ref_pct <- c(
  "Gray Matter"  = 100 * GM_MIN  / GM_total,
  "White Matter" = 100 * WM_MIN  / WM_total,
  "Outside Mask" = 100 * OUT_MAX / out_total
)

# formatter for axes
pct_lab <- scales::label_percent(scale = 1, accuracy = 1)  # since data already in 0-100

diag_fun <- function(data, mapping, ...) {
  p <- GGally::ggally_barDiag(data, mapping, ...)
  var <- as_label(mapping$x)
  thr <- ref_pct[[var]]
  if (!is.null(thr)) p <- p + geom_vline(xintercept = thr, linetype = "dashed")
  p + scale_x_continuous(labels = pct_lab)
}

lower_fun <- function(data, mapping, ...) {
  p <- GGally::ggally_points(data, mapping, ...)

  xvar <- as_label(mapping$x)
  yvar <- as_label(mapping$y)

  # Reference identity line y = x
  p <- p + geom_abline(
    slope = 1, intercept = 0,
    linetype = "solid", color = "black"
  )

  # Optional reference thresholds
  if (!is.null(ref_pct[[xvar]]))
    p <- p + geom_vline(xintercept = ref_pct[[xvar]], linetype = "dashed")

  if (!is.null(ref_pct[[yvar]]))
    p <- p + geom_hline(yintercept = ref_pct[[yvar]], linetype = "dashed")

  p +
    scale_x_continuous(labels = pct_lab) +
    scale_y_continuous(labels = pct_lab)
}

upper_fun <- function(data, mapping, ...) {
  GGally::ggally_cor(data, mapping, ...) +
    theme_void()
}

GGally::ggpairs(
  coverage_pct,
  title = "Figure 5: Spatial Coverage Percentage",
  upper = list(continuous = "cor"),
  diag  = list(continuous = diag_fun),
  lower = list(continuous = wrap(lower_fun, alpha = 0.35, size = 0.6))
)

We apply the following three exclusions:

Gray matter coverage above 55% — images with fewer than this proportion of nonzero voxels inside gray matter are likely missing substantial cortical signal, possibly due to partial brain coverage or failed registration.
White matter coverage above 32% — low white matter coverage similarly suggests the image does not span the full brain volume.
Non-brain signal below 15% — images with a high proportion of nonzero voxels outside the brain mask likely contain registration artifacts or were not properly skull-stripped prior to upload.

Together, these criteria remove images that are unlikely to represent valid whole-brain statistical maps in MNI space, regardless of whether they passed the earlier metadata-based filters.

# keep ONLY rows that pass all three checks (and are finite)
images$gm_exclude = images$nvox_gm <= GM_MIN
images$wm_exclude = images$nvox_wm <= WM_MIN
images$nvox_outside_exclude = (images$nvox_outside >= OUT_MAX)

# SNV: Hoazheng, to be consistent with the prior sections, lets use a table structure to show the exclusionary criteria
table3Form = update.formula(table2Form, ~ gm_exclude + wm_exclude + nvox_outside_exclude + .)
table1(table3Form, data=images,caption = "Table 5: Coverage-based exclusionary criteria prior to removal.")

Table 5: Coverage-based exclusionary criteria prior to removal.
	Overall (N=30307)
gm_exclude
Yes	213 (0.7%)
No	30094 (99.3%)
wm_exclude
Yes	78 (0.3%)
No	30229 (99.7%)
nvox_outside_exclude
Yes	787 (2.6%)
No	29520 (97.4%)
duplicated
Yes	0 (0%)
No	30307 (100%)
is_proportional
Yes	30307 (100%)
No	0 (0%)
range_exclusion
Yes	0 (0%)
No	30307 (100%)
dim_mm
182218182	18131 (59.8%)
195231196	2329 (7.7%)
194230194	2064 (6.8%)
195231195	1174 (3.9%)
192228192	1124 (3.7%)
158190158	840 (2.8%)
195232.5195	650 (2.1%)
159189138	536 (1.8%)
181.5217.5181.5	469 (1.5%)
195231198	460 (1.5%)
196.749230.755194.4	391 (1.3%)
159189156	299 (1.0%)
183219183	298 (1.0%)
194230194.4	288 (1.0%)
150186160	108 (0.4%)
184.25220184	77 (0.3%)
157.5192.5157.5	76 (0.3%)
162192150	76 (0.3%)
180216180	63 (0.2%)
194.4230.4194.4	59 (0.2%)
158190156	56 (0.2%)
158190138	45 (0.1%)
194230195	45 (0.1%)
196.218231.894196.42	33 (0.1%)
170210145	30 (0.1%)
197233189	28 (0.1%)
193.2229.2193.2	26 (0.1%)
181217181	24 (0.1%)
180217.5180	20 (0.1%)
188220184	19 (0.1%)
158190168	18 (0.1%)
182217182	18 (0.1%)
157.5190.5157.5	17 (0.1%)
169.5205.5169.5	17 (0.1%)
148184156	16 (0.1%)
202242212	16 (0.1%)
196232197.6	15 (0.0%)
157.5190157.5	14 (0.0%)
162.4191.4150.8	14 (0.0%)
194.247230.346195	14 (0.0%)
195.8231195.8	14 (0.0%)
157.5190137.5	13 (0.0%)
182.5220182.5	13 (0.0%)
159195156	12 (0.0%)
194230195.8	12 (0.0%)
192.5229.25192.5	11 (0.0%)
196.875231.25198	11 (0.0%)
153.9189160.38	10 (0.0%)
154190164	10 (0.0%)
164192140	10 (0.0%)
196.552231.384195.25	10 (0.0%)
152184152	9 (0.0%)
161190.75150.5	9 (0.0%)
195.816231.636198.75	9 (0.0%)
156186147	8 (0.0%)
158190136.5	8 (0.0%)
194.4230.4196.56	8 (0.0%)
153189147	7 (0.0%)
156.6189158.4	6 (0.0%)
158190174	6 (0.0%)
195231197.12	6 (0.0%)
157.5189136.5	5 (0.0%)
159195159	5 (0.0%)
196.472232.412194.4	5 (0.0%)
198.75232.5197.6	5 (0.0%)
198231195	5 (0.0%)
157.5195171	4 (0.0%)
182218183	4 (0.0%)
193.5229.5195	4 (0.0%)
195231196.56	4 (0.0%)
195231198.8	4 (0.0%)
139.01492173.76865151.65264	3 (0.0%)
183249183	3 (0.0%)
192.5227.5192.5	3 (0.0%)
196.735231.77197.1	3 (0.0%)
144180146	2 (0.0%)
157189156	2 (0.0%)
159192156	2 (0.0%)
160192140	2 (0.0%)
164192160	2 (0.0%)
169205169	2 (0.0%)
171207171	2 (0.0%)
178.5214.5178.5	2 (0.0%)
182.5217.5182.5	2 (0.0%)
192.5230192.5	2 (0.0%)
193.952229.216193.6	2 (0.0%)
195.468230.373195.5	2 (0.0%)
195.966233.784197.2	2 (0.0%)
195232.5198	2 (0.0%)
200235191.4	2 (0.0%)
202232192	2 (0.0%)
150.5185158	1 (0.0%)
153.6189.6158.4	1 (0.0%)
154192158	1 (0.0%)
159192138	1 (0.0%)
159192159	1 (0.0%)
161.25191.25136.8	1 (0.0%)
161.25191.25140	1 (0.0%)
161191151	1 (0.0%)
180219184	1 (0.0%)
183219174	1 (0.0%)
183219183.6	1 (0.0%)
184.8220.8182.4	1 (0.0%)
188224184	1 (0.0%)
190.9227.7190.9	1 (0.0%)
192.5230.3125192.5	1 (0.0%)
193.5229.5193.5	1 (0.0%)
193.6228.8193.6	1 (0.0%)
193229193	1 (0.0%)
198234190	1 (0.0%)
number_of_subjects
Mean (SD)	67.5 (207)
Median [Min, Max]	30.0 [1.00, 30000]
modality
Diffusion MRI	65 (0.2%)
fMRI-BOLD	15140 (50.0%)
fMRI-CBF	14 (0.0%)
fMRI-CBV	1 (0.0%)
MEG	39 (0.1%)
Other	14662 (48.4%)
PET FDG	26 (0.1%)
PET other	42 (0.1%)
Structural MRI	318 (1.0%)
image_type
statistic_map	30307 (100%)
NIDM results statistic map	0 (0%)
atlas	0 (0%)
map_type
univariate-beta map	0 (0%)
variance	0 (0%)
Z map	2505 (8.3%)
T map	27802 (91.7%)
other	0 (0%)
ROI/mask	0 (0%)
P map (given null hypothesis)	0 (0%)
parcellation	0 (0%)
F map	0 (0%)
1-P map ("inverted" probability)	0 (0%)
multivariate-beta map	0 (0%)
anatomical	0 (0%)
Chi squared map	0 (0%)
analysis_level
single-subject	0 (0%)
group	30307 (100%)
meta-analysis	0 (0%)
other	0 (0%)
	0 (0%)
is_thresholded
FALSE	30307 (100%)
TRUE	0 (0%)
not_mni
FALSE	30307 (100%)
TRUE	0 (0%)
perc_bad_voxels
Mean (SD)	75.9 (9.49)
Median [Min, Max]	78.9 [0, 86.1]
perc_voxels_outside
Mean (SD)	9.33 (19.3)
Median [Min, Max]	0 [0, 100]
is_meta_analysis
Yes	0 (0%)
No	30307 (100%)

## The following chunk finds out the ones that should be excluded, i.e. low coverage, and visualize them

par(mfrow = c(3, 4), mar = c(1, 1, 3.5, 1), pty = "m")

tempBG <- readNifti(template)

low_gm   <- images[which(images$gm_exclude == TRUE), ][seq_len(min(4, sum(images$gm_exclude,             na.rm = TRUE))), ]
low_wm   <- images[which(images$wm_exclude == TRUE), ][seq_len(min(4, sum(images$wm_exclude,             na.rm = TRUE))), ]
high_out <- images[which(images$nvox_outside_exclude == TRUE), ][seq_len(min(4, sum(images$nvox_outside_exclude, na.rm = TRUE))), ]

spot_check <- rbind(low_gm, low_wm, high_out)
spot_check$reason <- c(
  rep("Low gray matter",    nrow(low_gm)),
  rep("Low white matter",   nrow(low_wm)),
  rep("High outside mask",  nrow(high_out))
)

for (i in seq_len(nrow(spot_check))) {
  img <- readNifti(spot_check$resliced[i])
  
  pbj:::image.niftiImage(
    img, BGimg = tempBG, lo = FALSE, index = 35,
    limits = quantile(abs(img), probs = 0.01),
    plane  = 'sagittal'
  )
  
  title(
    main = paste0(
      spot_check$reason[i], "\n",
      "GM: ",      round(100 * spot_check$gm_perc[i],      1), "% | ",
      "WM: ",      round(100 * spot_check$wm_perc[i],      1), "% | ",
      "Outside: ", round(100 * spot_check$outside_perc[i], 1), "%"
    ),
    cex.main = 0.65,
    line     = 0.3
  )
}

Figure 6: Representative sagittal slices of excluded images grouped by exclusion reason. Top row: low gray matter coverage (images with sparse or misaligned cortical signal). Middle row: low white matter coverage (images with minimal subcortical signal, likely due to partial brain coverage or poor registration). Bottom row: excessive non-brain signal (images where signal extends well beyond the brain mask, indicating registration failure or absent skull-stripping).

n_before <- nrow(images)
images   <- images[which(!(images$gm_exclude | images$wm_exclude | images$nvox_outside_exclude)), ]
n_after  <- nrow(images)
cat(sprintf("Excluded %d images. Remaining: %d\n", n_before - n_after, n_after))

## Excluded 1000 images. Remaining: 29307

The following shows what the dataset looks like after excluding.

# SNV: Changed to coverage_pct to be consistent
# --- proportions already computed above ---
coverage_df <- images[, c("gm_perc", "wm_perc", "outside_perc")]
coverage_df <- na.omit(coverage_df)

colnames(coverage_df) <- c("Gray Matter", "White Matter", "Outside Mask")

# convert to percent for plotting
coverage_pct <- coverage_df * 100

GGally::ggpairs(
  coverage_pct,
  title = "Figure 7: Spatial Coverage Percentage",
  upper = list(continuous = "cor"),
  diag  = list(continuous = diag_fun),
  lower = list(continuous = wrap(lower_fun, alpha = 0.35, size = 0.6))
)

5 Effect size analysis and exclusionary criteria

5.1 Effect size estimation

Here, we convert nonzero values in the T and Z statistical images into the signed version of the robust effect size index (RESI; using RESI::t2S for T maps and RESI::z2S for Z maps), based on each metadata sample size reported with the images. We save the resulting effect size images and summary statistics within gray and white matter (mean and 0, 0.25, 0.5, 0.75, 1 quantiles).

library(RNifti)
library(RESI)

# Check if effect size computation was already done this month
rds_file_path <- get_monthly_rds_path("imageMetadataEffectSizes")
file_exists_this_month <- file.exists(rds_file_path)

if (!file_exists_this_month) {
  # Computes effect sizes, saves them internally, and returns the updated data frame
  compute_and_save_effect_sizes_fn <- function(images_data, gm_mask_path, wm_mask_path, 
                                              mask_2mm_path, n_cores, rds_dir) {
    library(RNifti)
    library(parallel)
    library(RESI)
    
    # Helper function to get monthly RDS path
    get_monthly_rds_path <- function(base_name, dir = rds_dir) {
      get_date_suffix <- function() {
        format(Sys.Date(), "%Y-%m")
      }
      date_suffix <- get_date_suffix()
      current_file <- file.path(dir, paste0(base_name, "_", date_suffix, ".rds"))
      
      pattern <- paste0("^", base_name, "_[0-9]{4}-[0-9]{2}\\.rds$")
      all_files <- list.files(dir, pattern = pattern, full.names = TRUE)
      
      old_files <- all_files[all_files != current_file]
      if (length(old_files) > 0) {
        invisible(lapply(old_files, file.remove))
      }
      
      return(current_file)
    }
    
    # Load masks once
    gm_mask_nifti <- RNifti::readNifti(gm_mask_path)
    wm_mask_nifti <- RNifti::readNifti(wm_mask_path)
    brain_mask_nifti <- RNifti::readNifti(mask_2mm_path)
    
    # Validate that required columns exist
    cat(sprintf("[%s] [BG] Starting effect size computation for %d images...\n", Sys.time(), nrow(images_data)))
    
    # Get indices of nonzero voxels once for efficiency
    gm_idx   <- which(gm_mask_nifti != 0)
    wm_idx   <- which(wm_mask_nifti != 0)
    brain_idx <- which(brain_mask_nifti != 0)
    
    # Worker function that takes ALL arguments explicitly (no lexical scoping)
    process_single_image <- function(idx, file_path, n_subj, map_type, gm_idx, wm_idx, brain_idx) {
      tryCatch({
        if (is.na(file_path) || !nzchar(file_path) || !file.exists(file_path)) return(NULL)
        if (is.na(n_subj) || is.na(map_type)) return(NULL)
        
        img <- RNifti::readNifti(file_path)
        
        # Only extract nonzero values
        nonzero_mask <- (img != 0) & !is.na(img)
        valid_idx <- which(nonzero_mask)
        
        if (length(valid_idx) == 0) return(NULL)
        
        valid_vals <- img[valid_idx]
        
        # Convert map_type to character to handle factors
        map_type_char <- as.character(map_type)
        
        # Convert effect sizes
        if (map_type_char == "T map") {
          converted <- RESI::t2S(valid_vals, rdf = n_subj - 1, n = n_subj)
        } else if (map_type_char == "Z map") {
          converted <- RESI::z2S(valid_vals, n = n_subj, unbiased = TRUE)
        } else {
          return(NULL)
        }
        
        # Build effect size image only with converted values
        effect_img <- img
        effect_img[valid_idx] <- converted
        effect_img[!is.finite(effect_img)] <- 0
        
        output_dir <- dirname(file_path)
        if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
        base_name  <- sub("\\.nii\\.gz$", "", basename(file_path))
        out_file   <- file.path(output_dir, paste0(base_name, "_effect_size.nii.gz"))
        RNifti::writeNifti(effect_img, out_file)
        
        # Get stats from converted values only (not from full image)
        effect_only <- converted[converted != 0 & is.finite(converted)]
        es_stats <- if (length(effect_only) > 0) {
          c(
            quantile(effect_only, type = 7),
            es_mean = mean(effect_only)
          )
        } else rep(NA_real_, 6)
        
        # Compute masked statistics using indices only (no intermediate arrays)
        grey_vals  <- effect_img[gm_idx]
        white_vals <- effect_img[wm_idx]
        
        # Write brain-masked version
        brain_file <- sub("\\.nii\\.gz$", "_masked_mni.nii.gz", out_file)
        brain_masked <- effect_img
        brain_masked[-brain_idx] <- 0  # Zero out non-brain voxels
        RNifti::writeNifti(brain_masked,  brain_file)
        
        list(
          index       = idx,
          effect_size_file = out_file,
          brain_file = brain_file,
          es_min      = es_stats[1],
          es_1qt      = es_stats[2],
          es_median   = es_stats[3],
          es_3qt      = es_stats[4],
          es_max      = es_stats[5],
          es_mean     = es_stats[6],
          grey_mean   = if (length(grey_vals)  > 0) mean(grey_vals[grey_vals != 0], na.rm = TRUE)  else NA_real_,
          grey_min    = if (length(grey_vals)  > 0) min(grey_vals[grey_vals != 0], na.rm = TRUE)   else NA_real_,
          grey_max    = if (length(grey_vals)  > 0) max(grey_vals[grey_vals != 0], na.rm = TRUE)   else NA_real_,
          white_mean  = if (length(white_vals) > 0) mean(white_vals[white_vals != 0], na.rm = TRUE) else NA_real_,
          white_min   = if (length(white_vals) > 0) min(white_vals[white_vals != 0], na.rm = TRUE)  else NA_real_,
          white_max   = if (length(white_vals) > 0) max(white_vals[white_vals != 0], na.rm = TRUE)  else NA_real_
        )
      }, error = function(e) {
        message(" [", Sys.time(), "] Skipping index ", idx, ": ", e$message)
        NULL
      })
    }
    
    # Use mcmapply to pass only necessary vectors to workers (not entire environment)
    cat(sprintf("[%s] [BG] Starting effect size computation for %d images...\n", Sys.time(), nrow(images_data)))
    results <- mcmapply(
      process_single_image,
      idx = 1:nrow(images_data),
      file_path = images_data$resliced,
      n_subj = images_data$number_of_subjects,
      map_type = images_data$map_type,
      MoreArgs = list(gm_idx = gm_idx, wm_idx = wm_idx, brain_idx = brain_idx),
      mc.cores = n_cores,
      SIMPLIFY = FALSE
    )
    
    # Filter and bind
    cat(sprintf("[%s] Processing results...\n", Sys.time()))
    results <- Filter(Negate(is.null), results)
    df_results <- do.call(rbind, lapply(results, as.data.frame))
    
    # Add index column to images_data for matching
    images_data$index <- 1:nrow(images_data)
    
    # merge keeping all of images, matched by index
    images_data = merge(images_data, df_results, all.x=TRUE)
    
    # Remove the temporary index column (not in original)
    images_data$index <- NULL
    
    # Save directly to versioned RDS file
    if (!dir.exists(rds_dir)) dir.create(rds_dir, showWarnings = FALSE, recursive = TRUE)
    rds_path <- get_monthly_rds_path("imageMetadataEffectSizes")
    saveRDS(images_data, file = rds_path)
    
    cat(sprintf("[%s] Effect size computation completed. Results saved to %s\n", Sys.time(), rds_path))
    return(images_data)
  }
  
  logdir <- "logfiles"
  if (!dir.exists(logdir)) dir.create(logdir, showWarnings = FALSE)
  logfiles = c('logfiles/effectSize_estimation.log', 'logfiles/effectSize_estimation.err')
  unlink(logfiles)

  images <- run_with_logging(
    compute_and_save_effect_sizes_fn(images, gm_mask, wm_mask, mask_2mm_file, n_cores, rds_dir),
    logfiles[1], logfiles[2]
  )
  
  message(sprintf("[%s] Effect size computation saved to %s", 
            Sys.time(), rds_file_path))
} else {
  # Load already computed effect sizes
  images <- load_monthly_rds("imageMetadataEffectSizes")
  message("Loaded previously computed effect sizes from ", rds_file_path)
}

5.2 Meta-analysis outlier detection model

To screen for outlier images, we model the expected variability of these mean effect sizes as a function of (i) study sample size and (ii) voxel coverage, and flag observations that fall outside the implied funnel. The mean value of the effect size, \(y_{ij}\), for image \(j\) in collection \(i\), is proportional to \(n_i^{-1} V_i^{-2}\), where \(n_i\) denotes the number of subjects (sample size) and \(V_{ij}\) denotes the number of nonzero voxels in image \(j\) (see preprint for details).

We fit an intercept-only linear mixed-effects model with a collection-specific random intercept:

\[y_{ij} = \mu + u_{i} + \varepsilon_{ij},\]

where \(\mu\) is the grand mean effect size and \(u_{i} \sim \mathcal{N}(0, \tau^2)\) captures between-collection heterogeneity.

To account for heteroskedasticity across images, we model the residual variance using a power-of-covariates structure implemented via the \(\texttt{varPower}\) variance function in \(\texttt{nlme}\), and combine multiple variance components using . Specifically, for study-level covariates \(N_{ij}\) (sample size) and \(V_{ij}\) (scaled voxel coverage), we assume

\[\mathrm{Var}(\varepsilon_{ij}) = \sigma^2 \, |N_{ij}|^{2\delta_1} \, |V_{ij}|^{2\delta_2},\]

so that residual variability decreases as \(N\) and/or \(V\) increase when \(\delta_1 < 0\) and/or \(\delta_2 < 0\).

The total per-observation variance incorporates both residual and between-collection variability:

\[\sigma_{ij}^2 = \tau^2 + \sigma^2 \, |N_{ij}|^{2\delta_1} \, |V_{ij}|^{2\delta_2}.\]

This quantity represents the model-implied variance of observed mean effect sizes as a function of sample size and voxel coverage.

We flag whole images as outliers if their deviation from the fitted grand mean exceeds a multiple of its modeled standard deviation:

\[\bigl| y_{ij} - \hat{\mu} \bigr| > k \, \sigma_{ij},\]

with \(k = 3\).

This criterion standardizes residuals across studies while accounting for heteroskedasticity and within-collection correlation. This outlier identification procedure was applied only to the gray matter mean effect sizes. To mitigate the effect of outliers on the model fit used for outlier identification, we windsorized the data at the 1st and 99th percentiles. We did this after noting that a few extreme image values had a huge effect on the residual variance. The same procedure was applied separately to gray matter and white matter mean effect sizes. These analyses serve both to validate the effect-size transformation and to provide a principled, model-based mechanism for identifying implausible images prior to downstream inference. While alternative filtering strategies were explored, we selected the mixed-effects variance modeling approach as the most statistically grounded and interpretable method for screening artifactual data points. A summary table consolidating key model estimates and variance components is provided below.

# 1. Remove missing es_mean first
images <- images[!is.na(images$es_mean), ]

# 2. Remove invalid nvox_new
images <- images[!is.na(images$nvox_new) & 
                 is.finite(images$nvox_new) & 
                 images$nvox_new > 0, ]

# 3. Compute voxel_scaled once
images$voxel_scaled <- images$nvox_new / 1000

cat("N images:", nrow(images), "\n")

## N images: 29272

cat("voxel_scaled range:", range(images$voxel_scaled), "\n")

## voxel_scaled range: 129.957 344.329

cat("N unique collections:", length(unique(images$collection_id)), "\n")

## N unique collections: 3554

# 4. Winsorize
winsorize_vec <- function(x, probs = c(0.01, 0.99)) {
  qs <- quantile(x, probs = probs, na.rm = TRUE, names = FALSE, type = 7)
  x[x < qs[1]] <- qs[1]
  x[x > qs[2]] <- qs[2]
  x
}
images$es_mean_winsorized <- winsorize_vec(images$es_mean)

# 5. Fit models
m_power_winsorized <- lme(
  es_mean_winsorized ~ 1,
  random  = ~ 1 | collection_id,
  weights = varComb(
    varPower(form = ~ number_of_subjects),
    varPower(form = ~ voxel_scaled)),
  data   = images,
  method = "REML"
)

m_power_original <- lme(
  es_mean ~ 1,
  random  = ~ 1 | collection_id,
  weights = varComb(
    varPower(form = ~ number_of_subjects),
    varPower(form = ~ voxel_scaled)),
  data   = images,
  method = "REML"
)

summary(m_power_winsorized)

## Linear mixed-effects model fit by REML
##   Data: images 
##         AIC       BIC   logLik
##   -28124.79 -28083.36 14067.39
## 
## Random effects:
##  Formula: ~1 | collection_id
##         (Intercept) Residual
## StdDev:  0.06295416 7844.239
## 
## Combination of variance functions: 
##  Structure: Power of variance covariate
##  Formula: ~number_of_subjects 
##  Parameter estimates:
##     power 
## -0.132612 
##  Structure: Power of variance covariate
##  Formula: ~voxel_scaled 
##  Parameter estimates:
##     power 
## -1.918934 
## Fixed effects:  es_mean_winsorized ~ 1 
##                  Value   Std.Error    DF  t-value p-value
## (Intercept) 0.00189248 0.001590739 25718 1.189686  0.2342
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -5.72721892 -0.44290170  0.02037559  0.48200506  5.51071923 
## 
## Number of Observations: 29272
## Number of Groups: 3554

summary(m_power_original)

## Linear mixed-effects model fit by REML
##   Data: images 
##        AIC      BIC    logLik
##   38797.67 38839.09 -19393.84
## 
## Random effects:
##  Formula: ~1 | collection_id
##         (Intercept) Residual
## StdDev:   0.6590946 96798495
## 
## Combination of variance functions: 
##  Structure: Power of variance covariate
##  Formula: ~number_of_subjects 
##  Parameter estimates:
##      power 
## -0.7718969 
##  Structure: Power of variance covariate
##  Formula: ~voxel_scaled 
##  Parameter estimates:
##     power 
## -3.028569 
## Fixed effects:  es_mean ~ 1 
##                  Value Std.Error    DF  t-value p-value
## (Intercept) 0.02511256 0.0121996 25718 2.058474  0.0396
## 
## Standardized Within-Group Residuals:
##           Min            Q1           Med            Q3           Max 
## -5.761098e+01 -1.413797e-01  1.273443e-04  1.531902e-01  5.570419e+01 
## 
## Number of Observations: 29272
## Number of Groups: 3554

m_power <- m_power_winsorized

The following chunk uses the outlier boundaries, as described above, from the winsorized model to make exclusions.

overall_effect <- as.numeric(fixef(m_power))                
sigma2 <- sigma(m_power)^2                                   
tau2   <- as.numeric(VarCorr(m_power)[1, "Variance"])    

# varPower exponents (delta_n, delta_v)
expParams <- intervals(m_power)$varStruct[, 2]
delta_n <- expParams[1]
delta_v <- expParams[2]

images$se_model <- sqrt(
  tau2 + sigma2 *
    (images$number_of_subjects ^ (2*delta_n)) *
    (images$voxel_scaled       ^ (2*delta_v))
)

# standardized residuals
images$res_std <- residuals(m_power, type = "pearson")

k <- 3
images$outlier <- abs(images$es_mean - overall_effect) > k * images$se_model

# SNV: These are all satisfied already in the dataset
# dat_plot <- subset(
#   images,
#   is.finite(number_of_subjects) & number_of_subjects > 0 &
#     is.finite(es_mean) & is.finite(se_model)
# )
n_out <- sum(images$outlier)
n_tot <- nrow(images)

grid <- data.frame(
  number_of_subjects = seq(1, 10000, length.out = 800)
)

# choose a constant voxel_scaled for the band (median is a sensible default)
grid$voxel_scaled <- median(images$voxel_scaled, na.rm = TRUE)

grid$se_model <- sqrt(
  tau2 + sigma2 *
    (grid$number_of_subjects ^ (2*delta_n)) *
    (grid$voxel_scaled       ^ (2*delta_v))
)

grid$upper <- overall_effect + k * grid$se_model
grid$lower <- overall_effect - k * grid$se_model
grid$log_n <- log10(grid$number_of_subjects)

library(plotly)

# SNV: dat_plot was the same as images, just replace it here
# dat_plot$se_model <- sqrt(
#   tau2 + sigma2 *
#     (dat_plot$number_of_subjects ^ (2 * delta_n)) *
#     (dat_plot$voxel_scaled       ^ (2 * delta_v))
# )
# dat_plot$upper_bound <- overall_effect + k * dat_plot$se_model
# dat_plot$lower_bound <- overall_effect - k * dat_plot$se_model
# 
# dat_plot$outlier <- with(dat_plot, es_mean > upper_bound | es_mean < lower_bound)

n_seq <- 10^seq(0, 4, length.out = 60)
v_seq <- seq(
  quantile(images$voxel_scaled, 0.05, na.rm = TRUE),
  quantile(images$voxel_scaled, 0.95, na.rm = TRUE),
  length.out = 60
)

grid3d <- expand.grid(number_of_subjects = n_seq, voxel_scaled = v_seq)

grid3d$se_model <- sqrt(
  tau2 + sigma2 *
    (grid3d$number_of_subjects ^ (2 * delta_n)) *
    (grid3d$voxel_scaled       ^ (2 * delta_v))
)

grid3d$upper <- overall_effect + k * grid3d$se_model
grid3d$lower <- overall_effect - k * grid3d$se_model


xg <- log10(n_seq) 
yg <- v_seq        

UP <- t(matrix(grid3d$upper, nrow = length(n_seq), ncol = length(v_seq)))
LO <- t(matrix(grid3d$lower, nrow = length(n_seq), ncol = length(v_seq)))
MU <- matrix(overall_effect, nrow = length(yg), ncol = length(xg))

fig8 <- plot_ly() |>
  add_markers(
    data   = subset(images, !outlier),
    x      = ~log10(number_of_subjects),
    y      = ~voxel_scaled,
    z      = ~es_mean,
    marker = list(size = 2, opacity = 0.25, color = "gray"),
    name   = "Inside funnel"
  ) |>
  add_markers(
    data   = subset(images, outlier),
    x      = ~log10(number_of_subjects),
    y      = ~voxel_scaled,
    z      = ~es_mean,
    marker = list(size = 4, color = "red"),
    name   = sprintf("Outliers (n = %d / %d)", n_out, n_tot)
  ) |>
  add_surface(
    x         = xg, y = yg, z = MU,
    opacity   = 0.35,
    colorscale = list(c(0,1), c("royalblue","royalblue")),
    showscale = FALSE,
    name      = "Overall effect (μ)"
  ) |>
  add_surface(
    x          = xg, y = yg, z = UP,
    opacity    = 0.20,
    colorscale = list(c(0,1), c("orange","orange")),
    showscale  = FALSE,
    name       = sprintf("Upper bound (μ + %d·SD)", k)
  ) |>
  add_surface(
    x          = xg, y = yg, z = LO,
    opacity    = 0.20,
    colorscale = list(c(0,1), c("orange","orange")),
    showscale  = FALSE,
    name       = sprintf("Lower bound (μ − %d·SD)", k)
  ) |>
  layout(
    title = list(
      text = paste0(
        "<b>Figure 8: Outlier Detection via Winsorized Variance-Power Meta-Analytic Model</b><br>",
        "<sup>Effect sizes flagged as outliers if |ES − μ| > ", k, "·SD(n, v); ",
        "SD modeled as a function of sample size and voxel coverage. ",
        "Red points (n = ", n_out, " of ", n_tot, " images) excluded from subsequent analyses.</sup>"
      ),
      x    = 0.5,
      xref = "paper",
      font = list(size = 13)
    ),
    scene = list(
      xaxis = list(
        title = "log\u2081\u2080(Sample size)",
        titlefont = list(size = 11)
      ),
      yaxis = list(
        title = "Voxel coverage (nonzero voxels / 1e5)",
        titlefont = list(size = 11)
      ),
      zaxis = list(
        title = "Mean effect size (ES)",
        titlefont = list(size = 11)
      ),
      camera = list(
        eye = list(x = 1.6, y = -1.6, z = 0.8)
      )
    ),
    legend = list(
      x = 0.01, y = 0.99,
      bgcolor = "rgba(255,255,255,0.7)",
      bordercolor = "gray",
      borderwidth = 1
    ),
    margin = list(t = 120)   # extra top margin for the title
  )

fig8

Figure \(\ref{fig:fig8}\) presents a three-dimensional visualization of the outlier detection procedure applied to the meta-analytic image database. Each point represents a single neuroimaging study, plotted along three axes: log-transformed sample size (\(\log_{10}(n)\)), voxel coverage scaled by \(1 \times 10^5\) (y-axis), and mean effect size (\(\bar{\theta}\), z-axis).

The blue surface represents the overall pooled effect estimate (\(\hat{\mu}\)) derived from the winsorized variance-power mixed-effects model, which accounts for heteroscedasticity in effect size variance as a joint function of sample size and voxel coverage. The orange surfaces define the upper and lower outlier boundaries, set at \(\hat{\mu} \pm k \cdot \hat{\sigma}(n, v)\) where \(k = 3\), and where \(\hat{\sigma}(n, v)\) is modeled as:

\[ \hat{\sigma}(n, v) = \sqrt{\tau^2 + \sigma^2 \cdot n^{2\delta_n} \cdot v^{2\delta_v}} \]

with \(\hat{\tau}^2 = 0.004\), \(\hat{\sigma}^2 = 6.1532079\times 10^{7}\), \(\hat{\delta}_n = -0.133\), and \(\hat{\delta}_v = -1.919\) estimated from the data.

Gray points fall within the acceptable range and are retained for subsequent analyses. Red points (\(n = 658\) of \(29272\) total images) fall outside the funnel boundaries and are flagged as outliers, indicating that their effect sizes deviate more than \(k = 3\) standard deviations from the pooled estimate given their sample size and voxel coverage characteristics. These images are excluded from all downstream analyses.

This approach is more principled than applying a fixed effect size threshold, as it explicitly acknowledges that smaller studies and those with lower voxel coverage are expected to show greater variability under the variance-power model, and calibrates the exclusion boundaries accordingly.

5.3 Model refit with cleaned dataset

With the above exclusion, a new dataset is created with the outliers as defined from the model above excluded. The following chunk refits the model with the new subset without Windsorization. Outliers are identified using the same procedure, but not excluded.

# Summary table for the final curated dataset using table1 for consistency
table1( ~ outlier + number_of_subjects + es_mean + voxel_scaled + gm_perc + wm_perc + outside_perc + grey_mean + white_mean, 
        data = images, 
        caption = "Table 6: Summary statistics for the final curated Neurovault dataset before outlier exclusions.")

Table 6: Summary statistics for the final curated Neurovault dataset before outlier exclusions.
	Overall (N=29272)
outlier
Yes	658 (2.2%)
No	28614 (97.8%)
number_of_subjects
Mean (SD)	66.6 (209)
Median [Min, Max]	30.0 [3.00, 30000]
es_mean
Mean (SD)	0.0206 (0.740)
Median [Min, Max]	0.00438 [-52.7, 45.9]
voxel_scaled
Mean (SD)	231 (40.7)
Median [Min, Max]	209 [130, 344]
gm_perc
Mean (SD)	0.865 (0.0922)
Median [Min, Max]	0.814 [0.555, 1.00]
wm_perc
Mean (SD)	0.913 (0.0813)
Median [Min, Max]	0.917 [0.342, 1.00]
outside_perc
Mean (SD)	0.0259 (0.0347)
Median [Min, Max]	0.00797 [0, 0.152]
grey_mean
Mean (SD)	0.0228 (0.848)
Median [Min, Max]	0.00614 [-64.9, 52.4]
white_mean
Mean (SD)	0.0164 (0.670)
Median [Min, Max]	0.00231 [-43.2, 38.2]

dat = images[ images$outlier==FALSE,]
m_power <- lme(
  es_mean ~ 1,
  random  = ~ 1 | collection_id,
  weights = varComb(
  varPower(form = ~ number_of_subjects),
  varPower(form = ~ voxel_scaled)),
  data    = dat,
  method  = "REML"
)

overall_effect <- as.numeric(fixef(m_power))                
sigma2 <- sigma(m_power)^2                                   
tau2   <- as.numeric(VarCorr(m_power)[1, "Variance"])    

# varPower exponents (delta_n, delta_v)
expParams <- intervals(m_power)$varStruct[, 2]
delta_n <- expParams[1]
delta_v <- expParams[2]

dat$se_model <- sqrt(
  tau2 + sigma2 *
    (dat$number_of_subjects ^ (2*delta_n)) *
    (dat$voxel_scaled       ^ (2*delta_v))
)

# standardized residuals (optional panel A)
dat$res_std <- residuals(m_power, type = "pearson")

k <- 3
dat$outlier <- abs(dat$es_mean - overall_effect) > k * dat$se_model

# dat_plot <- subset(
#   dat,
#   is.finite(number_of_subjects) & number_of_subjects > 0 &
#     is.finite(es_mean) & is.finite(se_model)
# )
# dat_sorted <- dat_plot[order(dat_plot$number_of_subjects), ]
# dat_clean  <- subset(dat_sorted, !outlier)
# n_out <- sum(dat_plot$outlier)
# n_tot <- nrow(dat_plot)
# dat_plot$se_model <- sqrt(
#   tau2 + sigma2 *
#     (dat_plot$number_of_subjects ^ (2 * delta_n)) *
#     (dat_plot$voxel_scaled       ^ (2 * delta_v))
# )
# dat_plot$upper_bound <- overall_effect + k * dat_plot$se_model
# dat_plot$lower_bound <- overall_effect - k * dat_plot$se_model
# dat_plot$outlier <- with(dat_plot, es_mean > upper_bound | es_mean < lower_bound)

grid <- data.frame(
  number_of_subjects = seq(1, 10000, length.out = 800)
)

# choose a constant voxel_scaled for the band (median is a sensible default)
grid$voxel_scaled <- median(dat$voxel_scaled, na.rm = TRUE)

grid$se_model <- sqrt(
  tau2 + sigma2 *
    (grid$number_of_subjects ^ (2*delta_n)) *
    (grid$voxel_scaled       ^ (2*delta_v))
)

grid$upper <- overall_effect + k * grid$se_model
grid$lower <- overall_effect - k * grid$se_model
grid$log_n <- log10(grid$number_of_subjects)

library(plotly)


n_seq <- 10^seq(0, 4, length.out = 60)
v_seq <- seq(
  quantile(dat$voxel_scaled, 0.05, na.rm = TRUE),
  quantile(dat$voxel_scaled, 0.95, na.rm = TRUE),
  length.out = 60
)

grid3d <- expand.grid(number_of_subjects = n_seq, voxel_scaled = v_seq)

grid3d$se_model <- sqrt(
  tau2 + sigma2 *
    (grid3d$number_of_subjects ^ (2 * delta_n)) *
    (grid3d$voxel_scaled       ^ (2 * delta_v))
)

grid3d$upper <- overall_effect + k * grid3d$se_model
grid3d$lower <- overall_effect - k * grid3d$se_model


xg <- log10(n_seq) 
yg <- v_seq        

UP <- t(matrix(grid3d$upper, nrow = length(n_seq), ncol = length(v_seq)))
LO <- t(matrix(grid3d$lower, nrow = length(n_seq), ncol = length(v_seq)))
MU <- matrix(overall_effect, nrow = length(yg), ncol = length(xg))

n_out_refit <- sum(dat$outlier)
n_tot_refit <- nrow(dat)

fig9 <- plot_ly() |>
  add_markers(
    data   = subset(dat, !outlier),
    x      = ~log10(number_of_subjects),
    y      = ~voxel_scaled,
    z      = ~es_mean,
    marker = list(size = 2, opacity = 0.25, color = "gray"),
    name   = "Inside funnel"
  ) |>
  add_markers(
    data   = subset(dat, outlier),
    x      = ~log10(number_of_subjects),
    y      = ~voxel_scaled,
    z      = ~es_mean,
    marker = list(size = 4, color = "red"),
    name   = sprintf("Outliers (n = %d / %d)", n_out_refit, n_tot_refit)
  ) |>
  add_surface(
    x          = xg, y = yg, z = MU,
    opacity    = 0.35,
    colorscale = list(c(0,1), c("royalblue","royalblue")),
    showscale  = FALSE,
    name       = "Overall effect (μ)"
  ) |>
  add_surface(
    x          = xg, y = yg, z = UP,
    opacity    = 0.20,
    colorscale = list(c(0,1), c("orange","orange")),
    showscale  = FALSE,
    name       = sprintf("Upper bound (μ + %d·SD)", k)
  ) |>
  add_surface(
    x          = xg, y = yg, z = LO,
    opacity    = 0.20,
    colorscale = list(c(0,1), c("orange","orange")),
    showscale  = FALSE,
    name       = sprintf("Lower bound (μ − %d·SD)", k)
  ) |>
  layout(
    title = list(
      text = paste0(
        "<b>Figure 9: Refitted Variance-Power Model After Outlier Exclusion</b><br>",
        "<sup>Model refit on cleaned dataset (n = ", n_tot_refit, " images). ",
        "Remaining outliers flagged in red (n = ", n_out_refit, "). ",
        "Final parameters: μ = ", round(overall_effect, 4),
        ", τ² = ", round(tau2, 4),
        ", σ² = ", round(sigma2, 4),
        ", δ_n = ", round(delta_n, 3),
        ", δ_v = ", round(delta_v, 3), ".</sup>"
      ),
      x    = 0.5,
      xref = "paper",
      font = list(size = 13)
    ),
    scene = list(
      xaxis = list(title = "log\u2081\u2080(Sample size)",       titlefont = list(size = 11)),
      yaxis = list(title = "Voxel coverage (nonzero voxels / 1e5)", titlefont = list(size = 11)),
      zaxis = list(title = "Mean effect size (ES)",              titlefont = list(size = 11))
    ),
    legend = list(
      x = 0.01, y = 0.99,
      bgcolor     = "rgba(255,255,255,0.7)",
      bordercolor = "gray",
      borderwidth = 1
    ),
    margin = list(t = 120)
  )

fig9

nv_doi <- nv[, c("id", "DOI", "paper_url")]
names(nv_doi)[1] <- "collection_id"
# SNV/copilot: merge DOI/paper_url into `dat` (the outlier-excluded final set), not the
# pre-exclusion `images` -- otherwise the merge is discarded by the `images <- dat` below.
dat$collection_id <- as.character(dat$collection_id)
nv_doi$collection_id <- as.character(nv_doi$collection_id)
dat <- merge(dat, nv_doi, by = "collection_id", all.x = TRUE)
cat("Has DOI:", sum(!is.na(dat$DOI)), "\n")

## Has DOI: 3479

cat("Has paper_url:", sum(!is.na(dat$paper_url)), "\n")

## Has paper_url: 6155

images <- dat
saveRDS(images, get_monthly_rds_path('curatedMetadataset'))
cat(sprintf("Saved curated dataset: %d images retained\n", nrow(images)))

## Saved curated dataset: 28614 images retained

6 Table Summary of the Three LMMs

To evaluate the stability of parameter estimates across data curation steps, we fit three variants of the same variance-power linear mixed-effects model:

Original: fit on the full raw dataset prior to any quality-based exclusion
Winsorized: fit after winsorizing extreme effect sizes to the 1st and 99th percentiles, used as the reference model for outlier boundary construction
Cleaned: refit on the final curated dataset after removing images flagged as outliers by the winsorized model

The model specification is identical across all three fits:

\[ \bar{\theta}_{ij} = \mu + b_i + \varepsilon_{ij}, \quad b_i \sim \mathcal{N}(0, \tau^2), \quad \varepsilon_{ij} \sim \mathcal{N}\!\left(0,\; \sigma^2 \cdot n_{ij}^{2\delta_n} \cdot v_{ij}^{2\delta_v}\right) \]

where \(\mu\) is the overall pooled effect, \(b_i\) is the collection-level random intercept, and the residual variance is modeled as a power function of sample size \(n_{ij}\) and voxel coverage \(v_{ij}\).

Table \(\ref{tab:tablesummary}\) presents the fixed effect intercept (\(\hat{\mu}\)), random effect standard deviation (\(\hat{\tau}\)), residual standard deviation (\(\hat{\sigma}\)), and the two variance-power exponents (\(\hat{\delta}_n\), \(\hat{\delta}_v\)) across all three model fits.

Across the three models, the fixed effect estimate \(\hat{\mu}\) is expected to remain stable if the overall pooled effect is robust to outliers and winsorization. Changes in \(\hat{\tau}\) across models reflect shifts in between-collection heterogeneity as extreme images are progressively removed. Similarly, reductions in \(\hat{\sigma}\) from the original to the cleaned model indicate that outlier exclusion successfully reduced residual noise. The variance-power exponents \(\hat{\delta}_n\) and \(\hat{\delta}_v\) characterize how strongly effect size variability depends on sample size and voxel coverage respectively — values close to zero indicate near-homoscedastic residuals, while larger values indicate stronger heteroscedasticity.

library(dplyr)
library(knitr)
library(kableExtra)
library(nlme)
library(tidyr)

# Extract fixed and random effects parameters from the three models
extract_model_params <- function(fit, model_name) {
  fixef_val <- as.numeric(fixef(fit))
  tau2 <- as.numeric(VarCorr(fit)[1, "Variance"])
  sigma2 <- sigma(fit)^2
  
  expParams <- intervals(fit)$varStruct[, 2]
  delta_n <- expParams[1]
  delta_v <- expParams[2]
  
  # Convert variances to standard deviations
  tau_sd <- sqrt(tau2)
  sigma_sd <- sqrt(sigma2)
  
  data.frame(
    Model = model_name,
    Parameter = c("Fixed effect (intercept)", "Random effect SD (tau)", "Residual SD (sigma)", "delta_n", "delta_v"),
    Estimate = c(fixef_val, tau_sd, sigma_sd, delta_n, delta_v),
    stringsAsFactors = FALSE
  )
}

# Extract parameters from all three models
params_combined <- bind_rows(
  extract_model_params(m_power_original, "Original"),
  extract_model_params(m_power_winsorized, "Winsorized"),
  extract_model_params(m_power, "Cleaned")
)

# Format the table
tab_display <- params_combined %>%
  pivot_wider(
    names_from = Model,
    values_from = Estimate
  ) %>%
  mutate(
    across(c(Original, Winsorized, Cleaned), function(x) {
      case_when(
        grepl("Fixed effect", Parameter) ~ sprintf("%.3f", x),
        grepl("Residual SD", Parameter) ~ sprintf("%.2e", x),
        grepl("delta", Parameter) ~ sprintf("%.3f", x),
        TRUE ~ sprintf("%.4f", x)
      )
    })
  )

kable(
  tab_display,
  format = "html",
  escape = FALSE,
  col.names = c("Parameter", "Original", "Winsorized", "Cleaned"),
  caption = "Fixed and random effects parameter estimates from the three mixed-effects models."
) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover"),
    full_width = FALSE,
    position = "center"
  )

Fixed and random effects parameter estimates from the three mixed-effects models.
Parameter	Original	Winsorized	Cleaned
Fixed effect (intercept)	0.025	0.002	0.002
Random effect SD (tau)	0.6591	0.0630	0.0361
Residual SD (sigma)	9.68e+07	7.84e+03	3.39e+05
delta_n	-0.772	-0.133	-0.171
delta_v	-3.029	-1.919	-2.608

7 Summary of the curated Neurovault dataset

Following a three-stage quality control pipeline, the final curated dataset comprises \(28614\) neuroimaging maps retained for downstream analysis.

The first stage examined tabular image properties available in the NeuroVault metadata, including sample size, voxel dimensions, and image type, removing maps that were duplicated, missing key descriptors, or outside the scope of the analysis (Section 3).

The second stage assessed spatial coverage by loading each image directly and computing gray matter coverage (\(86.4\%\) mean retained), white matter coverage (\(91.3\%\) mean retained), and outside-mask signal (\(2.6\%\) mean retained). Images with insufficient brain coverage or excessive non-brain signal were excluded at this stage (Section 4).

The third stage applied a model-based outlier detection procedure using a variance-power linear mixed-effects model fit on winsorized effect sizes. Images whose effect sizes deviated more than \(k = 3\) standard deviations from the pooled estimate — after accounting for heteroscedasticity as a function of sample size and voxel coverage — were flagged and removed (Section 5).

Table \(\ref{tab:curation_summary}\) summarizes the key image-level characteristics of the curated dataset. The distribution of sample sizes (\(n\)) reflects the typical range seen in neuroimaging meta-analyses, with substantial right skew driven by a small number of large studies. Mean effect size (\(\bar{\theta}\)) and voxel coverage metrics (gray matter percentage, white matter percentage, and outside-mask signal) characterize the spatial quality of the retained images.

Gray matter coverage (\(86.4\%\) mean) and white matter coverage (\(91.3\%\) mean) confirm that retained images provide adequate brain signal. The outside-mask percentage (\(2.6\%\) mean) confirms that non-brain signal is minimal across the curated set, consistent with successful skull-stripping and registration.

# Summary table for the final curated dataset using table1 for consistency
table1( ~ number_of_subjects + es_mean + voxel_scaled + gm_perc + wm_perc + outside_perc + grey_mean + white_mean, 
        data = images, 
        caption = "Table 6: Summary statistics for the final curated Neurovault dataset after outlier exclusions.")

Table 6: Summary statistics for the final curated Neurovault dataset after outlier exclusions.
	Overall (N=28614)
number_of_subjects
Mean (SD)	66.7 (212)
Median [Min, Max]	30.0 [3.00, 30000]
es_mean
Mean (SD)	0.0101 (0.149)
Median [Min, Max]	0.00386 [-0.727, 0.727]
voxel_scaled
Mean (SD)	231 (40.6)
Median [Min, Max]	209 [130, 343]
gm_perc
Mean (SD)	0.864 (0.0923)
Median [Min, Max]	0.814 [0.557, 1.00]
wm_perc
Mean (SD)	0.913 (0.0808)
Median [Min, Max]	0.917 [0.342, 1.00]
outside_perc
Mean (SD)	0.0257 (0.0345)
Median [Min, Max]	0.00796 [0, 0.151]
grey_mean
Mean (SD)	0.0117 (0.155)
Median [Min, Max]	0.00562 [-0.794, 0.794]
white_mean
Mean (SD)	0.00636 (0.154)
Median [Min, Max]	0.00199 [-0.823, 0.766]

8 Appendix

8.1 Collection-level variables

The following variables are available at the collection level in the NeuroVault metadata:

id, url, download_url, owner, contributors, owner_name, number_of_images, name, publication_status, DOI, preprint_DOI, authors, paper_url, journal_name, description, full_dataset_url, private, add_date, modify_date, doi_add_date, type_of_design, number_of_imaging_runs, number_of_experimental_units, length_of_runs, length_of_blocks, length_of_trials, optimization, optimization_method, subject_age_mean, subject_age_min, subject_age_max, handedness, proportion_male_subjects, inclusion_exclusion_criteria, number_of_rejected_subjects, group_comparison, group_description, scanner_make, scanner_model, field_strength, pulse_sequence, parallel_imaging, field_of_view, matrix_size, slice_thickness, skip_distance, acquisition_orientation, order_of_acquisition, repetition_time, echo_time, flip_angle, software_package, software_version, order_of_preprocessing_operations, quality_control, used_b0_unwarping, b0_unwarping_software, used_slice_timing_correction, slice_timing_correction_software, used_motion_correction, motion_correction_software, motion_correction_reference, motion_correction_metric, motion_correction_interpolation, used_motion_susceptibiity_correction, used_intersubject_registration, intersubject_registration_software, intersubject_transformation_type, nonlinear_transform_type, transform_similarity_metric, interpolation_method, object_image_type, functional_coregistered_to_structural, functional_coregistration_method, coordinate_space, target_template_image, target_resolution, used_smoothing, smoothing_type, smoothing_fwhm, resampled_voxel_size, intrasubject_model_type, intrasubject_estimation_type, intrasubject_modeling_software, hemodynamic_response_function, used_temporal_derivatives, used_dispersion_derivatives, used_motion_regressors, used_reaction_time_regressor, used_orthogonalization, orthogonalization_description, used_high_pass_filter, high_pass_filter_method, autocorrelation_model, group_model_type, group_estimation_type, group_modeling_software, group_inference_type, group_model_multilevel, group_repeated_measures, group_repeated_measures_method, nutbrain_hunger_state, nutbrain_food_viewing_conditions, nutbrain_food_choice_type, nutbrain_taste_conditions, nutbrain_odor_conditions, communities

8.2 Image-level variables

The following variables are available at the image level after preprocessing:

collection_id, add_date, collection, description, file, file_size, id, image_type, is_valid, modify_date, subject_species, target_template_image, url, name, brain_coverage, map_type, not_mni, perc_bad_voxels, data_origin, perc_voxels_outside, modality, is_thresholded, cognitive_paradigm_cogatlas, cognitive_paradigm_cogatlas_id, analysis_level, number_of_subjects, thumbnail, surface_left_file, surface_right_file, contrast_definition, binarized_group_level, binarized_map_type, unthresholded, mni, nonzero_sample, collection_title, collection_abstract, is_meta_analysis, localfile, download, voxdim, imgdim, nvox, orient, range, bounding_box, localfile_base, duplicated, dim_mm, log10_ratio_xy, log10_ratio_xz, log10_ratio_yz, ratio_x, ratio_y, ratio_z, is_proportional, range_low, range_high, range_diff, range_exclusion, localfileadj, resliced, nvox_new, voxdim_new, imgdim_new, orient_new, range_new, nvox_outside, nvox_gm, nvox_wm, q90_outside, gm_perc, wm_perc, outside_perc, gm_exclude, wm_exclude, nvox_outside_exclude, effect_size_file, brain_file, es_min, es_1qt, es_median, es_3qt, es_max, es_mean, grey_mean, grey_min, grey_max, white_mean, white_min, white_max, voxel_scaled, es_mean_winsorized, se_model, res_std, outlier, DOI, paper_url

8.3 Data Availability and Access

8.3.1 Curated Dataset

The curated NeuroVault dataset produced by this pipeline is publicly available on Zenodo and updated on the first of each month. Each monthly release contains:

Tabular metadata (curatedMetadataset_YYYY-MM.rds) — an R data frame with one row per retained image, containing all image-level metadata, coverage statistics, and effect size summaries computed during preprocessing.
Registered effect size images — brain-registered NIfTI images converted to the robust effect size index (RESI), stored as *_rigid_cost_step1_effect_size_masked_mni.nii.gz.

The dataset can be downloaded directly from Zenodo at: DOI pending publication

8.3.2 Image Naming Convention

Images in this pipeline pass through four processing stages, each appending a suffix to the original filename:

Table A1: Image file naming convention across preprocessing stages.
Stage	Column in dataset	File suffix	Description
Downloaded	localfile	.nii.gz	Original image downloaded from NeuroVault
NaN/epsilon cleaned	localfileadj	_eps.nii.gz	NaN replaced with 0; near-zero voxels (\|v\|
Rigid registration	resliced	_rigid_cost_step1.nii.gz	Registered to MNI 2mm template using FSL FLIRT (DOF=6, mutual information)
Effect size	effect_size_file	_rigid_cost_step1_effect_size.nii.gz	T/Z statistics converted to RESI effect sizes
Brain-masked effect size	brain_file	_rigid_cost_step1_effect_size_masked_mni.nii.gz	Effect size image masked to MNI brain template

All files are stored under /media/big/neurovault/{collection_id}/ on the processing server, mirroring the NeuroVault collection structure. The tar archive uploaded to Zenodo (neurovault_subset.tar) contains the registered images (resliced) for all retained collections.

8.3.3 Reproducing the Pipeline

This document is re-rendered automatically on the first of each month. To reproduce the full pipeline locally:

Clone the repository and set rds_dir and outdir to your local paths
Ensure FSL (flirt) and Convert3D (c3d) are installed and on PATH
Render the Rmd — each chunk checks for an existing monthly RDS before rerunning, so previously completed steps are skipped automatically
The final curated dataset is saved to {rds_dir}/curatedMetadataset_YYYY-MM.rds and uploaded to Zenodo

Table A2: Pipeline steps, outputs, and approximate runtimes on 24 cores.
Section	Step	Output RDS / file	Approx. runtime
1	Metadata download	neurovaultMetadata_YYYY-MM.rds	~30 min
2	Image download	localfile (NIfTI on disk)	~2 hours
3	Image-derived exclusions	imageMetadataWithDimensions_YYYY-MM.rds	~1 hour
4	Image value processing	localfileadj (*_eps.nii.gz)	~30 min
5	Registration	imageMetadataRegistered_YYYY-MM.rds	~4 hours
6	Effect size estimation	imageMetadataEffectSizes_YYYY-MM.rds	~6 hours

8.4 Session information

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: America/Chicago
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] XML_3.99-0.20      rentrez_1.2.4      gridExtra_2.3      sandwich_3.1-1    
##  [5] lmtest_0.9-40      zoo_1.8-12         papayaWidget_0.7.1 kableExtra_1.3.4  
##  [9] knitr_1.48         plotly_4.10.4      nlme_3.1-166       reactable_0.4.4   
## [13] stringr_1.6.0      RESI_1.2.4         pbapply_1.7-2      rlang_1.1.7       
## [17] scales_1.3.0       cowplot_1.1.3      gridGraphics_0.5-1 neurobase_1.32.4  
## [21] oro.nifti_0.11.4   forcats_1.0.0      tidyr_1.3.1        ComplexUpset_1.3.3
## [25] jsonlite_2.0.0     httr_1.4.8         GGally_2.2.1       ggplot2_3.5.1     
## [29] table1_1.4.3       dplyr_1.1.4        abind_1.4-8        pbj_0.1.6         
## [33] RNifti_1.7.0      
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-8         polynom_1.4-1        magrittr_2.0.4      
##  [4] multcomp_1.4-26      matrixStats_1.2.0    compiler_4.4.1      
##  [7] png_0.1-8            systemfonts_1.1.0    callr_3.7.6         
## [10] vctrs_0.7.1          rvest_1.0.5          pkgconfig_2.0.3     
## [13] fastmap_1.2.0        labeling_0.4.3       rmarkdown_2.28      
## [16] ps_1.8.1             ragg_1.2.5           purrr_1.0.2         
## [19] xfun_0.47            cachem_1.1.0         PDQutils_0.1.6      
## [22] clubSandwich_0.5.11  highr_0.11           R6_2.6.1            
## [25] bslib_0.8.0          stringi_1.8.7        RColorBrewer_1.1-3  
## [28] car_3.1-3            boot_1.3-31          jquerylib_0.1.4     
## [31] estimability_1.5.1   Rcpp_1.0.13-1        assertthat_0.2.1    
## [34] base64enc_0.1-3      R.utils_2.12.3       Matrix_1.6-3        
## [37] splines_4.4.1        tidyselect_1.2.1     rstudioapi_0.16.0   
## [40] yaml_2.3.10          codetools_0.2-20     curl_7.0.0          
## [43] processx_3.8.4       lattice_0.22-6       tibble_3.3.1        
## [46] plyr_1.8.9           withr_3.0.2          coda_0.19-4.1       
## [49] evaluate_0.24.0      orthopolynom_1.0-6.1 moments_0.14.1      
## [52] survival_3.7-0       ggstats_0.8.0        xml2_1.5.2          
## [55] pillar_1.11.1        carData_3.0-5        generics_0.1.3      
## [58] munsell_0.5.1        aod_1.3.3            xtable_1.8-4        
## [61] glue_1.8.0           emmeans_1.10.4       lazyeval_0.2.2      
## [64] tools_4.4.1          data.table_1.16.0    webshot_0.5.5       
## [67] mvtnorm_1.3-1        crosstalk_1.2.1      colorspace_2.1-1    
## [70] patchwork_1.3.0      pTFCE_0.2.2.1        mmand_1.6.3         
## [73] Formula_1.2-5        cli_3.6.5            textshaping_0.3.7   
## [76] viridisLite_0.4.2    svglite_2.1.1        gtable_0.3.6        
## [79] reactR_0.6.1         R.methodsS3_1.8.2    sass_0.4.9          
## [82] digest_0.6.37        TH.data_1.1-2        htmlwidgets_1.6.4   
## [85] farver_2.1.2         htmltools_0.5.8.1    R.oo_1.26.0         
## [88] lifecycle_1.0.5      MASS_7.3-61

Neurovault Image Preprocessing Pipeline

Haozheng Xu, Amber Hu, Simon Vandekar

2026-09-30