diff --git a/.Rbuildignore b/.Rbuildignore
index c541066..6258b8f 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -16,3 +16,4 @@
 ^Apptainer\.def$
 ^README\.Rmd$
 ^README\.html$
+^binder$
diff --git a/.lintr b/.lintr
index ed4ff18..49d559e 100644
--- a/.lintr
+++ b/.lintr
@@ -8,6 +8,7 @@ linters: linters_with_defaults(
   )
 exclusions: list(
     "inst/doc",
+    "inst/examples",
     "man",
     "vignettes",
     "data-raw",
diff --git a/DESCRIPTION b/DESCRIPTION
index 7c414aa..ca5e727 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -50,7 +50,9 @@ Suggests:
     lintr,
     stringr,
     tiledb,
-    TileDBArray
+    TileDBArray,
+    cowplot,
+    png
 Config/testthat/edition: 3
 VignetteBuilder: knitr
 URL: https://pennlinc.github.io/ModelArray, https://pennlinc.github.io/ModelArray/
diff --git a/R/utils.R b/R/utils.R
index 462ee6d..a64c474 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -73,6 +73,11 @@ flagAnalysisExistInh5 <- function(fn_h5, analysis_name) {
   if (grepl("\\.tdb/?$", filepath, ignore.case = TRUE)) {
     return("tiledb")
   }
+  # A remote URI (s3://, gcs://, azure://, ...) can only be a TileDB store;
+  # rhdf5 reads local files only. file:// is left to the local checks below.
+  if (.is_remote_uri(filepath)) {
+    return("tiledb")
+  }
   # Otherwise only claim a directory when it actually holds a store. Writers
   # delete the group they are about to replace, so a mistyped output directory
   # must not be mistaken for one.
@@ -269,6 +274,45 @@ flagAnalysisExistInh5 <- function(fn_h5, analysis_name) {
 }
 
 
+#' Is this a URI with a scheme other than file:// (for example s3://)?
+#' @noRd
+.is_remote_uri <- function(path) {
+  grepl("^[[:alpha:]][[:alnum:]+.-]*://", path) && !grepl("^file://", path, ignore.case = TRUE)
+}
+
+
+#' A TileDB VFS handle on the current default context
+#'
+#' tiledb::tiledb_get_vfs() caches a VFS built from the first context it saw,
+#' which would miss S3 settings a user applied later with tiledb_ctx().
+#' @noRd
+.tiledb_vfs <- function() {
+  tiledb::tiledb_vfs(ctx = tiledb::tiledb_get_context())
+}
+
+
+#' Does a TileDB store root exist, as a local directory or as a URI?
+#' @noRd
+.tiledb_store_exists <- function(uri) {
+  if (dir.exists(uri)) {
+    return(TRUE)
+  }
+  isTRUE(tryCatch(tiledb::tiledb_vfs_is_dir(uri, vfs = .tiledb_vfs()), error = function(e) FALSE))
+}
+
+
+#' Names of the entries directly under a TileDB URI (character(0) when absent)
+#'
+#' Uses TileDB's VFS, so it works for local paths, file:// and remote URIs
+#' alike. Plain directories are listed too (test fixtures store arrays under
+#' them), so callers still check for the array they need.
+#' @noRd
+.tiledb_list_children <- function(uri) {
+  children <- tryCatch(tiledb::tiledb_vfs_ls(uri, vfs = .tiledb_vfs()), error = function(e) character(0))
+  sort(unique(basename(sub("/+$", "", children))))
+}
+
+
 #' Convert TileDB metadata values to a character vector
 #' @noRd
 .parse_tiledb_names <- function(x) {
@@ -1126,12 +1170,28 @@ bind_cols_check_emptyTibble <- function(a, b) {
 #' into memory. The returned
 #' object has a \code{print} method that displays a formatted summary.
 #'
+#' A TileDB store can also be remote, for example in an S3 bucket. Configure
+#' the TileDB context before calling, exactly as you would before
+#' \code{\link{ModelArray}()}:
+#' \preformatted{
+#' cfg <- tiledb::tiledb_config()
+#' cfg["vfs.s3.region"] <- "us-east-1"
+#' cfg["vfs.s3.no_sign_request"] <- "true"   # public bucket, no credentials
+#' tiledb::tiledb_ctx(cfg)
+#' ModelArraySummary("s3://bucket/path/store.tdb")
+#' }
+#' The same settings can come from environment variables
+#' (\code{TILEDB_VFS_S3_REGION}, \code{TILEDB_VFS_S3_NO_SIGN_REQUEST}).
+#' Each scalar costs a few requests, so summarizing a store with many scalars
+#' over the network takes seconds rather than milliseconds.
+#'
 #' \code{h5summary()} is a backward-compatible alias for
 #' \code{ModelArraySummary()}. Despite the historical name, it supports both
 #' HDF5 files and TileDB stores.
 #'
 #' @param filepath Character. Path to an HDF5 (\code{.h5}) file or TileDB
-#'   (\code{.tdb}) store.
+#'   (\code{.tdb}) store. A TileDB store may be given as a URI
+#'   (\code{s3://...}, \code{file://...}); HDF5 files must be local paths.
 #' @param backend Character. Storage backend: \code{"auto"} (default),
 #'   \code{"hdf5"}, or \code{"tiledb"}.
 #'
@@ -1155,6 +1215,7 @@ bind_cols_check_emptyTibble <- function(a, b) {
 #' \dontrun{
 #' ModelArraySummary("path/to/data.h5")
 #' ModelArraySummary("path/to/store.tdb")
+#' ModelArraySummary("s3://bucket/path/store.tdb") # after configuring the TileDB context
 #'
 #' # Inspect before deciding which scalars to load
 #' info <- ModelArraySummary("path/to/data.h5")
@@ -1185,18 +1246,15 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
 
 #' @noRd
 .modelarray_summary_impl <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
-  if (!file.exists(filepath)) {
-    stop("File not found: ", filepath)
-  }
   backend <- .resolve_storage_backend(filepath, backend)
   if (identical(backend, "tiledb")) {
     .require_tiledb_support()
-    scalars_dir <- .modelarray_path(filepath, "scalars")
-    scalar_candidates <- if (dir.exists(scalars_dir)) {
-      basename(list.dirs(scalars_dir, recursive = FALSE, full.names = TRUE))
-    } else {
-      character(0)
+    # TileDB stores may live behind a URI (s3://, file://, ...), so existence and
+    # listing go through TileDB's VFS rather than file.exists()/list.dirs().
+    if (!.tiledb_store_exists(filepath)) {
+      stop("TileDB store not found: ", filepath)
     }
+    scalar_candidates <- .tiledb_list_children(.modelarray_path(filepath, "scalars"))
     scalar_names <- scalar_candidates[vapply(
       scalar_candidates,
       function(scalar_name) {
@@ -1240,12 +1298,7 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
       )
     }
 
-    results_dir <- .modelarray_path(filepath, "results")
-    analysis_candidates <- if (dir.exists(results_dir)) {
-      basename(list.dirs(results_dir, recursive = FALSE, full.names = TRUE))
-    } else {
-      character(0)
-    }
+    analysis_candidates <- .tiledb_list_children(.modelarray_path(filepath, "results"))
     analyses <- analysis_candidates[vapply(
       analysis_candidates,
       function(analysis_name) {
@@ -1267,6 +1320,9 @@ h5summary <- function(filepath, backend = c("auto", "hdf5", "tiledb")) {
     )
   }
 
+  if (!file.exists(filepath)) {
+    stop("File not found: ", filepath)
+  }
   listing <- rhdf5::h5ls(filepath)
   rhdf5::h5closeAll()
 
diff --git a/README.Rmd b/README.Rmd
index 403c9bc..41f9520 100644
--- a/README.Rmd
+++ b/README.Rmd
@@ -49,6 +49,7 @@ That's it.
 **It scales.** ModelArray uses [HDF5](https://www.hdfgroup.org/solutions/hdf5/) by default, and can also read TileDB stores when the optional [TileDB](https://tiledb.com/) R packages are installed. Both backends use [DelayedArray](https://bioconductor.org/packages/DelayedArray/) for lazy access.
 An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions - about 291 GB if loaded in its entirety - never enters RAM.
 ModelArray reads one element at a time, so memory usage stays flat regardless of dataset size.
+TileDB stores can also live in the cloud: `vignette("rbc-s3-tiledb")` fits GAMs to the 553 GB RBC cortical thickness store straight from a public S3 bucket, and you can run it in your browser with [Binder](https://mybinder.org/v2/gh/PennLINC/ModelArray/tiledb-example?urlpath=rstudio).
 
 **It's flexible.** Beyond linear models, ModelArray supports:
 
diff --git a/README.md b/README.md
index 30d9d41..5daa0f2 100644
--- a/README.md
+++ b/README.md
@@ -50,6 +50,10 @@ access. An ABCC dMRI dataset has ~350,000 voxels × ~26,000-sessions -
 about 291 GB if loaded in its entirety - never enters RAM. ModelArray
 reads one element at a time, so memory usage stays flat regardless of
 dataset size.
+TileDB stores can also live in the cloud: `vignette("rbc-s3-tiledb")`
+fits GAMs to the 553 GB RBC cortical thickness store straight from a
+public S3 bucket, and you can run it in your browser with
+[Binder](https://mybinder.org/v2/gh/PennLINC/ModelArray/tiledb-example?urlpath=rstudio).
 
 **It’s flexible.** Beyond linear models, ModelArray supports:
 
diff --git a/_pkgdown.yml b/_pkgdown.yml
index c1ad1ee..d9ba6f8 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -34,6 +34,8 @@ navbar:
           href: articles/modelling.html
         - text: "Exploring Data with Convenience Functions"
           href: articles/exploring-modelarray-data.html
+        - text: "Cloud TileDB: RBC Cortical Thickness from S3"
+          href: articles/rbc-s3-tiledb.html
         - text: "-------"
         - text: Resources
         - text: "Containers"
diff --git a/binder/apt.txt b/binder/apt.txt
new file mode 100644
index 0000000..7d32a8e
--- /dev/null
+++ b/binder/apt.txt
@@ -0,0 +1,6 @@
+libhdf5-dev
+libcurl4-openssl-dev
+libssl-dev
+libxml2-dev
+libgit2-dev
+bzip2
diff --git a/binder/install.R b/binder/install.R
new file mode 100644
index 0000000..5829ecc
--- /dev/null
+++ b/binder/install.R
@@ -0,0 +1,10 @@
+# Binder image for vignette("rbc-s3-tiledb"): everything the example report needs.
+# CRAN packages come from the Posit Package Manager snapshot named in runtime.txt.
+install.packages(c("remotes", "BiocManager", "tiledb", "mgcv", "ggplot2", "gratia",
+                   "gamm4", "lme4", "MASS", "rmarkdown", "knitr", "jsonlite", "cowplot", "png"))
+BiocManager::install(c("TileDBArray", "HDF5Array", "rhdf5", "DelayedArray"),
+                     update = FALSE, ask = FALSE)
+options(repos = BiocManager::repositories())
+remotes::install_github("andy1764/ComBatFamily", upgrade = "never")
+# This repository is ModelArray itself, at the branch Binder was launched from.
+remotes::install_local(".", upgrade = "never", force = TRUE)
diff --git a/binder/postBuild b/binder/postBuild
new file mode 100755
index 0000000..72659a3
--- /dev/null
+++ b/binder/postBuild
@@ -0,0 +1,13 @@
+#!/usr/bin/env bash
+# Runs once when the Binder image is built (as the notebook user, no root).
+set -euo pipefail
+
+# Connectome Workbench from conda-forge, installed into the home directory with micromamba.
+mkdir -p "$HOME/.local/bin"
+curl -Ls https://micro.mamba.pm/api/micromamba/linux-64/latest | tar -xj -C "$HOME/.local" bin/micromamba
+"$HOME/.local/bin/micromamba" create -y -p "$HOME/.wb" -c conda-forge connectome-workbench
+"$HOME/.local/bin/micromamba" clean -y --all
+
+# Put the report where the RStudio session opens, so "Knit" is the only click needed.
+mkdir -p "$HOME/rbc_s3_example"
+cp inst/examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd "$HOME/rbc_s3_example/"
diff --git a/binder/runtime.txt b/binder/runtime.txt
new file mode 100644
index 0000000..e47c783
--- /dev/null
+++ b/binder/runtime.txt
@@ -0,0 +1 @@
+r-4.4-2025-03-01
diff --git a/binder/start b/binder/start
new file mode 100755
index 0000000..be217c9
--- /dev/null
+++ b/binder/start
@@ -0,0 +1,7 @@
+#!/usr/bin/env bash
+# Runs before the server starts in every Binder session.
+export PATH="$HOME/.wb/bin:$PATH"
+# Anonymous reads from the public RBC bucket; the report also sets these itself.
+export TILEDB_VFS_S3_NO_SIGN_REQUEST=true
+export TILEDB_VFS_S3_REGION=us-east-1
+exec "$@"
diff --git a/inst/examples/rbc-s3-tiledb/README.md b/inst/examples/rbc-s3-tiledb/README.md
new file mode 100644
index 0000000..36d2dc6
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/README.md
@@ -0,0 +1,38 @@
+# RBC cortical thickness from S3: a runnable ModelArray example
+
+`rbc_s3_tiledb_example.Rmd` reads the RBC FreeSurferPost TileDB store straight from
+its public S3 bucket, fits vertex-wise GAMs of cortical thickness on age in a handful
+of Glasser parcels, and writes the results as CIFTI files you can open in Connectome
+Workbench. `vignette("rbc-s3-tiledb")` explains how to run it and what to expect.
+
+## Files
+
+| File | What it is |
+|---|---|
+| `rbc_s3_tiledb_example.Rmd` | The report. Render it with `rmarkdown::render()`; see the vignette for the exact call and parameters. |
+| `glasser360_fsLR164k_labels.rds` | One Glasser (HCP-MMP1.0) parcel key per column of the TileDB store (327,684 fsLR-164k vertices, left hemisphere then right, medial wall = 0), with the label table and each parcel's sensorimotor-association (S-A) axis rank. 26 KB. |
+| `sa_axis_glasser360_ranks.csv` | S-A axis average ranks for the 180 Glasser regions (`final.rank`: 1 = sensorimotor pole, 180 = association pole). |
+| `make_glasser_fsLR164k_labels.sh`, `make_glasser_fsLR164k_labels.R` | How the `.rds` was made: the 32k Glasser dlabel resampled to 164k with Connectome Workbench and TemplateFlow's fsLR spheres. Kept for provenance; you don't need to run them. |
+| `brain_sa_axis.png`, `brain_sa_3b.png`, `brain_sa_46.png` | Lateral views of the left inflated cortex coloured by the S-A axis in the colour scheme of Sydnor et al. 2021 Figure 2A (gold = sensorimotor, purple = association). `brain_sa_axis` outlines and labels both areas 3b and 46; the other two outline one parcel each. The figures use them as locators. |
+| `make_brain_sa_axis.R` | How the three brain PNGs were made (filled-mesh render of the 164k inflated surface). Provenance; you don't need to run it. |
+
+## Sources and attribution
+
+- **Parcellation:** Glasser MF et al. (2016). A multi-modal parcellation of human
+  cerebral cortex. *Nature* 536, 171–178. The 32k dlabel is the copy distributed with
+  [XCP-D](https://github.com/PennLINC/xcp_d) (`atlas-Glasser_space-fsLR_den-32k_dseg.dlabel.nii`);
+  the HCP-MMP1.0 parcellation is released under the HCP open-access data use terms.
+- **S-A axis ranks:** Sydnor VJ et al. (2021). Neurodevelopment of the association
+  cortices: patterns, mechanisms, and implications for psychopathology. *Neuron* 109,
+  2820–2846. File from [PennLINC/S-A_ArchetypalAxis](https://github.com/PennLINC/S-A_ArchetypalAxis)
+  (`Glasser360_MMP/Sensorimotor_Association_Axis_AverageRanks.csv`). Please cite
+  Sydnor et al. 2021 if you use the ranks. The S-A axis brain (`brain_sa_axis.png`)
+  reuses that paper's Figure 2A colour scheme (gold → white → purple, diverging at
+  rank 90).
+- **Surfaces for resampling and display:** [TemplateFlow](https://www.templateflow.org/)
+  `tpl-fsLR` (spheres, vertex-area maps and 32k display surfaces), downloaded by the
+  scripts when needed.
+- **Data:** Reproducible Brain Charts (RBC) FreeSurferPost release, published as a
+  TileDB store at `s3://fcp-indi/data/Projects/RBC/RBC_FreeSurfer-Post_ModelArray/`.
+  See the RBC documentation at <https://reprobrainchart.github.io/> for the data
+  release and its terms.
diff --git a/inst/examples/rbc-s3-tiledb/brain_sa_3b.png b/inst/examples/rbc-s3-tiledb/brain_sa_3b.png
new file mode 100644
index 0000000..572127b
Binary files /dev/null and b/inst/examples/rbc-s3-tiledb/brain_sa_3b.png differ
diff --git a/inst/examples/rbc-s3-tiledb/brain_sa_46.png b/inst/examples/rbc-s3-tiledb/brain_sa_46.png
new file mode 100644
index 0000000..0596436
Binary files /dev/null and b/inst/examples/rbc-s3-tiledb/brain_sa_46.png differ
diff --git a/inst/examples/rbc-s3-tiledb/brain_sa_axis.png b/inst/examples/rbc-s3-tiledb/brain_sa_axis.png
new file mode 100644
index 0000000..e1db676
Binary files /dev/null and b/inst/examples/rbc-s3-tiledb/brain_sa_axis.png differ
diff --git a/inst/examples/rbc-s3-tiledb/glasser360_fsLR164k_labels.rds b/inst/examples/rbc-s3-tiledb/glasser360_fsLR164k_labels.rds
new file mode 100644
index 0000000..b4ae788
Binary files /dev/null and b/inst/examples/rbc-s3-tiledb/glasser360_fsLR164k_labels.rds differ
diff --git a/inst/examples/rbc-s3-tiledb/make_brain_sa_axis.R b/inst/examples/rbc-s3-tiledb/make_brain_sa_axis.R
new file mode 100644
index 0000000..ea8d88b
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/make_brain_sa_axis.R
@@ -0,0 +1,95 @@
+#!/usr/bin/env Rscript
+# Build the S-A axis brain images that ship with the example: lateral views of the left
+# inflated fsLR-164k surface coloured by the Glasser sensorimotor-association (S-A) axis,
+# in the colour scheme of Sydnor et al. 2021 (Neuron), Figure 2A -- gold (sensorimotor)
+# through white to purple (association), diverging at rank 90.
+#
+#   make_brain_sa_axis.R <L inflated 164k surf.gii> <glasser360_fsLR164k_labels.rds> <out dir> [wb_command]
+#
+# Writes three PNGs into <out dir>:
+#   brain_sa_axis.png  areas 3b and 46 outlined and labelled (the main locator)
+#   brain_sa_3b.png    area 3b outlined only (for the 3b panel of the vertex figure)
+#   brain_sa_46.png    area 46 outlined only (for the 46 panel of the vertex figure)
+#
+# The 164k inflated surface is the one the example builds from TemplateFlow's 32k surfaces
+# (see the report's cifti display step, or make_glasser_fsLR164k_labels.sh). This is
+# provenance for files that ship with the example; users do not run it. Needs wb_command
+# (to convert the GIFTI to ASCII), xml2, ggplot2 and scales.
+#
+# S-A axis ranks: Sydnor VJ et al. 2021, Neuron (PennLINC/S-A_ArchetypalAxis). Cite Sydnor
+# et al. 2021 if you use them. The surface is rendered as shaded, back-face-culled filled
+# triangles; parcel outlines are the mesh edges between in- and out-of-parcel triangles.
+suppressMessages({library(xml2); library(ggplot2); library(scales)})
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) < 3) stop("usage: make_brain_sa_axis.R L_inflated.surf.gii labels.rds out_dir [wb_command]")
+surf_gii <- args[1]; labels_rds <- args[2]; out_dir <- args[3]
+wb <- if (length(args) >= 4) args[4] else "wb_command"
+dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)
+
+ascii <- tempfile(fileext = ".surf.gii")
+stopifnot(system2(wb, c("-gifti-convert", "ASCII", surf_gii, ascii)) == 0)
+doc <- read_xml(ascii)
+V <- matrix(as.numeric(scan(text = xml_text(xml_find_first(doc, ".//DataArray[@Intent='NIFTI_INTENT_POINTSET']//Data")), quiet = TRUE)), ncol = 3, byrow = TRUE)
+Tt <- matrix(as.integer(scan(text = xml_text(xml_find_first(doc, ".//DataArray[@Intent='NIFTI_INTENT_TRIANGLE']//Data")), quiet = TRUE)), ncol = 3, byrow = TRUE) + 1L
+labels <- readRDS(labels_rds)
+keyL <- labels$key[seq_len(nrow(V))]
+stopifnot(nrow(V) == 163842L)
+
+SA_COLS <- c("goldenrod1", "white", "#6f1282"); SA_MID <- 90; SA_RANGE <- c(1, 180)   # Sydnor 2021 Fig 2A
+rank_col <- function(r) gradient_n_pal(SA_COLS)(rescale_mid(r, to = c(0, 1), from = SA_RANGE, mid = SA_MID))
+WALL <- "grey92"
+rankv <- labels$table$sa_rank[match(keyL, labels$table$key)]
+sacol <- ifelse(is.na(rankv), WALL, rank_col(rankv))
+
+blend <- function(hex, f) { m <- grDevices::col2rgb(hex) / 255; m <- sweep(m, 2, f, `*`); m[m > 1] <- 1; m[m < 0] <- 0; grDevices::rgb(m[1,], m[2,], m[3,]) }
+# left lateral view: anterior (+Y) to the left, superior (+Z) up; keep lateral-facing triangles
+xi <- -V[, 2]; yi <- V[, 3]
+a <- V[Tt[,1],]; b <- V[Tt[,2],]; c <- V[Tt[,3],]
+n <- cbind((b[,2]-a[,2])*(c[,3]-a[,3])-(b[,3]-a[,3])*(c[,2]-a[,2]),
+           (b[,3]-a[,3])*(c[,1]-a[,1])-(b[,1]-a[,1])*(c[,3]-a[,3]),
+           (b[,1]-a[,1])*(c[,2]-a[,2])-(b[,2]-a[,2])*(c[,1]-a[,1]))
+nn <- sqrt(rowSums(n^2)); nn[nn == 0] <- 1; n <- n / nn
+ff <- n[, 1] < 0; Tf <- Tt[ff, , drop = FALSE]; nf <- n[ff, , drop = FALSE]
+light <- c(-0.5, 0.25, 0.85); light <- light / sqrt(sum(light^2))
+shade <- pmax(0.6, 0.78 + 0.22 * as.vector(nf %*% light))
+tcol <- blend(sacol[Tf[, 1]], shade)
+depth <- (V[Tf[,1],1] + V[Tf[,2],1] + V[Tf[,3],1]) / 3
+ord <- order(depth, decreasing = TRUE); Tf <- Tf[ord, ]; tcol <- tcol[ord]; m <- nrow(Tf)
+poly <- data.frame(x = c(xi[Tf[,1]], xi[Tf[,2]], xi[Tf[,3]]), y = c(yi[Tf[,1]], yi[Tf[,2]], yi[Tf[,3]]),
+                   grp = rep(seq_len(m), 3), fill = rep(tcol, 3))
+poly <- poly[order(poly$grp), ]
+base_plot <- ggplot(poly, aes(x, y, group = grp, fill = fill)) + geom_polygon(colour = NA) +
+  scale_fill_identity() + coord_fixed() + theme_void()
+
+outline_df <- function(nm) {
+  inpar <- keyL == labels$table$key[labels$table$name == nm]
+  tri_in <- rowSums(matrix(inpar[Tf], ncol = 3)) >= 2L
+  e <- rbind(Tf[,c(1,2)], Tf[,c(2,3)], Tf[,c(1,3)]); e <- cbind(pmin(e[,1],e[,2]), pmax(e[,1],e[,2]))
+  k <- paste(e[,1], e[,2])
+  keep <- ave(rep(1L, nrow(e)), k, FUN = sum) == 2L & ave(as.integer(rep(tri_in, 3)), k, FUN = sum) == 1L & !duplicated(k)
+  e <- e[keep, , drop = FALSE]
+  data.frame(x = xi[e[,1]], y = yi[e[,1]], xend = xi[e[,2]], yend = yi[e[,2]])
+}
+centroid <- function(nm) { v <- which(keyL == labels$table$key[labels$table$name == nm]); data.frame(x = mean(xi[v]), y = mean(yi[v])) }
+
+# region name -> (glasser parcel, label text, label nudge)
+PARCELS <- list("3b" = list(roi = "L_3b_ROI", dx = -14, dy = 16),
+                "46" = list(roi = "L_46_ROI", dx = -2,  dy = 18))
+render <- function(regions, label) {
+  g <- base_plot
+  for (rg in regions) g <- g + geom_segment(data = outline_df(PARCELS[[rg]]$roi),
+    aes(x, y, xend = xend, yend = yend), inherit.aes = FALSE, colour = "black", linewidth = 0.45)
+  if (label) {
+    labs <- do.call(rbind, lapply(regions, function(rg) cbind(centroid(PARCELS[[rg]]$roi), lab = rg,
+      dx = PARCELS[[rg]]$dx, dy = PARCELS[[rg]]$dy)))
+    g <- g + geom_text(data = labs, aes(x + dx, y + dy, label = lab), inherit.aes = FALSE, colour = "black", fontface = 2, size = 5.4)
+  }
+  g
+}
+save_png <- function(g, file) {
+  ggsave(file.path(out_dir, file), g, width = 3.4, height = 2.8, dpi = 200, bg = "transparent")
+  message(sprintf("wrote %s (%.0f KB)", file, file.size(file.path(out_dir, file)) / 1024))
+}
+save_png(render(c("3b", "46"), label = TRUE),  "brain_sa_axis.png")
+save_png(render("3b", label = FALSE),          "brain_sa_3b.png")
+save_png(render("46", label = FALSE),          "brain_sa_46.png")
diff --git a/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.R b/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.R
new file mode 100644
index 0000000..afad559
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.R
@@ -0,0 +1,66 @@
+#!/usr/bin/env Rscript
+# Second half of make_glasser_fsLR164k_labels.sh: turn the text export of the 164k
+# Glasser dlabel into glasser360_fsLR164k_labels.rds, with checks.
+#
+#   Rscript make_glasser_fsLR164k_labels.R <keys.txt> <labeltable.txt> <sa_ranks.csv> <out.rds> <input dlabel> <wb version>
+#
+# The .rds holds a list:
+#   key        integer vector, one entry per TileDB column (327,684): 0 = medial wall / unlabelled
+#   table      data.frame key, name, hemi, region, sa_rank: one row per parcel (360)
+#   provenance list: inputs, method, date, software
+args <- commandArgs(trailingOnly = TRUE)
+if (length(args) < 5) stop("usage: make_glasser_fsLR164k_labels.R keys.txt labeltable.txt sa_ranks.csv out.rds input.dlabel.nii [wb_version]")
+keys_txt <- args[1]; table_txt <- args[2]; sa_csv <- args[3]; out_rds <- args[4]; input_dlabel <- args[5]
+wb_version <- if (length(args) >= 6) args[6] else "unknown"
+
+key <- as.integer(scan(keys_txt, what = numeric(), quiet = TRUE))
+n_hemi <- 163842L
+stopifnot(length(key) == 2L * n_hemi)
+
+# -cifti-label-export-table writes pairs of lines: the name, then "key R G B A"
+lines <- readLines(table_txt)
+lines <- lines[nzchar(trimws(lines))]
+stopifnot(length(lines) %% 2 == 0)
+names_ <- lines[seq(1, length(lines), by = 2)]
+keys_ <- as.integer(vapply(strsplit(trimws(lines[seq(2, length(lines), by = 2)]), "[[:space:]]+"), `[`, "", 1))
+table <- data.frame(key = keys_, name = names_, stringsAsFactors = FALSE)
+table <- table[table$key != 0L, ]                       # drop the "???" background label
+table$hemi <- substr(table$name, 1, 1)
+table$region <- sub("^[LR]_(.*)_ROI$", "\\1", table$name)
+stopifnot(nrow(table) == 360L, all(table$hemi %in% c("L", "R")), length(unique(table$region)) == 180L)
+
+# Sensorimotor-association axis rank per region (Sydnor et al. 2021), same for both hemispheres
+sa <- read.csv(sa_csv, stringsAsFactors = FALSE)
+table$sa_rank <- sa$final.rank[match(table$region, sa$region)]
+missing_rank <- table$region[is.na(table$sa_rank)]
+if (length(missing_rank) > 0) stop("No S-A rank for: ", paste(unique(missing_rank), collapse = ", "))
+
+# Alignment with the TileDB column order: left vertices first, then right
+left_keys <- table$key[table$hemi == "L"]
+stopifnot(
+  all(key[seq_len(n_hemi)] %in% c(0L, left_keys)),
+  all(key[(n_hemi + 1):(2 * n_hemi)] %in% c(0L, table$key[table$hemi == "R"])),
+  all(table$key %in% key)
+)
+table <- table[order(table$key), ]
+rownames(table) <- NULL
+
+counts <- table(factor(key, levels = c(0L, table$key)))
+message(sprintf("medial wall / unlabelled: %d of %d vertices; parcel sizes %d-%d vertices (median %d)",
+                counts[["0"]], length(key), min(counts[-1]), max(counts[-1]), as.integer(median(counts[-1]))))
+
+out <- list(
+  key = key,
+  table = table,
+  provenance = list(
+    description = "Glasser et al. 2016 HCP-MMP1.0 parcellation on fsLR 164k, one key per grayordinate of the RBC FreeSurferPost dscalars (CORTEX_LEFT 0-163841, then CORTEX_RIGHT; medial wall = 0)",
+    input_dlabel = basename(input_dlabel),
+    input_md5 = unname(tools::md5sum(input_dlabel)),
+    method = "wb_command -cifti-separate; -label-resample ADAP_BARY_AREA with TemplateFlow tpl-fsLR 32k/164k spheres and desc-vaavg_midthickness vertex areas; -cifti-create-label; -cifti-convert -to-text",
+    sa_axis = "Sensorimotor_Association_Axis_AverageRanks.csv, PennLINC/S-A_ArchetypalAxis (Sydnor et al. 2021, Neuron), column final.rank",
+    workbench = wb_version,
+    date = format(Sys.Date())
+  )
+)
+saveRDS(out, out_rds, compress = "xz")
+message(sprintf("wrote %s (%.0f KB)", out_rds, file.size(out_rds) / 1024))
diff --git a/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.sh b/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.sh
new file mode 100755
index 0000000..684f1ab
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/make_glasser_fsLR164k_labels.sh
@@ -0,0 +1,57 @@
+#!/usr/bin/env bash
+# Build glasser360_fsLR164k_labels.rds: one Glasser (HCP-MMP1.0) parcel key for each of
+# the 327,684 columns of the RBC FreeSurferPost TileDB store (fsLR 164k, all vertices,
+# left hemisphere then right, medial wall = 0).
+#
+# This is provenance for a file that ships with the example; users of the example do
+# not run it. It needs Connectome Workbench (wb_command), curl and Rscript.
+#
+#   make_glasser_fsLR164k_labels.sh <atlas-Glasser_space-fsLR_den-32k_dseg.dlabel.nii> <workdir>
+#
+# The 32k Glasser dlabel is the one distributed with XCP-D
+# (xcp_d/data/atlases/atlas-Glasser/). Spheres and vertex-area maps come from
+# TemplateFlow's tpl-fsLR; the grayordinate layout is checked against one release
+# dscalar, so the text export is in the same order as the TileDB columns.
+set -euo pipefail
+IN=${1:?usage: make_glasser_fsLR164k_labels.sh <glasser 32k dlabel> <workdir>}
+W=${2:?usage: make_glasser_fsLR164k_labels.sh <glasser 32k dlabel> <workdir>}
+HERE=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)
+WB=${WB_COMMAND:-wb_command}
+TF=https://templateflow.s3.amazonaws.com/tpl-fsLR
+TEMPLATE_URL=https://fcp-indi.s3.amazonaws.com/data/Projects/RBC/CCNP_FreeSurfer-Post/SHA256E-s3385504--99c174365d91c22508a0de725b40e19f8e6b572c6efd88dfcd681e392e187267.nii
+TEMPLATE_SHA=99c174365d91c22508a0de725b40e19f8e6b572c6efd88dfcd681e392e187267
+
+mkdir -p "$W"
+for h in L R; do
+  for d in 32k 164k; do
+    for f in "tpl-fsLR_hemi-${h}_den-${d}_sphere.surf.gii" "tpl-fsLR_hemi-${h}_den-${d}_desc-vaavg_midthickness.shape.gii"; do
+      [ -s "$W/$f" ] || curl -sSfL -o "$W/$f" "$TF/$f"
+    done
+  done
+done
+[ -s "$W/template_fsLR164k.dscalar.nii" ] || curl -sSfL -o "$W/template_fsLR164k.dscalar.nii" "$TEMPLATE_URL"
+echo "$TEMPLATE_SHA  $W/template_fsLR164k.dscalar.nii" | shasum -a 256 -c - >/dev/null
+
+# 32k dlabel -> one label.gii per hemisphere -> 164k with area-corrected barycentric resampling
+"$WB" -cifti-separate "$IN" COLUMN -label CORTEX_LEFT "$W/L_32k.label.gii" -label CORTEX_RIGHT "$W/R_32k.label.gii"
+for h in L R; do
+  "$WB" -label-resample "$W/${h}_32k.label.gii" \
+    "$W/tpl-fsLR_hemi-${h}_den-32k_sphere.surf.gii" "$W/tpl-fsLR_hemi-${h}_den-164k_sphere.surf.gii" \
+    ADAP_BARY_AREA "$W/${h}_164k.label.gii" \
+    -area-metrics "$W/tpl-fsLR_hemi-${h}_den-32k_desc-vaavg_midthickness.shape.gii" \
+                  "$W/tpl-fsLR_hemi-${h}_den-164k_desc-vaavg_midthickness.shape.gii"
+done
+"$WB" -cifti-create-label "$W/glasser_fsLR164k.dlabel.nii" \
+  -left-label "$W/L_164k.label.gii" -right-label "$W/R_164k.label.gii"
+
+# Both files must list all 163,842 vertices per hemisphere, left before right
+for f in glasser_fsLR164k.dlabel.nii template_fsLR164k.dscalar.nii; do
+  echo "$f: $("$WB" -file-information "$W/$f" | grep -E 'Number of Rows|out of' | tr -s ' ' | paste -sd' ' -)"
+done
+
+"$WB" -cifti-convert -to-text "$W/glasser_fsLR164k.dlabel.nii" "$W/glasser_fsLR164k_keys.txt"
+"$WB" -cifti-label-export-table "$W/glasser_fsLR164k.dlabel.nii" 1 "$W/glasser_fsLR164k_labeltable.txt"
+Rscript "$HERE/make_glasser_fsLR164k_labels.R" \
+  "$W/glasser_fsLR164k_keys.txt" "$W/glasser_fsLR164k_labeltable.txt" \
+  "$HERE/sa_axis_glasser360_ranks.csv" "$HERE/glasser360_fsLR164k_labels.rds" \
+  "$IN" "$("$WB" -version | head -1)"
diff --git a/inst/examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd b/inst/examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd
new file mode 100644
index 0000000..1169afa
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd
@@ -0,0 +1,986 @@
+---
+title: "ModelArray on a TileDB store in S3: cortical thickness across the sensorimotor-association axis"
+date: "`r format(Sys.time(), '%Y-%m-%d %H:%M %Z')`"
+output:
+  html_document:
+    toc: true
+    toc_depth: 2
+    df_print: paged
+params:
+  store_uri: "s3://fcp-indi/data/Projects/RBC/RBC_FreeSurfer-Post_ModelArray/rbc_freesurferpost.tdb"
+  cohort_url: "https://fcp-indi.s3.amazonaws.com/data/Projects/RBC/RBC_FreeSurfer-Post_ModelArray/rbc_freesurferpost_cohort.csv"
+  s3_region: "us-east-1"
+  n_cores: 4
+  out_dir: "~/rbc_s3_example"
+  cache_dir: NULL
+  quick: TRUE
+  hemispheres: NULL
+  parcels: NULL
+  max_tiles_per_parcel: NULL
+  max_vertices_per_parcel: NULL
+  allocation_mb: 512
+  max_age: 23
+  harmonize: TRUE
+  write_cifti: TRUE
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = "#>")
+report_started <- Sys.time()
+
+# ---- presets -----------------------------------------------------------------
+# Binder gives one core and about 2 GB, so default to the smallest run there.
+on_binder <- nzchar(Sys.getenv("BINDER_REPO_URL"))
+p <- params
+p$out_dir <- path.expand(p$out_dir)
+if (is.null(p$cache_dir)) p$cache_dir <- file.path(p$out_dir, "cache")
+p$cache_dir <- path.expand(p$cache_dir)
+if (on_binder) {
+  p$quick <- TRUE
+  p$n_cores <- 1L
+  p$allocation_mb <- min(p$allocation_mb, 256)
+}
+if (is.null(p$hemispheres)) p$hemispheres <- if (p$quick) "L" else "LR"
+if (is.null(p$parcels)) p$parcels <- c("3b", "46")
+if (is.null(p$max_tiles_per_parcel)) p$max_tiles_per_parcel <- if (on_binder) 1L else if (p$quick) 2L else 6L
+if (is.null(p$max_vertices_per_parcel)) p$max_vertices_per_parcel <- if (on_binder) 150L else if (p$quick) 300L else 500L
+dir.create(p$out_dir, recursive = TRUE, showWarnings = FALSE)
+dir.create(p$cache_dir, recursive = TRUE, showWarnings = FALSE)
+
+# ---- packages ----------------------------------------------------------------
+needed <- c("ModelArray", "tiledb", "TileDBArray", "jsonlite", "rhdf5", "mgcv", "ggplot2", "MASS")
+missing_pkgs <- needed[!vapply(needed, requireNamespace, logical(1), quietly = TRUE)]
+if (length(missing_pkgs) > 0) {
+  stop("Please install these packages first: ", paste(missing_pkgs, collapse = ", "),
+       ". See vignette('rbc-s3-tiledb') for the install commands.", call. = FALSE)
+}
+have_combat <- requireNamespace("ComBatFamily", quietly = TRUE)
+have_lme4 <- requireNamespace("lme4", quietly = TRUE)
+have_gamm4 <- requireNamespace("gamm4", quietly = TRUE)
+have_gratia <- requireNamespace("gratia", quietly = TRUE)
+have_cowplot <- requireNamespace("cowplot", quietly = TRUE)   # optional: brain-locator insets on the figures
+harmonize <- isTRUE(p$harmonize) && have_combat
+suppressPackageStartupMessages({
+  library(ModelArray)
+  library(ggplot2)
+})
+if (!exists(".require_tiledb_support", envir = asNamespace("ModelArray"))) {
+  stop("This ModelArray predates TileDB support. Install the current GitHub main (see the vignette).", call. = FALSE)
+}
+
+# ---- Connectome Workbench ------------------------------------------------------
+wb <- Sys.which("wb_command")
+if (!nzchar(wb)) {
+  for (candidate in c("/Applications/wb_view.app/Contents/usr/bin/wb_command",
+                      "/Applications/ConnectomeWorkbench.app/Contents/usr/bin/wb_command",
+                      path.expand("~/.wb/bin/wb_command"))) {
+    if (file.exists(candidate)) { wb <- candidate; break }
+  }
+}
+have_wb <- nzchar(wb) && isTRUE(p$write_cifti)
+wb_run <- function(...) {
+  out <- system2(wb, c(...), stdout = TRUE, stderr = TRUE)
+  if (!is.null(attr(out, "status")) && attr(out, "status") != 0) {
+    stop("wb_command failed:\n", paste(out, collapse = "\n"), call. = FALSE)
+  }
+  invisible(out)
+}
+
+# ---- TileDB read settings and helpers ----------------------------------------
+# tiledb's R package reads a query in passes sized by this preference (10 MiB by
+# default). Each pass re-fetches the tiles it touches, so a read of a few tiles can
+# download them six times over. One pass per parcel read needs about 67 MB per tile.
+tiledb::set_allocation_size_preference(as.integer(p$allocation_mb) * 1024^2)
+TILE_ROWS <- 964L                      # vertices per tile in this store
+N_SESSIONS <- 8696L
+TILE_MB <- TILE_ROWS * N_SESSIONS * 4 / 1024^2
+tile_of <- function(i) (as.integer(i) - 1L) %/% TILE_ROWS   # 1-based element -> 0-based tile
+
+runtime <- data.frame(step = character(), secs = numeric(), tiles = numeric(),
+                      mb_uncompressed = numeric(), mb_wire = numeric(), stringsAsFactors = FALSE)
+stat_value <- function(raw, key) {
+  m <- regmatches(raw, regexpr(sprintf('"%s": *[0-9.]+', key), raw))
+  if (length(m) == 0) return(NA_real_)
+  as.numeric(sub(".*: *", "", m))
+}
+# Run expr with TileDB statistics on, and record time, tiles and bytes for the table.
+tdb_measure <- function(step, expr) {
+  tiledb::tiledb_stats_reset()
+  tiledb::tiledb_stats_enable()
+  t0 <- Sys.time()
+  value <- force(expr)
+  secs <- as.numeric(Sys.time() - t0, units = "secs")
+  raw <- paste(suppressMessages(suppressWarnings(tiledb::tiledb_stats_raw_get())), collapse = "")
+  tiledb::tiledb_stats_disable()
+  runtime[nrow(runtime) + 1, ] <<- list(step, secs,
+    stat_value(raw, "Context.Query.Reader.num_tiles_read"),
+    stat_value(raw, "Context.Query.Reader.read_unfiltered_byte_num") / 1024^2,
+    stat_value(raw, "Context.VFS.read_byte_num") / 1024^2)
+  value
+}
+# Time a local computation (no S3 traffic) into the same table.
+timed <- function(step, expr) {
+  t0 <- Sys.time()
+  value <- force(expr)
+  runtime[nrow(runtime) + 1, ] <<- list(step, as.numeric(Sys.time() - t0, units = "secs"), NA_real_, NA_real_, NA_real_)
+  value
+}
+last_runtime <- function() runtime[nrow(runtime), ]
+cached_download <- function(url, name = basename(url)) {
+  if (file.exists(url)) return(url)                # a local file was given instead of a URL
+  dest <- file.path(p$cache_dir, name)
+  if (!file.exists(dest) || file.size(dest) == 0) {
+    utils::download.file(url, dest, mode = "wb", quiet = TRUE)
+  }
+  dest
+}
+fmt_secs <- function(s) if (s < 90) sprintf("%.1f s", s) else sprintf("%.1f min", s / 60)
+```
+
+This report fits models to the [Reproducible Brain Charts](https://reprobrainchart.github.io/)
+FreeSurferPost cortical thickness data **without downloading the dataset**. The data
+live in a 553 GB TileDB store on a public S3 bucket; ModelArray reads only the pieces
+it needs. The question is one with a known answer, so you can judge the result:
+
+> The cortex thins through childhood and adolescence, but not everywhere at the same
+> time. Primary sensory and motor cortex matures first: its thinning is mostly over by
+> mid-adolescence. Association cortex matures last: it keeps thinning through the teens
+> and into the twenties. This example contrasts two parcels at opposite ends of the
+> sensorimotor–association (S-A) axis of Sydnor et al. (2021): area **3b**, primary
+> somatosensory cortex (S-A rank 7), and area **46**, dorsolateral prefrontal
+> association cortex (S-A rank 174). 3b should do most of its thinning early and flatten
+> by mid-adolescence; 46 should still be thinning in the late teens and twenties.
+
+Settings for this run: `quick = `r p$quick``, hemispheres `r p$hemispheres`, parcels
+`r paste(p$parcels, collapse = ", ")`, at most `r p$max_tiles_per_parcel` tile(s) and
+`r p$max_vertices_per_parcel` vertices per parcel, ages up to `r p$max_age` for the
+vertex-wise models, `r p$n_cores` core(s), site harmonization
+`r if (harmonize) "on (CovBat)" else "off (site as a covariate)"`.
+`r if (on_binder) "Running on Binder, so the smallest settings are used." else ""`
+
+## 1. Point TileDB at the bucket
+
+ModelArray talks to TileDB through the package's default *context*. Two settings make
+the public bucket readable: the region, and "do not sign requests" (no AWS account
+needed). Set them once, before `ModelArray()`.
+
+```{r tiledb-context}
+cfg <- tiledb::tiledb_config()
+cfg["vfs.s3.region"] <- p$s3_region
+cfg["vfs.s3.no_sign_request"] <- "true"
+tiledb::tiledb_ctx(cfg)
+
+reachable <- tryCatch(tiledb::tiledb_vfs_is_dir(p$store_uri, tiledb::tiledb_vfs()), error = function(e) e)
+if (!isTRUE(reachable)) {
+  cat("Could not reach", p$store_uri, "\n",
+      if (inherits(reachable, "error")) conditionMessage(reachable) else "", "\n",
+      "Check your network and proxy settings, or try again later.\n")
+  knitr::knit_exit()
+}
+```
+
+The same two settings can come from environment variables instead
+(`TILEDB_VFS_S3_REGION`, `TILEDB_VFS_S3_NO_SIGN_REQUEST`), which is what the Binder
+image does.
+
+## 2. What is in the store
+
+`ModelArraySummary()` lists the scalars without loading any data. Every scalar is a
+dense matrix of 327,684 fsLR-164k vertices by 8,696 sessions, stored in tiles of
+`r TILE_ROWS` vertices × all sessions (about `r round(TILE_MB)` MB each before
+compression). Reading one vertex means reading its whole tile; that is the unit all
+the costs below are counted in.
+
+```{r summary}
+info <- tdb_measure("ModelArraySummary (70 scalars)", ModelArraySummary(p$store_uri))
+cat(sprintf("%d scalars, each %d vertices x %d sessions; %s\n",
+            nrow(info$scalars), info$scalars$nElements[1], info$scalars$nInputFiles[1],
+            fmt_secs(last_runtime()$secs)))
+print(info$scalars$name)
+```
+
+## 3. The cohort
+
+The cohort CSV sits next to the store. It is long format (one row per session per
+scalar), so we keep the `thickness` rows: one per session, with age, sex and scanning
+site. Sessions without a sex or site (some HBN sessions) keep an empty value; models
+and harmonization skip them where they need the column.
+
+```{r cohort}
+cohort_file <- cached_download(p$cohort_url, "rbc_freesurferpost_cohort.csv")
+cohort <- read.csv(cohort_file, na.strings = c("", "NA"), stringsAsFactors = FALSE)
+phen <- cohort[cohort$scalar_name == "thickness", ]
+phen$subject <- paste(phen$dataset, phen$subject_id)
+have_covariates <- all(c("sex", "site") %in% names(phen))
+if (have_covariates) {
+  phen$sex <- factor(phen$sex)
+  phen$site <- factor(phen$site)
+} else {
+  cat("This cohort file has no sex/site columns yet; using the study (dataset) as the batch variable.\n")
+  phen$site <- factor(phen$dataset)
+  phen$sex <- factor(rep("unknown", nrow(phen)))
+}
+use_sex <- nlevels(phen$sex) > 1
+sex_term <- if (use_sex) " + sex" else ""
+site_term <- if (nlevels(phen$site) > 1) " + site" else ""
+harmonize_covariates <- c("Age", if (use_sex) "sex")
+
+by_study <- split(phen, phen$dataset)
+sessions <- data.frame(
+  dataset = names(by_study),
+  sessions = vapply(by_study, nrow, integer(1)),
+  subjects = vapply(by_study, function(d) length(unique(d$subject)), integer(1)),
+  age_min = vapply(by_study, function(d) round(min(d$Age), 1), numeric(1)),
+  age_max = vapply(by_study, function(d) round(max(d$Age), 1), numeric(1)),
+  no_sex = vapply(by_study, function(d) sum(is.na(d$sex)), integer(1)),
+  no_site = vapply(by_study, function(d) sum(is.na(d$site)), integer(1)),
+  row.names = NULL
+)
+sessions
+cat(sprintf("%d sessions, %d subjects, %d sites\n", nrow(phen), length(unique(phen$subject)), nlevels(phen$site)))
+```
+
+## 4. Open the store
+
+Only the scalar we name gets opened; nothing is read yet. The `curv` scalar would be a
+natural second covariate, but ModelArray attaches a second scalar only when it lists
+the same source files as the phenotypes, and in this store each scalar has its own
+file per session. So the models here use `thickness` alone.
+
+```{r open}
+ma <- tdb_measure("ModelArray() open", ModelArray(p$store_uri, scalar_types = "thickness"))
+ma
+src <- sources(ma)[["thickness"]]
+stopifnot(setequal(src, phen$source_file))
+phen <- phen[match(src, phen$source_file), ]           # same order as the store's columns
+cat(sprintf("Opened in %s\n", fmt_secs(last_runtime()$secs)))
+```
+
+## 5. Parcels and tiles
+
+A Glasser (HCP-MMP1.0) label for each of the 327,684 columns ships with the example,
+together with each parcel's S-A rank. Before reading anything we can count how many
+tiles a parcel touches. fsLR numbers its vertices in nested rings, so a parcel's
+vertices are spread over many tiles; the report keeps the vertices in each parcel's
+most populated tiles, which caps the download while keeping hundreds of vertices.
+
+```{r atlas}
+labels <- readRDS(system.file("examples/rbc-s3-tiledb/glasser360_fsLR164k_labels.rds", package = "ModelArray"))
+stopifnot(length(labels$key) == nrow(scalars(ma)[["thickness"]]))
+ladder <- labels$table[labels$table$region %in% p$parcels & labels$table$hemi %in% strsplit(p$hemispheres, "")[[1]], ]
+ladder <- ladder[order(ladder$sa_rank, ladder$hemi), ]
+
+pick_vertices <- function(key) {
+  idx <- which(labels$key == key)
+  tiles <- sort(table(tile_of(idx)), decreasing = TRUE)
+  keep_tiles <- as.integer(names(tiles))[seq_len(min(p$max_tiles_per_parcel, length(tiles)))]
+  idx <- idx[tile_of(idx) %in% keep_tiles]
+  if (length(idx) > p$max_vertices_per_parcel) {
+    set.seed(2026)
+    idx <- sort(sample(idx, p$max_vertices_per_parcel))
+  }
+  list(idx = idx, n_all = sum(tiles), n_tiles_all = length(tiles), tiles = sort(keep_tiles))
+}
+selection <- lapply(ladder$key, pick_vertices)
+names(selection) <- ladder$name
+ladder$vertices_in_parcel <- vapply(selection, function(s) s$n_all, numeric(1))
+ladder$tiles_in_parcel <- vapply(selection, function(s) s$n_tiles_all, numeric(1))
+ladder$tiles_read <- vapply(selection, function(s) length(s$tiles), numeric(1))
+ladder$vertices_used <- vapply(selection, function(s) length(s$idx), numeric(1))
+ladder$predicted_mb <- round(ladder$tiles_read * TILE_MB)
+ladder[, c("name", "sa_rank", "vertices_in_parcel", "tiles_in_parcel", "tiles_read", "vertices_used", "predicted_mb")]
+```
+
+## 6. Where the two parcels sit on the brain
+
+The two parcels sit at opposite ends of the sensorimotor–association (S-A) axis of
+Sydnor et al. (2021). The figure below is that axis on the left inflated cortex, in
+Sydnor's colour scheme (Figure 2A): gold at the sensorimotor pole (low rank), through
+white, to purple at the association pole (high rank). Area 3b (rank 7) is deep in the
+gold sensorimotor end; area 46 (rank 174) is deep in the purple association end. Every
+figure below colours the two parcels by these same S-A colours, so a parcel's colour
+tells you where it sits on the axis. The brain image ships with the example and is a
+left-lateral view (right-hemisphere parcels are shown at their left-hemisphere homolog).
+
+```{r fig-helpers}
+# Sydnor et al. 2021 (Neuron) Fig 2A S-A colour scheme: gold (sensorimotor) -> white -> purple (association)
+SA_COLS <- c("goldenrod1", "white", "#6f1282"); SA_MID <- 90; SA_RANGE <- c(1, 180)
+sa_rank_col <- function(r) scales::gradient_n_pal(SA_COLS)(scales::rescale_mid(r, to = c(0, 1), from = SA_RANGE, mid = SA_MID))
+parcel_regions <- unique(ladder$region[order(ladder$sa_rank)])
+parcel_ranks <- vapply(parcel_regions, function(rg) ladder$sa_rank[ladder$region == rg][1], numeric(1))
+parcel_cols <- stats::setNames(sa_rank_col(parcel_ranks), parcel_regions)
+
+# The shipped S-A brains, as overlay grobs: "axis" (3b and 46 outlined and labelled) and
+# one per parcel ("3b", "46", outlined only) for the per-panel locators.
+brain_files <- c(axis = "brain_sa_axis.png", "3b" = "brain_sa_3b.png", "46" = "brain_sa_46.png")
+brain_paths <- vapply(brain_files, function(f) system.file(file.path("examples/rbc-s3-tiledb", f), package = "ModelArray"), character(1))
+have_brain <- have_cowplot && all(nzchar(brain_paths)) && requireNamespace("png", quietly = TRUE)
+brain_grobs <- if (have_brain) lapply(brain_paths, function(p) grid::rasterGrob(png::readPNG(p), interpolate = TRUE)) else NULL
+add_brain <- function(plot, which = "axis", x = 0.0, y = 0.6, w = 0.3, h = 0.3) {
+  if (!have_brain) return(plot)
+  cowplot::ggdraw(plot) + cowplot::draw_grob(brain_grobs[[which]], x = x, y = y, width = w, height = h)
+}
+```
+
+```{r sa-axis-figure, fig.width = 7, fig.height = 3}
+if (have_brain) {
+  bar <- data.frame(rank = SA_RANGE[1]:SA_RANGE[2])
+  ticks <- data.frame(rank = parcel_ranks, lab = parcel_regions)
+  colourbar <- ggplot(bar, aes(rank, 1, fill = rank)) + geom_tile() +
+    scale_fill_gradient2(low = SA_COLS[1], mid = SA_COLS[2], high = SA_COLS[3], midpoint = SA_MID, guide = "none") +
+    geom_point(data = ticks, aes(rank, 1), inherit.aes = FALSE, shape = 25, size = 2.4, fill = "black") +
+    geom_text(data = ticks, aes(rank, 1.9, label = lab), inherit.aes = FALSE, fontface = 2, size = 3.4) +
+    annotate("text", x = SA_RANGE[1], y = -1.2, label = "sensorimotor", hjust = 0, size = 3.2) +
+    annotate("text", x = SA_RANGE[2], y = -1.2, label = "association", hjust = 1, size = 3.2) +
+    scale_y_continuous(limits = c(-2, 2.8)) + theme_void()
+  print(cowplot::plot_grid(cowplot::ggdraw() + cowplot::draw_grob(brain_grobs[["axis"]]),
+                           colourbar, ncol = 1, rel_heights = c(3, 1)))
+} else {
+  cat("Install cowplot and png to show the S-A axis brain figure.\n")
+}
+```
+
+## 7. One vertex at a time
+
+`exampleElementData()` returns the data frame a per-vertex model sees: the cohort plus
+that vertex's values. We take one vertex near the middle of the lower-ranked parcel
+(3b) and one near the middle of the higher-ranked parcel (46). The medial wall is
+stored as zeros, which is a handy check that the labels line up with the data.
+
+```{r vertex-data}
+mid_vertex <- function(sel) sel$idx[which.min(abs(sel$idx - median(sel$idx)))]
+low_name <- names(selection)[which.min(ladder$sa_rank)]
+top_name <- names(selection)[which.max(ladder$sa_rank)]
+v_low <- mid_vertex(selection[[low_name]])
+v_top <- mid_vertex(selection[[top_name]])
+low_region <- ladder$region[ladder$name == low_name]; low_rank <- ladder$sa_rank[ladder$name == low_name]
+top_region <- ladder$region[ladder$name == top_name]; top_rank <- ladder$sa_rank[ladder$name == top_name]
+
+# a medial-wall vertex in a tile we will read anyway, if there is one; otherwise the nearest one
+tiles_planned <- unique(unlist(lapply(selection, `[[`, "tiles")))
+wall_all <- which(labels$key == 0L)
+wall_in_plan <- wall_all[tile_of(wall_all) %in% tiles_planned]
+wall <- if (length(wall_in_plan)) wall_in_plan[1] else wall_all[which.min(abs(wall_all - v_low))]
+n_tiles <- length(unique(tile_of(c(v_low, wall))))
+
+two <- tdb_measure(sprintf("1 %s vertex + 1 medial-wall vertex (%d tile%s)", low_region, n_tiles, if (n_tiles > 1) "s" else ""),
+                   as.matrix(scalars(ma)[["thickness"]][c(v_low, wall), ]))
+cat(sprintf("%s vertex %d: mean thickness %.2f mm; medial-wall vertex %d: all zero = %s (%s, %.0f MB)\n",
+            low_region, v_low, mean(two[1, ]), wall, all(two[2, ] == 0), fmt_secs(last_runtime()$secs), last_runtime()$mb_wire))
+
+df_top <- exampleElementData(ma, scalar = "thickness", i_element = v_top, phenotypes = phen)
+df_low <- exampleElementData(ma, scalar = "thickness", i_element = v_low, phenotypes = phen)
+str(df_top[, c("dataset", "Age", "sex", "site", "thickness")])
+```
+
+```{r vertex-figure, fig.width = 9, fig.height = 4.5}
+low_lab <- sprintf("%s (S-A rank %s)", low_region, low_rank)
+top_lab <- sprintf("%s (S-A rank %s)", top_region, top_rank)
+both <- rbind(cbind(df_low, parcel = low_lab), cbind(df_top, parcel = top_lab))
+both$parcel <- factor(both$parcel, levels = c(low_lab, top_lab))
+fig1 <- ggplot(both, aes(Age, thickness)) +
+  geom_point(aes(colour = dataset), alpha = 0.15, size = 0.6) +
+  geom_smooth(method = "gam", formula = y ~ s(x, k = 4), colour = "black") +
+  facet_wrap(~parcel) +
+  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1))) +
+  labs(x = "Age (years)", y = "Cortical thickness (mm)", colour = "Study",
+       title = sprintf("Thickness against age at one vertex in each parcel (all %s sessions)", format(nrow(phen), big.mark = ","))) +
+  theme_minimal() + theme(plot.title = element_text(size = 11))
+# a small brain to the left of each panel's title: 3b by the 3b panel, 46 by the 46 panel
+if (have_brain && all(c(low_region, top_region) %in% names(brain_grobs))) {
+  fig1 <- cowplot::ggdraw(fig1) +
+    cowplot::draw_grob(brain_grobs[[low_region]], x = 0.10, y = 0.80, width = 0.11, height = 0.15) +
+    cowplot::draw_grob(brain_grobs[[top_region]], x = 0.47, y = 0.80, width = 0.11, height = 0.15)
+}
+fig1
+ggsave(file.path(p$out_dir, "fig1_vertex_thickness_vs_age.png"), fig1, width = 9, height = 4.5, dpi = 150)
+```
+
+A plain GAM at the area-`r top_region` vertex over the whole age range, with site as
+a covariate, shows the lifespan picture: thickness peaks in early childhood and
+declines from then on (Bethlehem et al., 2022), so the slope should be negative in
+old age too (ages above 23 come from NKI only). Two checks go with it: the derivative
+of the age curve computed two ways, and a mixed model that treats repeat sessions
+(two waves in BHRC, several in NKI) properly, which the vertex-wise GAM does not.
+
+```{r vertex-gam}
+AGE_GRID <- seq(5, 85, by = 0.5)
+
+# Slope of the age smooth by central finite differences, with a simultaneous 95%
+# interval from posterior simulation (the same construction as gratia::derivatives()).
+smooth_derivative <- function(fit, newdata, var = "Age", eps = 0.05, nsim = 1000) {
+  lo <- newdata; lo[[var]] <- newdata[[var]] - eps
+  hi <- newdata; hi[[var]] <- newdata[[var]] + eps
+  X <- (predict(fit, hi, type = "lpmatrix") - predict(fit, lo, type = "lpmatrix")) / (2 * eps)
+  X[, !grepl(paste0("s(", var, ")"), colnames(X), fixed = TRUE)] <- 0
+  b <- coef(fit)
+  Vb <- vcov(fit)
+  d <- drop(X %*% b)
+  se <- sqrt(rowSums((X %*% Vb) * X))
+  sims <- X %*% t(MASS::mvrnorm(nsim, b, Vb))
+  crit <- stats::quantile(apply(abs(sims - d) / se, 2, max), 0.95)
+  data.frame(Age = newdata[[var]], derivative = d, lower = d - crit * se, upper = d + crit * se)
+}
+# A data frame of typical covariate values (most common level of each factor) along the age grid.
+typical_newdata <- function(d, cov_cols, ages = AGE_GRID) {
+  ref <- d[1, cov_cols, drop = FALSE]
+  for (v in setdiff(cov_cols, "Age")) {
+    ref[[v]] <- factor(names(which.max(table(d[[v]]))), levels = levels(d[[v]]))
+  }
+  nd <- ref[rep(1, length(ages)), , drop = FALSE]
+  nd$Age <- ages
+  nd
+}
+
+d_demo <- df_top[stats::complete.cases(df_top[, c("thickness", "Age", "sex", "site")]), ]
+d_demo$site <- droplevels(d_demo$site)
+f_demo <- stats::as.formula(paste0("thickness ~ s(Age, k = 4)", sex_term, site_term))
+fit_demo <- mgcv::gam(f_demo, data = d_demo, method = "REML")
+summary(fit_demo)$s.table
+nd_demo <- typical_newdata(d_demo, c("Age", if (use_sex) "sex", if (nzchar(site_term)) "site"))
+deriv_demo <- smooth_derivative(fit_demo, nd_demo)
+cat(sprintf("Mean slope 8-22: %.1f microns/yr; 60-80: %.1f microns/yr\n",
+            1000 * mean(deriv_demo$derivative[deriv_demo$Age >= 8 & deriv_demo$Age <= 22]),
+            1000 * mean(deriv_demo$derivative[deriv_demo$Age >= 60 & deriv_demo$Age <= 80])))
+if (have_gratia) {
+  g <- gratia::derivatives(fit_demo, select = "s(Age)", n = 50, interval = "simultaneous")
+  ours <- smooth_derivative(fit_demo, typical_newdata(d_demo, names(nd_demo), ages = g$Age))
+  cat(sprintf("Derivative check against gratia::derivatives(): max |difference| = %.1e mm/yr\n", max(abs(ours$derivative - g$.derivative))))
+}
+if (have_lme4) {
+  f_mixed <- stats::as.formula(paste0("thickness ~ Age", sex_term, site_term, " + (1 | subject)"))
+  mixed <- lme4::lmer(f_mixed, data = d_demo)
+  cat(sprintf("Mixed model (%d subjects with repeat sessions): linear age slope = %.1f microns/yr (SE %.1f)\n",
+              sum(table(d_demo$subject) > 1), 1000 * lme4::fixef(mixed)[["Age"]], 1000 * sqrt(diag(as.matrix(vcov(mixed))))[["Age"]]))
+}
+```
+
+## 8. Pull each parcel from S3 and harmonize it
+
+Nine scanning sites contributed, and they differ in the mean and the spread of
+thickness, so a `site` covariate removes only part of the difference. CovBat
+(`ComBatFamily::covfam()`) harmonizes mean, variance and the covariance between
+features across sites while keeping the age and sex effects (`y ~ s(Age, k = 4) + sex`).
+Its features here are the vertices of a parcel, read from S3 in one go: that is
+what the method is built for, and it is far cheaper than harmonizing vertex by vertex.
+
+The raw and harmonized blocks go into a small local ModelArray store, one row per
+selected vertex, so that the per-vertex models in section 8 run with the usual
+ModelArray functions. Nothing is written back to S3.
+
+```{r parcel-harmonize}
+local_h5 <- file.path(p$out_dir, "rbc_parcels_local.h5")
+if (file.exists(local_h5)) unlink(local_h5)
+write_local_scalar <- function(file, name, mat, column_names) {
+  if (!file.exists(file)) rhdf5::h5createFile(file)
+  if (!"scalars" %in% rhdf5::h5ls(file, recursive = FALSE)$name) rhdf5::h5createGroup(file, "scalars")
+  rhdf5::h5createGroup(file, paste0("scalars/", name))
+  # one row (element) per chunk: ModelArray reads elements one at a time
+  rhdf5::h5createDataset(file, paste0("scalars/", name, "/values"), dims = dim(mat), storage.mode = "double",
+                         chunk = c(1L, ncol(mat)), level = 4L)
+  rhdf5::h5write(mat, file, paste0("scalars/", name, "/values"))
+  rhdf5::h5write(column_names, file, paste0("scalars/", name, "/column_names"))
+  rhdf5::h5closeAll()
+}
+
+# Read one parcel at a time, so each read can be checked against the tile prediction
+raw_blocks <- list(); elements <- list()
+for (nm in names(selection)) {
+  sel <- selection[[nm]]
+  block <- tdb_measure(sprintf("read %s (%d vertices, %d tile%s)", nm, length(sel$idx), length(sel$tiles), if (length(sel$tiles) > 1) "s" else ""),
+                       as.matrix(scalars(ma)[["thickness"]][sel$idx, ]))
+  r <- last_runtime()
+  cat(sprintf("%-12s %4d vertices  %5.1f s  tiles read %2.0f (predicted %2d) %s  %4.0f MB on the wire\n",
+              nm, nrow(block), r$secs, r$tiles, length(sel$tiles),
+              if (is.na(r$tiles) || r$tiles <= length(sel$tiles)) "PASS" else "FAIL", r$mb_wire))
+  usable <- apply(block, 1, function(v) all(is.finite(v)) && stats::sd(v) > 0)
+  raw_blocks[[nm]] <- block[usable, , drop = FALSE]
+  elements[[nm]] <- data.frame(vertex = sel$idx[usable], name = nm, stringsAsFactors = FALSE)
+}
+```
+
+CovBat fits one GAM per vertex and keeps every fit. Each fit carries copies of the
+data (about 5 MB at this sample size) that `covfam()` never uses again, so
+`slim_gam()` drops them; the harmonized values are identical. Parcels are harmonized
+independently, so they can run at the same time, one per core (forked processes;
+on Windows they run one after another). Each worker needs about 1 GB per 300 vertices.
+
+```{r parcel-covbat}
+ok_batch <- stats::complete.cases(phen[, c("site", harmonize_covariates)])
+d_batch <- droplevels(phen[ok_batch, , drop = FALSE])
+slim_gam <- function(...) {
+  fit <- mgcv::gam(...)
+  fit[c("call", "model", "y", "fitted.values", "residuals", "linear.predictors", "weights", "prior.weights",
+        "working.weights", "hat", "offset", "R", "qrx", "db.drho", "dw.drho", "F")] <- NULL
+  fit
+}
+covbat_block <- function(block) {
+  h <- suppressWarnings(ComBatFamily::covfam(
+    data = t(block[, ok_batch, drop = FALSE]), bat = d_batch$site,
+    covar = d_batch[, harmonize_covariates, drop = FALSE], model = slim_gam,
+    formula = stats::as.formula(paste0("y ~ s(Age, k = 4)", sex_term)), eb = TRUE))
+  harm <- matrix(NA_real_, nrow(block), ncol(block))
+  harm[, ok_batch] <- t(as.matrix(h$dat.covbat))
+  harm
+}
+covbat_cores <- if (.Platform$OS.type == "windows") 1L else max(1L, min(as.integer(p$n_cores), length(raw_blocks)))
+if (harmonize) {
+  harm_blocks <- timed(sprintf("CovBat, %d parcels on %d core%s", length(raw_blocks), covbat_cores, if (covbat_cores > 1) "s" else ""),
+    if (covbat_cores > 1) {
+      parallel::mclapply(raw_blocks, covbat_block, mc.cores = covbat_cores, mc.preschedule = FALSE)
+    } else {
+      lapply(raw_blocks, covbat_block)
+    })
+  failed <- vapply(harm_blocks, function(x) inherits(x, "try-error") || is.null(x), logical(1))
+  if (any(failed)) {
+    stop("CovBat failed for ", paste(names(harm_blocks)[failed], collapse = ", "), ": ",
+         paste(unique(vapply(harm_blocks[failed], function(x) if (is.null(x)) "worker returned nothing (out of memory?)" else as.character(x), "")), collapse = "; "),
+         "\nTry a smaller n_cores or max_vertices_per_parcel.", call. = FALSE)
+  }
+  cat(sprintf("CovBat on %d sessions x %d vertices (%d parcels, %d core%s) in %s\n", sum(ok_batch), sum(vapply(raw_blocks, nrow, integer(1))),
+              length(raw_blocks), covbat_cores, if (covbat_cores > 1) "s" else "", fmt_secs(last_runtime()$secs)))
+} else {
+  harm_blocks <- lapply(raw_blocks, function(b) matrix(NA_real_, nrow(b), ncol(b)))
+}
+
+elements <- do.call(rbind, elements)
+elements$element_id <- seq_len(nrow(elements)) - 1L               # ModelArray's 0-based ids in the local store
+elements <- merge(elements, ladder[, c("name", "hemi", "region", "sa_rank")], by = "name", sort = FALSE)
+elements <- elements[order(elements$element_id), c("element_id", "vertex", "name", "hemi", "region", "sa_rank")]
+rownames(elements) <- NULL
+write_local_scalar(local_h5, "thickness_raw", do.call(rbind, raw_blocks), src)
+if (harmonize) write_local_scalar(local_h5, "thickness_covbat", do.call(rbind, harm_blocks), src)
+utils::write.csv(elements, file.path(p$out_dir, "elements.csv"), row.names = FALSE)
+cat(sprintf("\nLocal store %s: %d vertices x %d sessions, %.0f MB\n", basename(local_h5), nrow(elements), length(src), file.size(local_h5) / 1024^2))
+```
+
+```{r site-figure, fig.width = 9, fig.height = 4}
+if (harmonize) {
+  row_top <- match(v_top, elements$vertex[elements$name == top_name])
+  dh <- data.frame(site = phen$site, Age = phen$Age,
+                   raw = raw_blocks[[top_name]][row_top, ],
+                   harmonized = harm_blocks[[top_name]][row_top, ])
+  dh <- dh[ok_batch, ]
+  dh$site <- droplevels(dh$site)
+  site_stats <- data.frame(site = levels(dh$site),
+                           n = as.vector(table(dh$site)),
+                           mean_raw = round(tapply(dh$raw, dh$site, mean), 3),
+                           sd_raw = round(tapply(dh$raw, dh$site, sd), 3),
+                           mean_harmonized = round(tapply(dh$harmonized, dh$site, mean), 3),
+                           sd_harmonized = round(tapply(dh$harmonized, dh$site, sd), 3), row.names = NULL)
+  site_stats
+  long <- rbind(data.frame(site = dh$site, thickness = dh$raw, stage = "raw"),
+                data.frame(site = dh$site, thickness = dh$harmonized, stage = "CovBat-harmonized"))
+  long$stage <- factor(long$stage, levels = c("raw", "CovBat-harmonized"))
+  fig_sites <- ggplot(long, aes(site, thickness, fill = site)) +
+    geom_boxplot(outlier.size = 0.3, linewidth = 0.3) + facet_wrap(~stage) +
+    scale_fill_brewer(palette = "Set3", guide = "none") +
+    labs(x = NULL, y = "Thickness (mm)",
+         title = sprintf("Area %s vertex %d: thickness by site, before and after harmonization", top_region, v_top)) +
+    theme_minimal() + theme(axis.text.x = element_text(angle = 45, hjust = 1), plot.title = element_text(size = 11))
+  fig_sites
+  ggsave(file.path(p$out_dir, "fig_site_boxplots.png"), fig_sites, width = 9, height = 4, dpi = 150)
+} else {
+  cat("Harmonization is off (ComBatFamily not installed or harmonize = FALSE); site enters the models as a covariate.\n")
+}
+```
+
+Site means still differ a little after harmonization because the age and sex mix
+differs between sites, and those effects are kept on purpose.
+
+## 9. Models at every vertex
+
+The developmental question is about ages 5–`r p$max_age`, which is where four of the
+five studies sit. ModelArray needs the phenotype table to list every session in the
+store, so instead of dropping the older sessions we blank their age; `mgcv` then
+leaves them out of every fit.
+
+### The built-in GAM runner, from S3 and from the local store
+
+`ModelArray.gam()` fits a GAM at every element. First straight from S3 on one parcel,
+with site as a covariate (no harmonization); then on the local store for every parcel,
+on the raw and on the harmonized values. `changed.rsq.term.index = 1` adds the partial
+R² of the age smooth. The three runs give a feel for what harmonization changes.
+(The S3 run reads more tiles than section 7 did for the same vertices: each worker
+fetches its own block of elements, so a tile shared between workers is fetched more
+than once. When a block fits in memory, reading it once yourself is cheaper.)
+
+```{r parcel-gam}
+phen_dev <- phen
+phen_dev$Age[phen_dev$Age > p$max_age] <- NA
+cat(sprintf("%d of %d sessions are aged %s or younger\n", sum(!is.na(phen_dev$Age)), nrow(phen_dev), p$max_age))
+
+gam_args <- list(changed.rsq.term.index = list(1), method = "REML", n_cores = p$n_cores,
+                 on_error = "skip", pbar = FALSE, verbose = FALSE)
+demo_name <- names(selection)[1]
+demo_idx <- elements$vertex[elements$name == demo_name]
+gam_s3 <- tdb_measure(sprintf("ModelArray.gam from S3: %s (%d vertices)", demo_name, length(demo_idx)),
+  do.call(ModelArray.gam, c(list(stats::as.formula(paste0("thickness ~ s(Age, k = 4)", sex_term, site_term)),
+                                 ma, phen_dev, "thickness", element.subset = as.integer(demo_idx)), gam_args)))
+cat(sprintf("S3, %s: %d vertices in %s (%.0f tiles)\n", demo_name, nrow(gam_s3), fmt_secs(last_runtime()$secs), last_runtime()$tiles))
+
+local_scalars <- c("thickness_raw", if (harmonize) "thickness_covbat")
+ma_local <- ModelArray(local_h5, scalar_types = local_scalars)
+ma_local
+gam_local <- list()
+for (sc in local_scalars) {
+  f <- stats::as.formula(paste0(sc, " ~ s(Age, k = 4)", sex_term, if (sc == "thickness_raw") site_term))
+  gam_local[[sc]] <- timed(sprintf("ModelArray.gam local: %s (%d vertices)", sc, nrow(elements)),
+    do.call(ModelArray.gam, c(list(f, ma_local, phen_dev, sc), gam_args)))
+  cat(sprintf("local, %s: %d vertices in %s\n", sc, nrow(gam_local[[sc]]), fmt_secs(last_runtime()$secs)))
+}
+
+gam_summary <- function(res, label) {
+  rsq <- intersect(c("s_Age.partial.rsq", "s_Age.delta.adj.rsq"), names(res))[1]
+  data.frame(run = label, vertices = nrow(res),
+             pct_fdr_sig = round(100 * mean(res$s_Age.p.value.fdr < 0.05, na.rm = TRUE)),
+             median_F = round(median(res$s_Age.statistic, na.rm = TRUE), 1),
+             median_partial_rsq = round(median(res[[rsq]], na.rm = TRUE), 4))
+}
+demo_rows <- elements$element_id[elements$name == demo_name]
+rbind(gam_summary(gam_s3, sprintf("S3, raw + site covariate (%s)", demo_name)),
+      do.call(rbind, lapply(local_scalars, function(sc)
+        gam_summary(gam_local[[sc]][gam_local[[sc]]$element_id %in% demo_rows, ],
+                    sprintf("local, %s%s (%s)", sc, if (sc == "thickness_raw") " + site covariate" else "", demo_name)))))
+```
+
+### A custom per-vertex function with `ModelArray.wrap()`
+
+`fit_vertex()` adds what `ModelArray.gam()` does not report: a test statistic from an
+`fx = TRUE` refit (penalized p-values can be optimistic), ΔR²adj from a reduced model,
+predicted thickness every two years, and the slope of the age curve with a simultaneous
+95% interval. From the slope it takes the mean thinning rate over 8–22, the age at
+which thinning is fastest, the last age at which the slope is still clearly negative,
+and the **late fraction**: the share of the 8–22 thinning that happens after age 14.
+One row of numbers per vertex.
+
+```{r fit-fun}
+response_scalar <- if (harmonize) "thickness_covbat" else "thickness_raw"
+DEV_GRID <- seq(5, p$max_age, by = 0.5)
+PRED_AGES <- seq(6, 22, by = 2)
+fit_vertex <- function(data) {
+  d <- data
+  d$y <- d[[response_scalar]]
+  cov_cols <- c("Age", if (use_sex) "sex", if (!harmonize && nzchar(site_term)) "site")
+  d <- d[stats::complete.cases(d[, c("y", cov_cols)]), ]
+  for (v in setdiff(cov_cols, "Age")) d[[v]] <- droplevels(d[[v]])
+  rhs <- if (length(cov_cols) > 1) paste(cov_cols[-1], collapse = " + ") else "1"
+  fit <- mgcv::gam(stats::as.formula(paste("y ~ s(Age, k = 4) +", rhs)), data = d, method = "REML")
+  fit_fx <- mgcv::gam(stats::as.formula(paste("y ~ s(Age, k = 4, fx = TRUE) +", rhs)), data = d, method = "REML")
+  fit_red <- mgcv::gam(stats::as.formula(paste("y ~", rhs)), data = d, method = "REML")
+  s_pen <- summary(fit)$s.table
+  s_fx <- summary(fit_fx)$s.table
+
+  newdata <- typical_newdata(d, cov_cols, ages = DEV_GRID)
+  pred <- predict(fit, newdata[match(PRED_AGES, DEV_GRID), , drop = FALSE])
+  names(pred) <- PRED_AGES
+  deriv <- smooth_derivative(fit, newdata)
+  in_range <- function(a, b) deriv$Age >= a & deriv$Age <= b
+  sig_ages <- deriv$Age[deriv$upper < 0]
+  out <- data.frame(
+    n = nrow(d),
+    edf = unname(s_pen[1, "edf"]),
+    F_fx = unname(s_fx[1, "F"]),
+    p_fx = unname(s_fx[1, "p-value"]),
+    p_pen = unname(s_pen[1, "p-value"]),
+    delta_adj_rsq = summary(fit)$r.sq - summary(fit_red)$r.sq,
+    rate_8_22 = mean(deriv$derivative[in_range(8, 22)]),
+    rate_8_14 = mean(deriv$derivative[in_range(8, 14)]),
+    rate_14_22 = mean(deriv$derivative[in_range(14, 22)]),
+    late_fraction = unname((pred["14"] - pred["22"]) / (pred["8"] - pred["22"])),
+    age_fastest = deriv$Age[which.min(deriv$derivative)],
+    thinning_until = if (length(sig_ages)) max(sig_ages) else NA_real_,
+    sex_estimate = if (use_sex) unname(coef(fit)[grep("^sex", names(coef(fit)))[1]]) else NA_real_
+  )
+  out[paste0("pred_", PRED_AGES)] <- as.list(unname(pred))
+  out[paste0("slope_", PRED_AGES)] <- as.list(deriv$derivative[match(PRED_AGES, deriv$Age)])
+  out
+}
+
+# Try it on one element the way ModelArray.wrap() will call it
+one <- exampleElementData(ma_local, scalar = response_scalar, i_element = which(elements$vertex == v_top), phenotypes = phen_dev)
+round(t(fit_vertex(one)), 4)
+check <- analyseOneElement.wrap(which(elements$vertex == v_top), fit_vertex, ma_local, phen_dev, response_scalar,
+                                num.subj.lthr = 10, flag_initiate = TRUE, on_error = "stop")
+check$column_names
+```
+
+```{r parcel-wrap}
+vertex_results <- timed(sprintf("ModelArray.wrap fit_vertex (%d vertices)", nrow(elements)),
+  ModelArray.wrap(fit_vertex, ma_local, phen_dev, response_scalar, n_cores = p$n_cores,
+                  on_error = "skip", pbar = FALSE, verbose = FALSE))
+cat(sprintf("%d vertices in %s (%.2f s per vertex per core)\n", nrow(vertex_results), fmt_secs(last_runtime()$secs),
+            last_runtime()$secs * p$n_cores / nrow(vertex_results)))
+vertex_results <- merge(elements, vertex_results, by = "element_id", sort = TRUE)
+vertex_results$p_fx_fdr <- stats::p.adjust(vertex_results$p_fx, method = "fdr")
+head(vertex_results[, c("element_id", "vertex", "name", "edf", "F_fx", "p_fx_fdr", "delta_adj_rsq", "rate_8_22", "late_fraction", "age_fastest")])
+```
+
+### Repeat sessions, at one vertex
+
+CovBat and the vertex-wise GAM treat repeat sessions as independent. A GAMM with a
+subject random intercept (`gamm4`) is the honest alternative for the model, and
+longitudinal ComBat (the `longCombat` package) for the harmonization. The GAMM is far
+slower than a GAM, so the report fits it at the area-`r top_region` vertex only, with
+site as a covariate, and gives up after two minutes.
+
+```{r vertex-gamm}
+if (have_gamm4) {
+  d_long <- d_demo
+  d_long$subject <- factor(d_long$subject)
+  t0 <- Sys.time()
+  gamm_fit <- tryCatch({
+    setTimeLimit(elapsed = 120, transient = TRUE)
+    gamm4::gamm4(f_demo, random = ~ (1 | subject), data = d_long, REML = TRUE)
+  }, error = function(e) e)
+  setTimeLimit(elapsed = Inf)
+  if (inherits(gamm_fit, "error")) {
+    cat("The GAMM did not finish within 2 minutes at one vertex; skipped. (", conditionMessage(gamm_fit), ")\n")
+  } else {
+    nd_gamm <- typical_newdata(d_long, names(nd_demo))
+    deriv_gamm <- smooth_derivative(gamm_fit$gam, nd_gamm)
+    cat(sprintf("GAMM at one vertex took %s; smooth edf = %.2f (GAM: %.2f); mean slope 8-22 = %.1f microns/yr (GAM: %.1f)\n",
+                fmt_secs(as.numeric(Sys.time() - t0, units = "secs")),
+                summary(gamm_fit$gam)$s.table[1, "edf"], summary(fit_demo)$s.table[1, "edf"],
+                1000 * mean(deriv_gamm$derivative[deriv_gamm$Age >= 8 & deriv_gamm$Age <= 22]),
+                1000 * mean(deriv_demo$derivative[deriv_demo$Age >= 8 & deriv_demo$Age <= 22])))
+  }
+}
+```
+
+## 10. Does the result match the literature?
+
+```{r parcel-summary}
+q <- function(x, pr) unname(stats::quantile(x, pr, na.rm = TRUE))
+summ <- do.call(rbind, lapply(split(vertex_results, vertex_results$name), function(v) {
+  data.frame(name = v$name[1], hemi = v$hemi[1], region = v$region[1], sa_rank = v$sa_rank[1],
+             vertices = sum(!is.na(v$late_fraction)),
+             pct_fdr_sig = round(100 * mean(v$p_fx_fdr < 0.05, na.rm = TRUE)),
+             thickness_8 = median(v$pred_8, na.rm = TRUE),
+             rate_8_22 = median(v$rate_8_22, na.rm = TRUE),
+             pct_per_yr_8_22 = median(100 * v$rate_8_22 / v$pred_8, na.rm = TRUE),
+             rate_8_14 = median(v$rate_8_14, na.rm = TRUE),
+             rate_14_22 = median(v$rate_14_22, na.rm = TRUE),
+             late_fraction = median(v$late_fraction, na.rm = TRUE),
+             late_fraction_q25 = q(v$late_fraction, 0.25),
+             late_fraction_q75 = q(v$late_fraction, 0.75),
+             age_fastest = median(v$age_fastest, na.rm = TRUE),
+             thinning_until = median(v$thinning_until, na.rm = TRUE),
+             delta_adj_rsq = median(v$delta_adj_rsq, na.rm = TRUE))
+}))
+summ <- summ[order(summ$sa_rank, summ$hemi), ]
+rownames(summ) <- NULL
+summ_print <- summ
+for (col in c("rate_8_22", "rate_8_14", "rate_14_22")) summ_print[[col]] <- round(1000 * summ_print[[col]], 1)
+names(summ_print) <- sub("^rate_", "um_per_yr_", names(summ_print))
+for (col in c("thickness_8", "pct_per_yr_8_22", "late_fraction", "late_fraction_q25", "late_fraction_q75")) summ_print[[col]] <- round(summ_print[[col]], 2)
+summ_print$delta_adj_rsq <- round(summ_print$delta_adj_rsq, 4)
+summ_print
+
+if (nrow(summ) > 2) {
+  rho_late <- stats::cor(summ$sa_rank, summ$late_fraction, method = "spearman")
+  rho_age <- stats::cor(summ$sa_rank, summ$age_fastest, method = "spearman")
+  cat(sprintf("\nSpearman correlation with S-A rank across %d parcel(s): late fraction rho = %.2f, age of fastest thinning rho = %.2f (both expected positive)\n",
+              nrow(summ), rho_late, rho_age))
+}
+cat(sprintf("All parcels thin between 8 and 22: %s. Earliest-maturing parcel (smallest late fraction): %s (expected: 3b). Latest: %s (expected: the highest rank, %s).\n",
+            all(summ$rate_8_22 < 0), summ$region[which.min(summ$late_fraction)], summ$region[which.max(summ$late_fraction)], summ$region[which.max(summ$sa_rank)]))
+```
+
+```{r parcel-figures, fig.width = 8, fig.height = 4.5}
+# parcels coloured by their S-A position (gold = sensorimotor, purple = association); L/R nudged apart
+summ$xpos <- summ$sa_rank + ifelse(summ$hemi == "L", -4, 4)
+fig2 <- ggplot(summ, aes(xpos, late_fraction, colour = region, shape = hemi)) +
+  geom_errorbar(aes(ymin = late_fraction_q25, ymax = late_fraction_q75), width = 5, linewidth = 0.5) +
+  geom_point(size = 3.4, stroke = 0.7) +
+  geom_text(aes(label = region), vjust = -1.4, size = 3.4, show.legend = FALSE) +
+  scale_colour_manual(values = parcel_cols, guide = "none") +
+  scale_shape_manual(values = c(L = 16, R = 17)) +
+  expand_limits(x = c(0, 180)) +
+  labs(x = "Sensorimotor-association axis rank (1 = sensorimotor, 180 = association)",
+       y = "Share of 8-22 thinning after age 14", shape = "Hemisphere",
+       title = "Later-maturing parcels do more of their thinning late (median, IQR across vertices)") +
+  theme_minimal() + theme(plot.title = element_text(size = 11))
+fig2 <- add_brain(fig2, x = 0.36, y = 0.6, w = 0.22, h = 0.32)
+fig2
+ggsave(file.path(p$out_dir, "fig2_late_fraction_vs_sa_rank.png"), fig2, width = 8, height = 4.5, dpi = 150)
+
+curves <- do.call(rbind, lapply(split(vertex_results, vertex_results$name), function(v) {
+  data.frame(name = v$name[1], region = v$region[1], hemi = v$hemi[1], sa_rank = v$sa_rank[1], Age = PRED_AGES,
+             thickness = vapply(paste0("pred_", PRED_AGES), function(col) median(v[[col]], na.rm = TRUE), numeric(1)),
+             slope = vapply(paste0("slope_", PRED_AGES), function(col) median(v[[col]], na.rm = TRUE), numeric(1)))
+}))
+curves$parcel <- factor(curves$region, levels = parcel_regions)
+fig3 <- ggplot(curves, aes(Age, thickness, colour = parcel, linetype = hemi)) +
+  geom_line(linewidth = 0.9) + geom_point(aes(fill = parcel), shape = 21, colour = "black", size = 2, stroke = 0.4) +
+  scale_colour_manual(values = parcel_cols) + scale_fill_manual(values = parcel_cols, guide = "none") +
+  labs(x = "Age (years)", y = "Predicted thickness (mm), parcel median", colour = "Parcel", linetype = "Hemisphere",
+       title = "Fitted age trajectories by parcel") + theme_minimal() + theme(plot.title = element_text(size = 11))
+fig3
+ggsave(file.path(p$out_dir, "fig3_parcel_trajectories.png"), fig3, width = 8, height = 4.5, dpi = 150)
+fig4 <- ggplot(curves, aes(Age, 1000 * slope, colour = parcel, linetype = hemi)) +
+  geom_hline(yintercept = 0, linetype = 2, colour = "grey50") +
+  geom_line(linewidth = 1) + geom_point(aes(fill = parcel), shape = 21, colour = "black", size = 2.2, stroke = 0.5) +
+  scale_colour_manual(values = parcel_cols) + scale_fill_manual(values = parcel_cols, guide = "none") +
+  labs(x = "Age (years)", y = "Slope of the age curve (microns/year), parcel median", colour = "Parcel", linetype = "Hemisphere",
+       title = "When each parcel thins: slope of the thickness-age curve") +
+  theme_minimal() + theme(plot.title = element_text(size = 11))
+fig4 <- add_brain(fig4, x = 0.36, y = 0.58, w = 0.22, h = 0.32)
+fig4
+ggsave(file.path(p$out_dir, "fig4_parcel_slopes.png"), fig4, width = 8, height = 4.5, dpi = 150)
+```
+
+What to expect, from the literature. Shaw et al. (2008, *J Neurosci*) found that
+primary sensory and motor cortex follow simple, early trajectories while association
+cortex follows later and more drawn-out ones. Tamnes et al. (2017, *J Neurosci*),
+studying 8–29-year-olds, found about 1% thinning per year, fastest in prefrontal and
+parietal association cortex and slow in sensorimotor cortex. Sydnor et al. (2021,
+*Neuron*) put these regional differences on one S-A axis along which association
+cortex matures last. Together they predict what the slope-by-age figure should show:
+area 3b (sensorimotor, low rank) thinning fastest in childhood and flattening by
+mid-adolescence, prefrontal area 46 (association, high rank) still thinning in the
+twenties, and a higher late fraction for 46 than for 3b.
+
+## 11. Save the results locally
+
+Results never go back to S3. They are written as CSVs, and the per-vertex table also
+goes into the local store in ModelArray's own results layout, next to the scalars it
+came from.
+
+```{r write-results}
+utils::write.csv(vertex_results, file.path(p$out_dir, "vertex_results.csv"), row.names = FALSE)
+utils::write.csv(summ, file.path(p$out_dir, "parcel_summary.csv"), row.names = FALSE)
+numeric_cols <- vapply(vertex_results, is.numeric, logical(1))
+writeResults(local_h5, df.output = vertex_results[, numeric_cols], analysis_name = "gam_thickness_sa_ladder")
+ModelArray(local_h5, scalar_types = response_scalar, analysis_names = "gam_thickness_sa_ladder")
+cat("Wrote", file.path(p$out_dir, "vertex_results.csv"), "\n      ", file.path(p$out_dir, "parcel_summary.csv"), "\n      ", local_h5, "\n")
+```
+
+## 12. Look at the maps in Connectome Workbench
+
+The per-vertex results become a CIFTI dscalar on the fsLR-164k mesh: one map per
+statistic at the analysed vertices, and parcel-wide maps that paint each parcel with
+its median thinning rate and its S-A rank. Vertices that were not analysed are NaN.
+Workbench needs surfaces to draw on; TemplateFlow publishes fsLR surfaces at 32k, which
+are upsampled to 164k here, and a `.spec` file ties everything together.
+
+```{r cifti-export}
+if (!have_wb) {
+  cat("Connectome Workbench (wb_command) was not found, so no CIFTI files were written.",
+      "Install Workbench and put wb_command on your PATH (or set write_cifti = TRUE) to get them.\n")
+} else {
+  n_vertices <- length(labels$key)
+  maps <- c("sAge_F", "sAge_neglog10_fdr_p", "rate_8_22_mm_per_yr", "late_fraction", "age_fastest",
+            "parcel_median_late_fraction", "parcel_median_rate_8_22", "parcel_SA_rank")
+  m <- matrix(NA_real_, nrow = n_vertices, ncol = length(maps), dimnames = list(NULL, maps))
+  i <- vertex_results$vertex
+  m[i, "sAge_F"] <- vertex_results$F_fx
+  m[i, "sAge_neglog10_fdr_p"] <- -log10(vertex_results$p_fx_fdr)
+  m[i, "rate_8_22_mm_per_yr"] <- vertex_results$rate_8_22
+  m[i, "late_fraction"] <- vertex_results$late_fraction
+  m[i, "age_fastest"] <- vertex_results$age_fastest
+  for (k in seq_len(nrow(summ))) {
+    in_parcel <- labels$key == labels$table$key[labels$table$name == summ$name[k]]
+    m[in_parcel, "parcel_median_late_fraction"] <- summ$late_fraction[k]
+    m[in_parcel, "parcel_median_rate_8_22"] <- summ$rate_8_22[k]
+    m[in_parcel, "parcel_SA_rank"] <- summ$sa_rank[k]
+  }
+  txt <- file.path(p$out_dir, "results_fsLR164k.txt")
+  utils::write.table(m, txt, sep = "\t", row.names = FALSE, col.names = FALSE, na = "NaN")
+  names_file <- file.path(p$out_dir, "map_names.txt")
+  writeLines(maps, names_file)
+
+  template <- cached_download("https://fcp-indi.s3.amazonaws.com/data/Projects/RBC/CCNP_FreeSurfer-Post/SHA256E-s3385504--99c174365d91c22508a0de725b40e19f8e6b572c6efd88dfcd681e392e187267.nii",
+                              "template_fsLR164k.dscalar.nii")
+  stopifnot(file.size(template) == 3385504)
+  dscalar <- file.path(p$out_dir, "results_fsLR164k.dscalar.nii")
+  converted <- tryCatch({
+    wb_run("-cifti-convert", "-from-text", txt, template, dscalar, "-reset-scalars"); TRUE
+  }, error = function(e) FALSE)
+  if (!converted) {
+    # This Workbench build did not accept "NaN" tokens: use a sentinel and turn it into NaN with 0/0
+    utils::write.table(m, txt, sep = "\t", row.names = FALSE, col.names = FALSE, na = "-99999")
+    tmp <- file.path(p$out_dir, "results_fsLR164k_sentinel.dscalar.nii")
+    wb_run("-cifti-convert", "-from-text", txt, template, tmp, "-reset-scalars")
+    wb_run("-cifti-math", "(x > -90000) * x / (x > -90000)", dscalar, "-var", "x", tmp)
+    unlink(tmp)
+  }
+  wb_run("-set-map-names", dscalar, "-name-file", names_file)
+  cat("Wrote", dscalar, "with maps:", paste(maps, collapse = ", "), "\n")
+}
+```
+
+```{r cifti-display}
+if (have_wb) {
+  tf <- "https://templateflow.s3.amazonaws.com/tpl-fsLR"
+  spheres <- list()
+  for (h in c("L", "R")) {
+    for (d in c("32k", "164k")) {
+      spheres[[paste(h, d, "sphere")]] <- cached_download(sprintf("%s/tpl-fsLR_hemi-%s_den-%s_sphere.surf.gii", tf, h, d))
+      spheres[[paste(h, d, "area")]] <- cached_download(sprintf("%s/tpl-fsLR_hemi-%s_den-%s_desc-vaavg_midthickness.shape.gii", tf, h, d))
+    }
+  }
+  spec164 <- file.path(p$out_dir, "results_fsLR164k.spec")
+  if (file.exists(spec164)) unlink(spec164)
+  for (h in c("L", "R")) {
+    structure <- if (h == "L") "CORTEX_LEFT" else "CORTEX_RIGHT"
+    for (surf in c("inflated", "midthickness")) {
+      s32 <- cached_download(sprintf("%s/tpl-fsLR_den-32k_hemi-%s_%s.surf.gii", tf, h, surf))
+      s164 <- file.path(p$out_dir, sprintf("fsLR164k_hemi-%s_%s.surf.gii", h, surf))
+      if (!file.exists(s164)) {
+        wb_run("-surface-resample", s32, spheres[[paste(h, "32k", "sphere")]], spheres[[paste(h, "164k", "sphere")]], "BARYCENTRIC", s164)
+      }
+      wb_run("-add-to-spec-file", spec164, structure, s164)
+    }
+  }
+  wb_run("-add-to-spec-file", spec164, "INVALID", dscalar)
+
+  # A 32k copy for people who have their own 32k surfaces (parcel borders lose one vertex ring to NaN)
+  template32 <- file.path(p$cache_dir, "template_fsLR32k.dscalar.nii")
+  if (!file.exists(template32)) {
+    wb_run("-cifti-create-dense-scalar", template32, "-left-metric", spheres[["L 32k area"]], "-right-metric", spheres[["R 32k area"]])
+  }
+  dscalar32 <- file.path(p$out_dir, "results_fsLR32k.dscalar.nii")
+  wb_run("-cifti-resample", dscalar, "COLUMN", template32, "COLUMN", "ADAP_BARY_AREA", "CUBIC", dscalar32,
+         "-left-spheres", spheres[["L 164k sphere"]], spheres[["L 32k sphere"]],
+         "-left-area-metrics", spheres[["L 164k area"]], spheres[["L 32k area"]],
+         "-right-spheres", spheres[["R 164k sphere"]], spheres[["R 32k sphere"]],
+         "-right-area-metrics", spheres[["R 164k area"]], spheres[["R 32k area"]])
+  spec32 <- file.path(p$out_dir, "results_fsLR32k.spec")
+  if (file.exists(spec32)) unlink(spec32)
+  for (h in c("L", "R")) {
+    structure <- if (h == "L") "CORTEX_LEFT" else "CORTEX_RIGHT"
+    for (surf in c("inflated", "midthickness")) {
+      wb_run("-add-to-spec-file", spec32, structure, cached_download(sprintf("%s/tpl-fsLR_den-32k_hemi-%s_%s.surf.gii", tf, h, surf)))
+    }
+  }
+  wb_run("-add-to-spec-file", spec32, "INVALID", dscalar32)
+
+  cat("Open the results in Workbench with either of:\n\n")
+  cat("  wb_view", spec164, "\n")
+  cat("  wb_view", spec32, "\n\n")
+  cat("Pick the map 'parcel_median_late_fraction' and the inflated surfaces.",
+      "To render a saved scene to PNG without the GUI: wb_command -scene-capture-image <scene> 1 <out.png>\n")
+}
+```
+
+## 13. How long did it take, and how much was read?
+
+```{r runtime-table}
+runtime$secs <- round(runtime$secs, 1)
+runtime$mb_uncompressed <- round(runtime$mb_uncompressed)
+runtime$mb_wire <- round(runtime$mb_wire)
+runtime
+cat(sprintf("\nRead from S3: %.0f MB on the wire (%.0f MB uncompressed) in %s; whole report: %s.\n",
+            sum(runtime$mb_wire, na.rm = TRUE), sum(runtime$mb_uncompressed, na.rm = TRUE),
+            fmt_secs(sum(runtime$secs[!is.na(runtime$tiles)])), fmt_secs(as.numeric(Sys.time() - report_started, units = "secs"))))
+cat(sprintf("For comparison, the thickness scalar alone is %.1f GB; the whole store is 553 GB.\n", 340 * TILE_MB / 1024))
+utils::write.csv(runtime, file.path(p$out_dir, "runtime.csv"), row.names = FALSE)
+
+cat("\nSoftware: ModelArray", as.character(utils::packageVersion("ModelArray")),
+    "| tiledb", as.character(utils::packageVersion("tiledb")), "(libtiledb", paste(tiledb::tiledb_version(), collapse = "."), ")",
+    "| TileDBArray", as.character(utils::packageVersion("TileDBArray")),
+    "| mgcv", as.character(utils::packageVersion("mgcv")),
+    if (have_combat) paste("| ComBatFamily", as.character(utils::packageVersion("ComBatFamily"))) else "| ComBatFamily not installed",
+    if (have_wb) paste("| Workbench", sub("Version: *", "", grep("Version", system2(wb, "-version", stdout = TRUE), value = TRUE)[1])) else "", "\n")
+```
diff --git a/inst/examples/rbc-s3-tiledb/sa_axis_glasser360_ranks.csv b/inst/examples/rbc-s3-tiledb/sa_axis_glasser360_ranks.csv
new file mode 100644
index 0000000..78d379f
--- /dev/null
+++ b/inst/examples/rbc-s3-tiledb/sa_axis_glasser360_ranks.csv
@@ -0,0 +1,181 @@
+"region","brains.average.rank","final.rank"
+"V1",29.3,2
+"MST",65.6,39
+"V6",49.9,17.5
+"V2",30.4,3
+"V3",36.2,4
+"V4",48.05,14
+"V8",46.55,12
+"4",46.8,13
+"3b",44.5,7
+"FEF",86.8,85
+"PEF",104.1,117
+"55b",103.3,114.5
+"V3A",45.55,9
+"RSC",76.35,63
+"POS2",89.45,89
+"V7",58.8,26
+"IPS1",57.4,24
+"FFC",67.5,46
+"V3B",49,15
+"LO1",59.35,28
+"LO2",60.7,31
+"PIT",55.9,22
+"MT",69.2,52
+"A1",68.65,50
+"PSL",92.8,96
+"SFL",124.3,154.5
+"PCV",80.8,70
+"STV",85.7,82.5
+"7Pm",83,74
+"7m",114.6,136
+"POS1",91.8,94
+"23d",115.6,139
+"v23ab",119.8,149
+"d23ab",124.2,153
+"31pv",113.4,134
+"5m",42.3,5
+"5mv",51.8,19
+"23c",79.5,68
+"5L",60.6,30
+"24dd",62.7,33
+"24dv",75.95,62
+"7AL",69.9,54
+"SCEF",104.3,118
+"6ma",103.4,116
+"7Am",72.3,57
+"7PL",61.2,32
+"7PC",64.3,36
+"LIPv",81,72.5
+"VIP",69.3,53
+"MIP",73.2,59
+"1",49.05,16
+"2",58,25
+"3a",45.6,10
+"6d",46.4,11
+"6mp",56.8,23
+"6v",84.7,78
+"p24pr",99.8,106
+"33pr",86.2,84
+"a24pr",106.4,120
+"p32pr",104.5,119
+"a24",139.2,177
+"d32",129.2,166
+"8BM",130.2,167
+"p32",126.5,160
+"10r",115.1,137
+"47m",112.6,129
+"8Av",140.8,179
+"8Ad",138.9,176
+"9m",137.5,175
+"8BL",139.45,178
+"9p",132.4,170
+"10d",124.7,156
+"8C",133.4,172
+"44",115.3,138
+"45",124.3,154.5
+"47l",128.4,163.5
+"a47r",132.5,171
+"6r",101.3,109
+"IFJa",118.1,145
+"IFJp",110.3,124
+"IFSp",112.8,132
+"IFSa",109.7,123
+"p9-46v",127,161
+"46",135.65,174
+"a9-46v",133.85,173
+"9-46d",126.4,159
+"9a",123.9,152
+"10v",112.2,126
+"a10p",125.6,157
+"10pp",117.3,142
+"11l",118.3,146
+"13l",121.1,151
+"OFC",94.2,99
+"47s",118.5,147
+"LIPd",96.9,101
+"6a",88.5,88
+"i6-8",126.2,158
+"s6-8",128.4,163.5
+"43",78.65,66
+"OP4",66.55,41
+"OP1",67.25,44.5
+"OP2-3",54.95,21
+"52",79.35,67
+"RI",68.9,51
+"PFcm",87.55,86
+"PoI2",99.15,105
+"TA2",84.8,79
+"FOP4",98.85,104
+"MI",112.75,131
+"Pir",71.1,56
+"AVI",111,125
+"AAIC",97.3,102
+"FOP1",77.95,65
+"FOP3",92.65,95
+"FOP2",67.25,44.5
+"PFt",75.8,61
+"AIP",85.7,82.5
+"EC",81,72.5
+"PreS",54.5,20
+"H",80.45,69
+"ProS",59.1,27
+"PeEc",93.85,97
+"STGa",103.3,114.5
+"PBelt",65.8,40
+"A5",96.2,100
+"PHA1",68.25,49
+"PHA3",68.1,48
+"STSda",116.95,141
+"STSdp",100.5,108
+"STSvp",120.05,150
+"TGd",103.05,113
+"TE1a",132,169
+"TE1p",128.6,165
+"TE2a",112.5,128
+"TF",87.65,87
+"TE2p",98.3,103
+"PHT",103,112
+"PH",83.8,76.5
+"TPOJ1",89.8,90
+"TPOJ2",85.6,81
+"TPOJ3",67.8,47
+"DVT",67.2,43
+"PGp",83.7,75
+"IP2",102.9,111
+"IP1",85.4,80
+"IP0",72.55,58
+"PFop",83.8,76.5
+"PF",100.3,107
+"PFm",113.8,135
+"PGi",112.3,127
+"PGs",112.9,133
+"V6A",49.9,17.5
+"VMV1",29.1,1
+"VMV3",44.45,6
+"PHA2",73.5,60
+"V4t",64.1,35
+"FST",76.7,64
+"V3CD",63.65,34
+"LO3",66.85,42
+"VMV2",45.2,8
+"31pd",107.7,122
+"31a",102,110
+"VVC",65.1,37
+"25",112.7,130
+"s32",117.8,143
+"pOFC",90,91
+"PoI1",93.95,98
+"Ig",60.25,29
+"FOP5",106.6,121
+"p10p",131.5,168
+"p47r",127.7,162
+"TGv",91,92
+"MBelt",65.25,38
+"LBelt",70.4,55
+"A4",80.9,71
+"STSva",117.9,144
+"TE1m",144.6,180
+"PI",91.7,93
+"a32pr",116.8,140
+"p24",119,148
diff --git a/man/ModelArraySummary.Rd b/man/ModelArraySummary.Rd
index f1d073d..d958a42 100644
--- a/man/ModelArraySummary.Rd
+++ b/man/ModelArraySummary.Rd
@@ -17,7 +17,8 @@ h5summary(filepath, backend = c("auto", "hdf5", "tiledb"))
 }
 \arguments{
 \item{filepath}{Character. Path to an HDF5 (\code{.h5}) file or TileDB
-(\code{.tdb}) store.}
+(\code{.tdb}) store. A TileDB store may be given as a URI
+(\code{s3://...}, \code{file://...}); HDF5 files must be local paths.}
 
 \item{backend}{Character. Storage backend: \code{"auto"} (default),
 \code{"hdf5"}, or \code{"tiledb"}.}
@@ -55,6 +56,21 @@ For HDF5, this function opens the file read-only via
 into memory. The returned
 object has a \code{print} method that displays a formatted summary.
 
+A TileDB store can also be remote, for example in an S3 bucket. Configure
+the TileDB context before calling, exactly as you would before
+\code{\link{ModelArray}()}:
+\preformatted{
+cfg <- tiledb::tiledb_config()
+cfg["vfs.s3.region"] <- "us-east-1"
+cfg["vfs.s3.no_sign_request"] <- "true"   # public bucket, no credentials
+tiledb::tiledb_ctx(cfg)
+ModelArraySummary("s3://bucket/path/store.tdb")
+}
+The same settings can come from environment variables
+(\code{TILEDB_VFS_S3_REGION}, \code{TILEDB_VFS_S3_NO_SIGN_REQUEST}).
+Each scalar costs a few requests, so summarizing a store with many scalars
+over the network takes seconds rather than milliseconds.
+
 \code{h5summary()} is a backward-compatible alias for
 \code{ModelArraySummary()}. Despite the historical name, it supports both
 HDF5 files and TileDB stores.
@@ -63,6 +79,7 @@ HDF5 files and TileDB stores.
 \dontrun{
 ModelArraySummary("path/to/data.h5")
 ModelArraySummary("path/to/store.tdb")
+ModelArraySummary("s3://bucket/path/store.tdb") # after configuring the TileDB context
 
 # Inspect before deciding which scalars to load
 info <- ModelArraySummary("path/to/data.h5")
diff --git a/tests/testthat/test-convenience_accessors.R b/tests/testthat/test-convenience_accessors.R
index 7871710..805a366 100644
--- a/tests/testthat/test-convenience_accessors.R
+++ b/tests/testthat/test-convenience_accessors.R
@@ -75,3 +75,8 @@ test_that("h5summary errors on missing file", {
 test_that("ModelArraySummary errors on missing file", {
   expect_error(ModelArraySummary("/nonexistent/file.h5"), "not found")
 })
+
+test_that("ModelArraySummary still needs a local path for HDF5", {
+  h5_path <- system.file("extdata", "n50_fixels.h5", package = "ModelArray")
+  expect_error(ModelArraySummary(paste0("file://", normalizePath(h5_path))), "not found")
+})
diff --git a/tests/testthat/test-tiledb_backend.R b/tests/testthat/test-tiledb_backend.R
index 7b5b6d5..0f0c58b 100644
--- a/tests/testthat/test-tiledb_backend.R
+++ b/tests/testthat/test-tiledb_backend.R
@@ -98,6 +98,68 @@ test_that("ModelArraySummary summarizes TileDB scalar arrays", {
   expect_equal(legacy$backend, "tiledb")
 })
 
+test_that("ModelArraySummary reads a TileDB store through a file:// URI", {
+  skip_if_no_tiledb()
+
+  store <- tempfile(fileext = ".tdb")
+  on.exit(unlink(store, recursive = TRUE), add = TRUE)
+  write_tiledb_scalar(store, "FA", matrix(seq_len(12), nrow = 3), paste0("sub-", 1:3))
+  writeResults(store, data.frame(element_id = 0:3, estimate = 1:4), analysis_name = "lm_uri", backend = "tiledb")
+
+  # file.exists() is FALSE for a URI, which is what a remote (s3://) store looks like to R
+  uri <- paste0("file://", normalizePath(store))
+  expect_false(file.exists(uri))
+
+  info <- ModelArraySummary(uri)
+  expect_equal(info$backend, "tiledb")
+  expect_identical(info$filepath, uri)
+  expect_equal(info$scalars$name, "FA")
+  expect_equal(info$scalars$nElements, 4L)
+  expect_equal(info$scalars$nInputFiles, 3L)
+  expect_equal(info$analyses, "lm_uri")
+
+  local_info <- ModelArraySummary(store)
+  expect_equal(info[c("scalars", "analyses")], local_info[c("scalars", "analyses")])
+})
+
+test_that("ModelArraySummary ignores group internals and stray files in a TileDB store", {
+  skip_if_no_tiledb()
+
+  store <- tempfile(fileext = ".tdb")
+  on.exit(unlink(store, recursive = TRUE), add = TRUE)
+  write_tiledb_scalar(store, "FD", matrix(1:4, nrow = 2), c("s1", "s2"), metadata_location = "group")
+  # Make scalars/ a real TileDB group (adds __tiledb_group.tdb, __group, __meta) and drop stray files
+  write_tiledb_group_metadata(file.path(store, "scalars"), "note", "x")
+  dir.create(file.path(store, "results"))
+  writeLines("x", file.path(store, "scalars", "README.txt"))
+  writeLines("x", file.path(store, "results", "README.txt"))
+
+  info <- ModelArraySummary(paste0("file://", normalizePath(store)))
+  expect_equal(info$scalars$name, "FD")
+  expect_equal(info$analyses, character(0))
+})
+
+test_that("ModelArraySummary errors clearly when a TileDB store is missing", {
+  skip_if_no_tiledb()
+
+  expect_error(ModelArraySummary("file:///nonexistent/store.tdb"), "not found")
+  expect_error(ModelArraySummary(tempfile(fileext = ".tdb")), "not found")
+})
+
+test_that("ModelArray constructs from a file:// TileDB URI", {
+  skip_if_no_tiledb()
+
+  store <- tempfile(fileext = ".tdb")
+  on.exit(unlink(store, recursive = TRUE), add = TRUE)
+  values <- matrix(seq_len(8), nrow = 2, byrow = TRUE)
+  write_tiledb_scalar(store, "FD", values, c("sub-1.nii.gz", "sub-2.nii.gz"))
+
+  ma <- ModelArray(paste0("file://", normalizePath(store)), scalar_types = "FD")
+  loaded <- as.matrix(scalars(ma)[["FD"]])
+  dimnames(loaded) <- NULL
+  expect_equal(loaded, t(values))
+})
+
 test_that("ModelArraySummary agrees for fake HDF5 and TileDB stores", {
   skip_if_no_tiledb()
 
diff --git a/vignettes/installations.Rmd b/vignettes/installations.Rmd
index 71bbf28..6ac5f71 100644
--- a/vignettes/installations.Rmd
+++ b/vignettes/installations.Rmd
@@ -126,6 +126,8 @@ library(ModelArray)
 ModelArraySummary("path/to/store.tdb", backend = "tiledb")
 ```
 
+`ModelArraySummary()` and `ModelArray()` also take a remote TileDB URI such as `s3://bucket/store.tdb`, once the TileDB context knows how to reach it. `vignette("rbc-s3-tiledb")` runs a full analysis on a public S3 store that way.
+
 If installing the `tiledb` R package fails, install the TileDB system library first. On macOS:
 
 ``` {.console}
diff --git a/vignettes/modelling.Rmd b/vignettes/modelling.Rmd
index cc3904f..4cc77d7 100644
--- a/vignettes/modelling.Rmd
+++ b/vignettes/modelling.Rmd
@@ -275,14 +275,18 @@ head(result)
 writeResults(h5_path, df.output = result, analysis_name = "site_analysis")
 ```
 
-### Case study: harmonize with `covfam` and stream to a new scalar
+### Case study: harmonize with ComBat and stream to a new scalar
 
 You can also use `ModelArray.wrap()` as a transformation engine rather than a test runner.
 The example below applies a harmonization step per element, streams the harmonized values
 directly into `/scalars`, and then loads that harmonized scalar in a new `ModelArray`.
+It uses [ComBatFamily](https://github.com/andy1764/ComBatFamily)
+(`remotes::install_github("andy1764/ComBatFamily")`), whose `comfam()` is ComBat and
+`covfam()` is CovBat; both take the same arguments. A worked example on real data,
+with a GAM age term and CovBat, is in `vignette("rbc-s3-tiledb")`.
 
 ```{r wrap_covfam_case_study, eval=FALSE}
-library(covfam)
+library(ComBatFamily)
 
 # Assume phenotypes has harmonization variables, e.g.:
 #   source_file, site, age, sex
@@ -297,44 +301,48 @@ phenotypes <- read.csv("phenotypes.csv")
 # Copy metadata/layout into a new file so we can append new scalars there.
 file.copy(h5_in, h5_harmonized, overwrite = TRUE)
 
-covfam_harmonize_element <- function(data) {
-  # Harmonize one element across subjects.
-  # Adjust arguments here to match your covfam configuration.
-  out <- covfam::covfam(
-    y = data$thickness,
-    batch = data$site,
-    mod = data.frame(age = data$age, sex = data$sex)
+combat_harmonize_element <- function(data) {
+  # Harmonize one element across subjects, keeping the age and sex effects.
+  # ModelArray hands the function one element (one feature) at a time, so the
+  # empirical-Bayes step has nothing to pool across: eb = FALSE.
+  out <- ComBatFamily::comfam(
+    data = data[, "thickness", drop = FALSE],
+    bat = factor(data$site),
+    covar = data[, c("age", "sex")],
+    model = stats::lm,
+    formula = y ~ age + sex,
+    eb = FALSE
   )
 
   # Return named subject-level vector so wrap columns map to source_file order.
-  vals <- out$y_harmonized
+  vals <- out$dat.combat[, "thickness"]
   names(vals) <- data$source_file
   vals
 }
 
-# Stream harmonized subject-level values into scalars/thickness_covfam/values.
+# Stream harmonized subject-level values into scalars/thickness_combat/values.
 # Set return_output = FALSE to avoid keeping the full output table in memory.
 ModelArray.wrap(
-  FUN = covfam_harmonize_element,
+  FUN = combat_harmonize_element,
   data = ma_raw,
   phenotypes = phenotypes,
   scalar = "thickness",
   n_cores = 8,
-  write_scalar_name = "thickness_covfam",
+  write_scalar_name = "thickness_combat",
   write_scalar_file = h5_harmonized,
   write_scalar_flush_every = 2000L,
   return_output = FALSE
 )
 
 # Load the harmonized scalar as a new modelling input
-ma_harmonized <- ModelArray(h5_harmonized, scalar_types = c("thickness_covfam"))
+ma_harmonized <- ModelArray(h5_harmonized, scalar_types = c("thickness_combat"))
 
 # Use harmonized data in downstream models
 fit_harmonized <- ModelArray.lm(
-  thickness_covfam ~ age + sex,
+  thickness_combat ~ age + sex,
   data = ma_harmonized,
   phenotypes = phenotypes,
-  scalar = "thickness_covfam",
+  scalar = "thickness_combat",
   n_cores = 8
 )
 ```
diff --git a/vignettes/rbc-s3-tiledb.Rmd b/vignettes/rbc-s3-tiledb.Rmd
new file mode 100644
index 0000000..53b9a39
--- /dev/null
+++ b/vignettes/rbc-s3-tiledb.Rmd
@@ -0,0 +1,238 @@
+---
+title: "Cloud TileDB: RBC Cortical Thickness from S3"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Cloud TileDB: RBC Cortical Thickness from S3}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+[![Launch in Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/PennLINC/ModelArray/tiledb-example?urlpath=rstudio)
+
+The [Reproducible Brain Charts](https://reprobrainchart.github.io/) (RBC) FreeSurferPost
+release is available as a TileDB store on a public S3 bucket: 70 surface measures,
+each 327,684 fsLR-164k vertices by 8,696 sessions from five studies (BHRC, CCNP, HBN,
+NKI, PNC; ages 5–85). That is 553 GB. This article is about a report you can run
+yourself that fits models to those data **without downloading them**: ModelArray reads
+the handful of tiles it needs, straight from S3, and leaves the rest where it is.
+
+The report asks a question with a known answer, so you can tell whether it worked.
+The cortex thins through childhood and adolescence, but not everywhere at the same
+time: primary sensory cortex matures first and has little thinning left after
+mid-adolescence, while association cortex keeps thinning through the teens and into
+the twenties (Shaw et al., 2008; Tamnes et al., 2017). Sydnor et al. (2021) order every
+region along a sensorimotor–association (S-A) axis, from sensory and motor cortex at one
+end to prefrontal association cortex at the other. Put together: **a parcel near the
+sensorimotor end should do most of its thinning early, and one near the association end
+should do most of it late.** The report tests this on two Glasser parcels at opposite
+ends of the axis — area 3b (primary somatosensory, rank 7) and area 46 (dorsolateral
+prefrontal, rank 174) — one vertex at a time, and never on the whole cortex.
+
+## What the report does
+
+1. Points TileDB at the bucket (two config settings, no AWS account).
+2. Lists the store with `ModelArraySummary()`, opens `thickness` with `ModelArray()`,
+   and downloads the cohort CSV (age, sex, site).
+3. Works out, before reading anything, how many tiles each parcel touches and how
+   many megabytes that is.
+4. Looks at single vertices with `exampleElementData()`: thickness against age, a GAM
+   with site as a covariate, its derivative, and a mixed model for repeat sessions.
+5. Reads each parcel's vertices from S3 in one go and harmonizes them across the nine
+   scanning sites with CovBat (`ComBatFamily::covfam()`, keeping `s(Age, k = 4) + sex`),
+   then writes the raw and harmonized blocks into a small local ModelArray store.
+6. Runs `ModelArray.gam()` straight from S3 on one parcel (site as a covariate) and on
+   the local store for every parcel, raw and harmonized, to show what harmonization
+   changes; then `ModelArray.wrap()` with a custom `fit_vertex()` that adds an
+   `fx = TRUE` test statistic, ΔR²adj, predicted thickness every two years and the
+   derivative of the age curve with a simultaneous interval, for ages 5–23.
+7. Summarizes per parcel, plots the late fraction of thinning and the slope-by-age
+   curves against S-A rank, and compares them with the literature.
+8. Writes CSVs, saves the results into the local store with `writeResults()`, and
+   exports CIFTI dscalars plus `.spec` files for Connectome Workbench.
+9. Prints how long each step took and how many bytes came over the wire.
+
+## Run it
+
+### In the browser
+
+Click the Binder badge at the top. After a minute or two RStudio opens with
+`~/rbc_s3_example/rbc_s3_tiledb_example.Rmd` ready; open it and press **Knit**. Binder
+gives one core and about 2 GB of memory, so the report detects it and uses its smallest
+settings (three left-hemisphere parcels, one tile each). Expect 10–15 minutes. Binder
+sessions are public and disappear after ten idle minutes, so download anything you
+want to keep (the `rbc_s3_example` folder has the HTML, CSVs and CIFTI files).
+
+### On your own machine
+
+You need ModelArray from the GitHub main branch (TileDB support), the TileDB packages,
+and a few modelling packages:
+
+```{r install, eval=FALSE}
+install.packages(c("tiledb", "jsonlite", "mgcv", "ggplot2", "gratia", "gamm4", "lme4", "rmarkdown", "cowplot", "png"))
+BiocManager::install("TileDBArray")
+remotes::install_github("andy1764/ComBatFamily")
+remotes::install_github("PennLINC/ModelArray")
+```
+
+[Connectome Workbench](https://www.humanconnectome.org/software/get-connectome-workbench)
+is optional; with `wb_command` on your PATH the report writes the CIFTI files and a
+`.spec` you can open in `wb_view`. Then render:
+
+```{r render, eval=FALSE}
+example <- system.file("examples/rbc-s3-tiledb/rbc_s3_tiledb_example.Rmd", package = "ModelArray")
+out_dir <- path.expand("~/rbc_s3_example")
+dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
+
+# quick: areas 3b and 46 in the left hemisphere, two tiles per parcel
+rmarkdown::render(example, output_dir = out_dir, knit_root_dir = out_dir, intermediates_dir = tempdir(),
+                  params = list(quick = TRUE, n_cores = 4, out_dir = out_dir))
+
+# full: areas 3b and 46 in both hemispheres, six tiles per parcel
+rmarkdown::render(example, output_dir = out_dir, knit_root_dir = out_dir, intermediates_dir = tempdir(),
+                  params = list(quick = FALSE, n_cores = 8, out_dir = out_dir))
+```
+
+Other parameters (`parcels`, `hemispheres`, `max_tiles_per_parcel`,
+`max_vertices_per_parcel`, `harmonize`, `write_cifti`, `allocation_mb`) are listed at
+the top of the Rmd. The report never writes to S3; everything lands in `out_dir`.
+
+## What to expect
+
+Numbers from runs on an Apple-silicon Mac (14 cores, 36 GB; the report used 8 cores).
+
+| | Quick | Full |
+|---|---|---|
+| Parcels | 3b, 46 (left) | 3b, 46 (both hemispheres) |
+| Vertices analysed | 600 | 2,000 |
+| Read from S3 | 296 MB (2 tiles per parcel) | 996 MB (6 tiles per parcel) |
+| CovBat | 2.5 min on 2 cores | 4.1 min on 4 cores |
+| `fit_vertex()` at every vertex | 42 s | 3.6 min |
+| Whole report | 8 min | 20 min |
+
+The full thickness scalar is 10.6 GB and the whole store 553 GB, so the full run
+reads about 0.2% of the store.
+
+Every tile check passed: each parcel read exactly the tiles predicted before
+reading. The derivative of the age curve matches `gratia::derivatives()` to about
+1e-9 mm/yr.
+
+The literature check passes. Area 3b (sensorimotor) does most of its thinning early
+and area 46 (association) does most of its late: the share of the 8–22 thinning that
+happens after age 14 is about 0.4 in 3b and about 0.7 in 46 (Spearman ρ = 0.89 across
+the four parcel-hemispheres; age of fastest thinning ρ = 1.00):
+
+| Parcel | S-A rank | Late fraction (L / R) | Age of fastest thinning | Slope no longer clearly negative after |
+|---|---|---|---|---|
+| 3b | 7 | 0.41 / 0.48 | 5 y | 19–23 y |
+| 46 | 174 | 0.74 / 0.71 | 23 y | still thinning at 23 |
+
+Area 3b thins fastest in early childhood and is nearly flat by the late teens; area 46
+thins slowly in childhood and fastest in the late teens and twenties — exactly the
+sensorimotor-first, association-last ordering of the S-A axis.
+
+The S-A axis and the two parcels, in the colour scheme of Sydnor et al. (2021,
+Figure 2A) — gold at the sensorimotor pole, purple at the association pole:
+
+```{r fig-sa, echo=FALSE, out.width="70%", fig.align="center"}
+knitr::include_graphics("rbc_s3_brain_sa_axis.png")
+```
+
+```{r fig-vertex, echo=FALSE, out.width="100%"}
+knitr::include_graphics("rbc_s3_fig1_vertex_thickness_vs_age.png")
+```
+
+```{r fig-sites, echo=FALSE, out.width="100%"}
+knitr::include_graphics("rbc_s3_fig_site_boxplots.png")
+```
+
+```{r fig-rate, echo=FALSE, out.width="100%"}
+knitr::include_graphics("rbc_s3_fig2_late_fraction_vs_sa_rank.png")
+```
+
+```{r fig-slopes, echo=FALSE, out.width="100%"}
+knitr::include_graphics("rbc_s3_fig4_parcel_slopes.png")
+```
+
+
+## How it works
+
+**Reaching the bucket.** ModelArray uses the `tiledb` package's default context. Set
+the region and anonymous access on it once, before `ModelArray()`:
+
+```{r ctx, eval=FALSE}
+cfg <- tiledb::tiledb_config()
+cfg["vfs.s3.region"] <- "us-east-1"
+cfg["vfs.s3.no_sign_request"] <- "true"
+tiledb::tiledb_ctx(cfg)
+
+uri <- "s3://fcp-indi/data/Projects/RBC/RBC_FreeSurfer-Post_ModelArray/rbc_freesurferpost.tdb"
+ModelArraySummary(uri)
+ma <- ModelArray(uri, scalar_types = c("thickness", "curv"))
+```
+
+The environment variables `TILEDB_VFS_S3_REGION` and `TILEDB_VFS_S3_NO_SIGN_REQUEST`
+do the same thing. For a private bucket, use your usual AWS credentials instead of
+`no_sign_request`.
+
+**Tiles are the unit of cost.** The store is dense float32, tiled as 964 vertices by
+all 8,696 sessions (33.5 MB before compression, 4–8 MB on the wire). Reading one vertex
+reads its tile, so `element.subset` ranges that fall in the same tile are nearly free,
+and a parcel whose vertices are spread over 20 tiles costs 20 tiles. The report prints
+this arithmetic before it reads, and caps each parcel at its most populated tiles.
+
+**One read pass per query.** The `tiledb` R package reads a query in passes whose size
+is set by `tiledb::set_allocation_size_preference()` (10 MiB by default). Each pass
+re-fetches the tiles it touches, so with the default a 3-tile read can download 18
+tiles. The report sets 512 MB; the report's runtime table checks that the number of
+tiles read matches the prediction. This is separate from
+`options(ModelArray.tiledb_read_block_mb)`, which controls how many elements ModelArray
+asks for at once (see `vignette("large-scale-analyses")`).
+
+**Results stay local.** `writeResults()` writes to a local `.h5` (or `.tdb`); the S3
+store is never written to.
+
+**Why CovBat on a parcel at a time.** Nine sites differ in both the mean and the
+spread of thickness, so a `site` covariate is not enough. CovBat harmonizes mean,
+variance and the covariance between features across sites. Harmonizing vertex by
+vertex inside the per-element function would work, but it is slow (one GAM fit per
+vertex just for the correction) and leaves nothing for the empirical-Bayes step to pool
+across. Reading a parcel's vertices as one block and harmonizing them together is the
+method's intended use, costs one `covfam()` call per parcel, and the harmonized block
+goes into a local HDF5 store that ModelArray reads like any other. Parcels are
+harmonized independently, so the report runs them in parallel, one per core
+(`parallel::mclapply()`; serially on Windows and on Binder). CovBat keeps every
+per-vertex GAM fit, so the report drops the copies of the data each fit carries,
+which `covfam()` never uses again; the harmonized values are unchanged, and each
+worker needs about 1 GB per 300 vertices. Repeat sessions
+(BHRC, NKI) are treated as independent, which the report says plainly and contrasts
+with a one-vertex longitudinal ComBat + GAMM fit.
+
+**One scalar at a time.** ModelArray can attach a second scalar as an element-wise
+covariate when it lists the same source files as the phenotypes. In this store every
+scalar has its own file per session, so `curv` cannot ride along with `thickness`; the
+report uses thickness alone.
+
+**GAM choices.** `s(Age, k = 4)` with REML, following PennLINC practice: a developmental
+curve needs at most two or three bends. Penalized p-values can be optimistic, so the
+test statistic comes from an `fx = TRUE` refit. Effect size is the change in adjusted
+R² when the smooth is dropped (what `changed.rsq.term.index` gives in
+`ModelArray.gam()`). Direction and timing come from the derivative of the fitted curve
+with a simultaneous interval (the same construction as `gratia::derivatives()`).
+
+## See also
+
+- `vignette("large-scale-analyses")` for memory and block-size settings on big stores.
+- `vignette("modelling")` for `ModelArray.gam()`, `ModelArray.wrap()` and the
+  streaming-harmonization case study.
+- `vignette("modelarray-storage")` for the HDF5 and TileDB layouts.
+- The example folder: `system.file("examples/rbc-s3-tiledb", package = "ModelArray")`,
+  including the provenance of the vendored Glasser-164k labels and S-A ranks.
+
+## References
+
+- Bethlehem RAI et al. (2022). Brain charts for the human lifespan. *Nature* 604, 525–533.
+- Glasser MF et al. (2016). A multi-modal parcellation of human cerebral cortex. *Nature* 536, 171–178.
+- Sydnor VJ et al. (2021). Neurodevelopment of the association cortices. *Neuron* 109, 2820–2846.
+- Shaw P et al. (2008). Neurodevelopmental trajectories of the human cerebral cortex. *J Neurosci* 28, 3586–3594.
+- Tamnes CK et al. (2017). Development of the cerebral cortex across adolescence. *J Neurosci* 37, 3402–3412.
+- Chen AA et al. (2022). Mitigating site effects in covariance for machine learning in neuroimaging data. *Hum Brain Mapp* 43, 1179–1195. (CovBat)
diff --git a/vignettes/rbc_s3_brain_sa_axis.png b/vignettes/rbc_s3_brain_sa_axis.png
new file mode 100644
index 0000000..e1db676
Binary files /dev/null and b/vignettes/rbc_s3_brain_sa_axis.png differ
diff --git a/vignettes/rbc_s3_fig1_vertex_thickness_vs_age.png b/vignettes/rbc_s3_fig1_vertex_thickness_vs_age.png
new file mode 100644
index 0000000..527f4b3
Binary files /dev/null and b/vignettes/rbc_s3_fig1_vertex_thickness_vs_age.png differ
diff --git a/vignettes/rbc_s3_fig2_late_fraction_vs_sa_rank.png b/vignettes/rbc_s3_fig2_late_fraction_vs_sa_rank.png
new file mode 100644
index 0000000..b5d94ba
Binary files /dev/null and b/vignettes/rbc_s3_fig2_late_fraction_vs_sa_rank.png differ
diff --git a/vignettes/rbc_s3_fig4_parcel_slopes.png b/vignettes/rbc_s3_fig4_parcel_slopes.png
new file mode 100644
index 0000000..47e5d2c
Binary files /dev/null and b/vignettes/rbc_s3_fig4_parcel_slopes.png differ
diff --git a/vignettes/rbc_s3_fig_site_boxplots.png b/vignettes/rbc_s3_fig_site_boxplots.png
new file mode 100644
index 0000000..14dabfa
Binary files /dev/null and b/vignettes/rbc_s3_fig_site_boxplots.png differ