Matched section-wise 2D networks and a 3D comparison

This vignette uses a small registered synthetic cell map to run the matched section-wise planar comparator. The 2D model fits a Matérn-3/2 GP independently within each zone and section using $(x,y)$ coordinates; its adjusted rows are then passed to the same Gaussian-likelihood shared-plus-zone covariance fitter as the 3D model. The matched comparison keeps the selected cell roster and downstream estimator fixed. It changes both coordinate dimension and the way sections are pooled, so a 3D-minus-2D difference does not isolate a causal effect of depth.

The planar workflow follows the cell-density, spatial-adjustment, factor-covariance, and conditional-network architecture of the published ISPat method (Bhadury et al., 2026). Its scalable Vecchia and Gaussian-likelihood estimators are those used for the ISPAT-3D comparison described by Bhadury and Rao (2026); it is not an execution of the original Bayesian ISPat package.

1. Prepare a small serial-section image

library(ISPAT3D)
image <- ispat3d_example_image(n_per_section = 30, n_sections = 3,
                               bandwidth = 0.12, seed = 2026)
table(image$zones, image$sections)

##       
##         1  2  3
##   Low  15 15 15
##   High 15 15 15

The example has three registered sections, three annotated source cell types, a cell-type-specific KDE at every cell, and Low/High zones defined by each section’s tumor-cell KDE median. The model matrix image$Y is log1p(1e9 * image$kde). This miniature image is for execution and teaching. Real images need checked cell classification, registration, section IDs, KDE surfaces, and a prespecified tumor-burden score before fitting.

The section-zone group sizes matter because the planar GP is fitted inside each group. A group with fewer than ten cells is mean-centered. With the selection below, each group has 12 cells and exercises the actual GP code.

old <- par(mfrow = c(1, 3), mar = c(3, 3, 2, 1))
zone_color <- c(Low = "#4C72B0", High = "#C44E52")
for (s in sort(unique(image$sections))) {
  at <- image$sections == s
  plot(image$coords[at, 1:2], xlim = c(0, 1), ylim = c(0, 1),
       xlab = "x", ylab = "y", main = paste("Section", s), pch = 19,
       cex = 0.6, col = zone_color[as.character(image$zones[at])])
}

Synthetic planar cell locations, colored by relative tumor-density zone.

par(old)

2. Fix the cell roster before either fit

Call ispat3d_sample() once and pass its row indices to both analyses. This is the key matching rule; selecting separate random subsets would conflate sampling with the spatial-model comparison.

selected <- ispat3d_sample(image$coords, image$zones, image$sections,
                           budget = 36L, budget_kind = "per_zone",
                           seed = 2027L)
rows <- unlist(selected, use.names = FALSE)
table(image$zones[rows], image$sections[rows])

##       
##         1  2  3
##   Low  19 13  1
##   High 29  6  4

stopifnot(length(rows) == 72L)
Y <- image$Y[rows, , drop = FALSE]
xyz <- image$coords[rows, , drop = FALSE]
zone <- image$zones[rows]
section <- image$sections[rows]

For experimental CRC-style analysis, the budget denotes cells per zone. For the manuscript’s breast specimen, budget_kind = "total" allocates one budget across the five zones. Save selected with saveRDS() when comparing fit variants or repeating an analysis.

3. Fit the section-wise planar model

ispat3d_fit_2d() takes the same three-column registered coordinate matrix as the 3D function but uses only the first two columns for each section-zone GP. GP predictions and residuals are made for all selected cells in that group. The complete adjusted covariance for each zone then enters the same factor-covariance likelihood as the volumetric fit.

This example uses a small rank and short optimizer budget for a fast executable vignette. Defaults for large applications are 15 Vecchia neighbors and up to 5,000 anchors per fit; the manuscript used rank five for CRC and breast.

common <- list(
  Y = Y, coords = xyz, zones = zone, sections = section,
  rank = 2L, anchor_min = 12L, anchor_max = 24L,
  neighbors = 5L, gp_maxit = 3L, factor_maxit = 60L,
  threads = 1L, return_residuals = TRUE
)
fit2d <- do.call(ispat3d_fit_2d, common)
table(fit2d$gp_log$status)

## 
##            ok small_section 
##             9             9

fit2d$counts

##  Low High 
##   33   39

Inspect fit2d$gp_log before using the result. A small_section or constant status means mean-centering was used as specified; fallback:... means an attempted GP failed and that particular variable-group was mean-centered. The returned residual matrices are available because return_residuals = TRUE was set.

4. Compute and plot conditional networks

The fitted covariance in zone $q$ is the sum of shared factors, zone factors, and positive diagonal uniqueness. The partial correlation for a pair is obtained by inverting this full covariance and normalizing its precision entries. The shared factor product alone is low rank and should not be inverted as a separate network.

round(fit2d$full[["High"]], 3)

##             Tumor T_cell Macrophage
## Tumor       0.058 -0.038     -0.033
## T_cell     -0.038  0.487      0.033
## Macrophage -0.033  0.033      0.376

round(fit2d$partial[["High"]], 3)

##             Tumor T_cell Macrophage
## Tumor       1.000 -0.214     -0.209
## T_cell     -0.214  1.000      0.029
## Macrophage -0.209  0.029      1.000

stopifnot(isTRUE(all.equal(
  ispat3d_partial_correlation(fit2d$full[["High"]]),
  fit2d$partial[["High"]]
)))
ispat3d_edge_table(fit2d, zone = "High", threshold = 0.02)

##     from         to partial_correlation     sign absolute_effect
## 1  Tumor     T_cell         -0.21412730 negative      0.21412730
## 2  Tumor Macrophage         -0.20938917 negative      0.20938917
## 3 T_cell Macrophage          0.02929854 positive      0.02929854

ispat3d_plot_zones(fit2d, columns = 2, threshold = 0.02,
                   label_cex = 0.8)

Matched planar conditional-density networks in the two simulated zones.

The edge table and base-R circular plots use the fitted partial correlations. Red and blue edges are positive and negative conditional density associations. A display threshold declutters the drawing; it is not a test of significance.

5. Compare with the volumetric fit on identical cells

fit3d <- do.call(ispat3d_fit, common)
stopifnot(identical(fit3d$counts, fit2d$counts))
delta_high <- fit3d$partial[["High"]] - fit2d$partial[["High"]]
round(delta_high, 3)

##             Tumor T_cell Macrophage
## Tumor       0.000 -0.133     -0.150
## T_cell     -0.133  0.000     -0.133
## Macrophage -0.150 -0.133      0.000

old <- par(mfrow = c(1, 2), mar = c(1, 1, 3, 1))
ispat3d_plot_network(fit3d, zone = "High", threshold = 0.02,
                     main = "Pooled 3D", label_cex = 0.8)
ispat3d_plot_network(fit2d, zone = "High", threshold = 0.02,
                     main = "Section-wise 2D", label_cex = 0.8)

The same selected cells analyzed with pooled 3D and section-wise planar GP adjustment.

par(old)

A 3D-minus-2D contrast is a difference between two complete spatial-adjustment strategies. It should be reported with effect magnitudes, GP fit diagnostics, sampling details, and the caution that a smooth field can contain meaningful tissue biology. Neither fit identifies physical contacts or cell-cell signaling from density measurements alone.

6. Use processed serial-section data

# dat contains registered X,Y,Z, section, zone and cell-type KDE columns.
coords <- as.matrix(dat[, c("X", "Y", "Z")])
kde <- as.matrix(dat[, c("kde_Tumor", "kde_T_cell", "kde_Macrophage")])
colnames(kde) <- c("Tumor", "T_cell", "Macrophage")
Y <- log1p(1e9 * kde)
selected <- ispat3d_sample(coords, dat$zone, dat$section,
                           budget = 1000L, budget_kind = "per_zone")
rows <- unlist(selected, use.names = FALSE)
args <- list(Y = Y[rows, , drop = FALSE],
             coords = coords[rows, , drop = FALSE],
             zones = dat$zone[rows], sections = dat$section[rows],
             rank = 2L)
planar <- do.call(ispat3d_fit_2d, args)
volumetric <- do.call(ispat3d_fit, args)
ispat3d_plot_zones(planar, threshold = 0.05)

See the README for the exact processed CRC and breast input schemas and full sample budgets.

References

Bhadury, S., Peruzzi, M., Acharyya, S., et al. (2026). “Informed spatially aware patterns for multiplexed immunofluorescence data.” Scientific Reports 16, 5015. doi:10.1038/s41598-026-35341-8.

Bhadury, S. and Rao, A. (2026). “Estimating Conditional Cell Population Associations Across Tumor Density Zones in Three Dimensional Tissue Images.” Manuscript submitted to Scientific Reports. Accompanying software archive: doi:10.5281/zenodo.22798429. This is the software DOI, not a publication DOI for the manuscript.