Contents

1 Introduction

This vignette demonstrates how to use spatialFDA to perform functional data analysis on spatial statistics metrics. The main aim of this package is to detect differential spatial arrangements within and between celltypes given several samples/conditions. It does so by calculating spatial statistics metrics via the spatstat package and comparing differences using functional additive mixed models as implemented in the refund package (Baddeley and Turner 2005; Goldsmith et al. 2024).

The use case is a dataset from the imcdatasets package, which contains images from 12 human donors (Damond et al. 2019).

This package is similar to other packages in python and R. The following table shows the main differences in terms of functionality (Ali et al. 2024; Canete et al. 2022; Wrobel et al. 2024).

Package name Foundation Testing procedure
spicyR \(L\)-function Scalar comparison
GraphCompass Graph-based Graph and scalar comparison
mxfda K-, G- and L- function Function as input for survival modelling
spatialFDA most spatstat functions Functional comparison over domain

2 Installation

spatialFDA can be installed and loaded from Bioconductor as follows

if (!requireNamespace("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install("spatialFDA")
library("spatialFDA")
library("dplyr")
library("ggplot2")
library("tidyr")
library("stringr")
library("dplyr")
library("patchwork")

3 Getting started

In this vignette we will analyse a diabetes dataset acquired by imaging mass cytometry (IMC) as developed by Damond et al. (Damond et al. 2019). The dataset contains images from 12 human donors, 4 healthy and 8 with type 1 diabetes (T1D). With IMC, 35 markers were measured at single cell resolution (Damond et al. 2019).

3.1 Loading the data

The Damond et al. (Damond et al. 2019) dataset is easily loaded via the imcdatasets package. The entire dataset can be loaded by setting full_dataset = TRUE. For computational reasons, one can reduce to three patients as well by setting this flag to FALSE. We will use the small version to show the functionality of the package even though the statistical power of results is greatly reduced by considering only 100 out of 845 images. For a more interesting use case we encourage readers to load the full dataset for a multi-condition/multi-sample example. The package offers multiple datatypes, we will use the SpatialExperiment (SPE) object.

# load the dataset as SpatialExperiment object
spe <- imcdatasets::Damond_2019_Pancreas("spe", full_dataset = TRUE)
# subset the example for computational reasons
spe <- subset(spe, ,patient_id %in% c(6089,6180,6126,6134,6228,6362))

3.2 Visualising the raw data

We can look at the fields of view (FOVs) of the diabetes dataset. To do so we extract the spatial coordinates, store them as a dataframe and add the colData from the SPE to this. We will look only at the first four FOVs of the healthy sample. We plot both the cell categories of all cells and then the cell types of secretory cells (\(\alpha, \beta\) and \(\delta\) cells) and T-cells (CD8+ and CD4+ T-cells).

df <- data.frame(spatialCoords(spe), colData(spe))

dfSub <- df %>%
    subset(image_name %in% c("E02", "E03", "E04", "E05"))

p <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_category)) +
    geom_point(size= 0.5) +
    facet_wrap(~image_name) +
    theme(legend.title.size = 20, legend.text.size = 20) +
    xlab("x") +
    ylab("y") +
    labs(color = "cell category")+
    coord_equal() +
    theme_light()

dfSub <- dfSub %>%
    subset(cell_type %in% c("alpha", "beta", "delta", "Th", "Tc"))

q <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_type)) +
    geom_point(size= 0.5) +
    facet_wrap(~image_name) +
    theme(legend.title.size = 20, legend.text.size = 20) +
    xlab("x") +
    ylab("y") +
    labs(color = "cell type") +
    coord_equal() +
    theme_light()
wrap_plots(list(p,q), widths = c(1,1), heights = c(1,1), nrow = 2, ncol = 1)

4 Calculating Spatial Statistics Metrics

In a first step, we calculate a spatial statistic curve as implemented by spatstat. One can choose from a range of metrics for discrete marks and calculate these within a mark or between two marks. Common metrics are:

Note that all of these functions have different implementations to correct for inhomogeneity and for comparison between two marks (cross functions) (Baddeley, Rubak, and Turner 2015).

4.1 Correlation

With correlation metrics, we assess the distances of all points to one another while normalising for density effects and of the window size \(|W|\). Furthermore, spatial metrics are corrected for edge effects, due to the fact that points at the border of a FOV do not have a fully-observed neighborhood (Baddeley, Rubak, and Turner 2015, 203 ff.).

A well-known metric is Ripley’s \(K\) function or its variance-stabilised transformation, the \(L\) function. We can calculate a variant of the \(L\) function with the function calcMetricPerFov between e.g \(\alpha\) and cytotoxic T cells. The output is a dataframe with the following most important columns:

  • r: the radius at which the spatial metric is evaluated

  • theo: the theoretical value of a homogeneous (Poisson) realisation of a point process

  • iso: an isotropic edge corrected value of the \(L\) function

metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
                              subsetby = "image_number", fun = "Lcross", 
                              marks = "cell_type",
                              rSeq = seq(0, 50, length.out = 50), 
                              by = c("patient_stage", "patient_id",
                                     "image_number"),
                              ncores = 1)
metricRes %>% head(3)
#>          r     theo border trans iso patient_stage patient_id image_number
#> 1 0.000000 0.000000      0     0   0         Onset       6362            1
#> 2 1.020408 1.020408      0     0   0         Onset       6362            1
#> 3 2.040816 2.040816      0     0   0         Onset       6362            1
#>   npoints centroidx centroidy    fun    selection
#> 1    1000  396.5893  327.1143 Lcross alpha and Tc
#> 2    1000  396.5893  327.1143 Lcross alpha and Tc
#> 3    1000  396.5893  327.1143 Lcross alpha and Tc

We can visualise this metric with plotMetricPerFov function. Here, we need to specify which border correction we want to plot and what the x-axis is. Both can vary from function to function.

# create a unique plotting ID
metricRes$ID <- paste0(
    metricRes$patient_stage, "|", metricRes$patient_id
)
# change levels for plotting
metricRes$ID <- factor(metricRes$ID, levels = c("Non-diabetic|6126", 
                                                "Non-diabetic|6134", 
                                                "Onset|6228","Onset|6362",
                                                "Long-duration|6089", 
                                                "Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "iso", x = "r",
                 imageId = "image_number", ID = "ID", ncol = 2)

By eye, we see no visible difference between the conditions in terms of correlation of \(\alpha\) and cytotoxic T cells.

4.2 Spacing

Another important aspect of spatial analysis is spacing. Here, the shortest distances or empty space to the next neighbor is calculated. This quantifies a different aspect of a point pattern than correlation or intensity of points. Two well-known functions are (Baddeley, Rubak, and Turner 2015, 255–66):

  • nearest-neighbor distance distribution \(G\)

  • empty space function \(F\)

For spacing metrics, we get different border corrections but otherwise the output stays the same:

metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
                              subsetby = "image_number", fun = "Gcross", 
                              marks = "cell_type",
                              rSeq = seq(0, 50, length.out = 50), 
                              by = c("patient_stage", "patient_id",
                                     "image_number"),
                              ncores = 1)
metricRes %>% head(3)
#>          r        theo han rs km hazard patient_stage patient_id image_number
#> 1 0.000000 0.000000000   0  0  0      0         Onset       6362            1
#> 2 1.020408 0.006057893   0  0  0      0         Onset       6362            1
#> 3 2.040816 0.024012273   0  0  0      0         Onset       6362            1
#>   npoints centroidx centroidy    fun    selection
#> 1    1000  396.5893  327.1143 Gcross alpha and Tc
#> 2    1000  396.5893  327.1143 Gcross alpha and Tc
#> 3    1000  396.5893  327.1143 Gcross alpha and Tc
# create a unique plotting ID
metricRes$ID <- paste0(
    metricRes$patient_stage, "|", metricRes$patient_id
)
# change levels for plotting
metricRes$ID <- factor(metricRes$ID, levels = c("Non-diabetic|6126", 
                                                "Non-diabetic|6134", 
                                                "Onset|6228","Onset|6362",
                                                "Long-duration|6089", 
                                                "Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "km", x = "r",
                 imageId = "image_number", ID = "ID", ncol = 2)

In the nearest-neighbor distance function, we see a strong difference between onset T1D, long-duration T1D and non-diabetic controls in terms of spacing of \(\alpha\) and cytotoxic T cells.

5 Functional boxplot

Looking at raw spatial statistics curves can be challenging. In order to summarise this information, we can plot functional boxplots by aggregating the curves into boxplots via a user-defined variable aggregate_by. We use the fbplot function from the fda package (Sun and Genton 2011; Ramsay 2024).

# create a unique ID per row in the dataframe
metricRes$ID <- paste0(
    metricRes$patient_stage, "x", metricRes$patient_id,
    "x", metricRes$image_number
)
#removing field of views that have as a curve only zeros - these are cases where
#there is no cells of one type
metricRes <- metricRes %>% group_by(ID) %>% filter(sum(km) >= 1)

collector <- plotFbPlot(metricRes, "r", "km", "patient_stage")

The functional boxplot shows that onset \(G\)-curves are more variable than the corresponding long-duration and non-diabetic curves. We note as well, that the variability is heteroscedastic along the domain (i.e., variance increases with radius), which is undesirable for our statistical modelling. Therefore, we can e.g. apply a variance stabilising transformation to our data or model this variance in the statistical model.

6 Functional principal component analysis

Another analysis that can be performed is functional principal componentent analysis (fPCA). This is a method to capture the main modes of variation in functional data (Ramsay and Silverman 2005). We use the refund implementation of fPCA.

# filter out all rows that have a constant zero part - all r<10
metricRes <- metricRes %>% filter(r > 10)

# prepare dataframe from calcMetricRes to be in the correct format for pffr
dat <- prepData(metricRes, "r", "km")

# create meta info of the IDs
splitData <- dat$ID %>%
  str_replace("-","_") %>%
  str_split_fixed("x", 3) %>% 
  data.frame(stringsAsFactors = TRUE) %>%
  setNames(c("condition", "patient_id", "imageId")) %>%
  mutate(condition = relevel(condition,"Non_diabetic"))
dat <- cbind(dat, splitData)

# drop rows with NA
dat <- dat |> drop_na()
# calculate the fPCA
pca <- functionalPCA(dat = dat, r = metricRes$r |> unique())
evalues <- pca$evalues
efunctions <- pca$efunctions
# plot the mean curve and the two first eigenfunctions
p_mu <- ggplot(data.frame(r = unique(metricRes$r), mu = pca$mu), 
               aes(x = r, y = mu)) +
    geom_line() +
    theme_light() +
    xlab("r [µm]")

p_efunction1 <- ggplot(data.frame(r = unique(metricRes$r), 
                                  phi1 = pca$efunctions[,1]), 
                       aes(x = r, y = phi1)) +
    geom_line() +
    theme_light() +
    ylim(-0.3,0.3) +
    xlab("r [µm]")

p_efunction2 <- ggplot(data.frame(r = unique(metricRes$r), 
                                  phi2 = pca$efunctions[,2]), 
                       aes(x = r, y = phi2)) +
    geom_line() +
    theme_light() +
    ylim(-0.3,0.3) +
    xlab("r [µm]")

wrap_plots(list(p_mu, p_efunction1, p_efunction2), ncol = 3)

# plot the biplot of the first two PCs
plotFpca(dat = dat, res = pca, colourby = "condition")

In the biplot above we get a very basic differentiation of the \(G\) curves. Onset T1D shows most variability along the first fPC. The second fPC describes less variation.

7 Functional additive mixed models

The \(L\) function above showed no clear difference between the three conditions whereas the \(G\) function showed a strong difference between onset T1D and the two other conditions. In order to test these differences we will use generalised functional additive mixed models. These are generalisations of standard additive mixed models to compare functions over their entire domain. The package that we use is the refund package (Scheipl, Staicu, and Greven 2015; Scheipl, Gertheiss, and Greven 2016).

The model implemented here is of the form:

\[ \mathbb{E}[y_i(r)] = g(\alpha(r) + \beta_{0,g(i)}(r) + \sum_{j=1}^J f_j(X_{ji},r) + \epsilon_i(r)) \]

With the following terms:

For the family we will use a quasi-likelihood distribution where the variance is modeled as quadratic to the mean. We do this to account for the heteroscedastic variance along the domain.

In this context we need to specify a design matrix and contrasts.

library('refund')
# create a design matrix
mm <- model.matrix(~condition, data = dat)
colnames(mm)[1] <- "Intercept"
mm %>% head()
#>   Intercept conditionLong_duration conditionOnset
#> 1         1                      1              0
#> 2         1                      1              0
#> 3         1                      1              0
#> 4         1                      1              0
#> 5         1                      1              0
#> 6         1                      1              0

r <- metricRes$r |> unique()
# fit the model
mdl <- functionalGam(
    data = dat, x = r,
    designmat = mm, weights = dat$npoints,
    formula = formula(Y ~ 1 + conditionLong_duration +
                          conditionOnset + s(patient_id, bs = "re")),
    family = quasi(link = "identity", variance = "mu^2"),
    algorithm = "gam"
)
summary(mdl)
#> 
#> Family: quasi 
#> Link function: identity 
#> 
#> Formula:
#> Y ~ 1 + conditionLong_duration + conditionOnset + s(patient_id, 
#>     bs = "re")
#> 
#> Constant coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 0.207228   0.008404   24.66   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Smooth terms & functional coefficients:
#>                              edf Ref.df      F p-value    
#> Intercept(x)              12.697 19.000 61.891  <2e-16 ***
#> conditionLong_duration(x)  1.304  1.386  2.412  0.0768 .  
#> conditionOnset(x)          1.030  1.036 11.253  0.0006 ***
#> s(patient_id)             15.282 25.000 23.344  <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =  0.349   Deviance explained = 24.7%
#> -REML score = 2445.5  Scale est. = 0.59218   n = 13960 (349 x 40)

plotLs <- lapply(colnames(mm), plotMdl, mdl = mdl,
                 shift = mdl$coefficients[["(Intercept)"]])
#> using seWithMean for  s(x.vec) .
#> using seWithMean for  s(x.vec) .
#> using seWithMean for  s(x.vec) .
wrap_plots(plotLs, nrow = 3, axes = 'collect')

We note that there is a small difference in the \(G\) function between non-diabetic and long-duration T1D samples, but a strong difference between non-diabetic and onset T1D according to the model summary. The point wise confidence bands are a limitation of this method and could be improved with either bootstrapping or continuous confidence bands (Liebl and Reimherr 2023). Thus, we see not only that a spatial difference in co-localisation of \(\alpha\) and cytotoxic T cells is statistically significant but also at which spatial scale this difference occurs.

7.1 Model evaluation

One open problem is the implementation of confidence bands that reflect the non-independently and non-identically distributed residuals. To visualise how much of a problem this is, we can plot the contours of the correlation/covariance and look at some model diagnostics.

resid(mdl) |> cor() |> filled.contour()

resid(mdl) |> cov() |> filled.contour()

try(refund::pffr.check(mdl))

#> 
#> Method: REML   Optimizer: outer newton
#> step failed after 3 iterations.
#> Gradient range [-213.4066,31.77544]
#> (score 2445.476 & scale 0.5921764).
#> Hessian positive definite, eigenvalue range [0.05723549,7192.14].
#> Model rank =  57 / 57 
#> 
#> Basis dimension (k) checking results. Low p-value (k-index<1) may
#> indicate that k is too low, especially if edf is close to k'.
#> 
#>                                    k'   edf k-index p-value  
#> s(x.vec)                        19.00 12.70    0.97   0.025 *
#> s(x.vec):conditionLong_duration  5.00  1.30    0.97   0.025 *
#> s(x.vec):conditionOnset          5.00  1.03    0.97   0.040 *
#> ti(patient_id,x.vec)            27.00 15.28      NA      NA  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In these model diagnostics, we note that there is still some variability in the residuals that is not considered by the model. The QQ plot indicates a good model fit. The residuals show a considerable structure that is in line with the structure in the auto-covariance / correlation plots.

In the functional additive mixed model, we have a specified global intercept varying over the domain \(r\) as well as functional random intercepts varying over the domain \(r\) per grouping variable patient_id. We can plot these smooth estimates of the random intercepts.

# look at the smooth random intercepts per patient
data <- coef(mdl)$smterms$`s(patient_id)`$coef
#> using seWithMean for  s(x.vec) .
data <- data %>% left_join(splitData %>% 
                             select(patient_id, condition) %>% unique)
#> Joining with `by = join_by(patient_id)`

p <- ggplot(data, aes(x.vec, value, colour = patient_id)) +
  geom_point(aes(shape=factor(condition))) +
  theme_light() + 
  geom_smooth(aes(group = 1), col = 'black') +
  xlab("r [µm]")

p
#> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

We note that there is still some structure in these intercepts, which is most likely due to identifiability problems between the global and random intercepts.

sessionInfo()
#> R Under development (unstable) (2025-01-20 r87609)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.21-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] refund_0.1-37               imcdatasets_1.15.0         
#>  [3] cytomapper_1.19.0           EBImage_4.49.0             
#>  [5] SpatialExperiment_1.17.0    SingleCellExperiment_1.29.1
#>  [7] SummarizedExperiment_1.37.0 Biobase_2.67.0             
#>  [9] GenomicRanges_1.59.1        GenomeInfoDb_1.43.4        
#> [11] IRanges_2.41.2              S4Vectors_0.45.2           
#> [13] BiocGenerics_0.53.6         generics_0.1.3             
#> [15] MatrixGenerics_1.19.1       matrixStats_1.5.0          
#> [17] patchwork_1.3.0             stringr_1.5.1              
#> [19] tidyr_1.3.1                 ggplot2_3.5.1              
#> [21] dplyr_1.1.4                 spatialFDA_0.99.6          
#> [23] BiocStyle_2.35.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] gamm4_0.2-6             splines_4.5.0           later_1.4.1            
#>   [4] bitops_1.0-9            filelock_1.0.3          tibble_3.2.1           
#>   [7] svgPanZoom_0.3.4        polyclip_1.10-7         deSolve_1.40           
#>  [10] lifecycle_1.0.4         Rdpack_2.6.2            lattice_0.22-6         
#>  [13] MASS_7.3-64             magrittr_2.0.3          sass_0.4.9             
#>  [16] rmarkdown_2.29          jquerylib_0.1.4         yaml_2.3.10            
#>  [19] httpuv_1.6.15           sp_2.1-4                spatstat.sparse_3.1-0  
#>  [22] minqa_1.2.8             DBI_1.2.3               RColorBrewer_1.1-3     
#>  [25] abind_1.4-8             purrr_1.0.2             RCurl_1.98-1.16        
#>  [28] pracma_2.4.4            rappdirs_0.3.3          GenomeInfoDbData_1.2.13
#>  [31] spatstat.utils_3.1-2    terra_1.8-15            goftest_1.2-3          
#>  [34] spatstat.random_3.3-2   svglite_2.1.3           codetools_0.2-20       
#>  [37] DelayedArray_0.33.4     tidyselect_1.2.1        raster_3.6-31          
#>  [40] UCSC.utils_1.3.1        farver_2.1.2            lme4_1.1-36            
#>  [43] viridis_0.6.5           BiocFileCache_2.15.1    spatstat.explore_3.3-4 
#>  [46] jsonlite_1.8.9          ks_1.14.3               systemfonts_1.2.1      
#>  [49] tools_4.5.0             Rcpp_1.0.14             glue_1.8.0             
#>  [52] gridExtra_2.3           SparseArray_1.7.4       mgcv_1.9-1             
#>  [55] xfun_0.50               HDF5Array_1.35.11       shinydashboard_0.7.2   
#>  [58] withr_3.0.2             BiocManager_1.30.25     fastmap_1.2.0          
#>  [61] boot_1.3-31             rhdf5filters_1.19.0     digest_0.6.37          
#>  [64] R6_2.5.1                mime_0.12               colorspace_2.1-1       
#>  [67] tensor_1.5              jpeg_0.1-10             spatstat.data_3.1-4    
#>  [70] RSQLite_2.3.9           h5mread_0.99.4          pbs_1.1                
#>  [73] httr_1.4.7              htmlwidgets_1.6.4       S4Arrays_1.7.1         
#>  [76] pkgconfig_2.0.3         gtable_0.3.6            blob_1.2.4             
#>  [79] XVector_0.47.2          pcaPP_2.0-5             htmltools_0.5.8.1      
#>  [82] bookdown_0.42           fftwtools_0.9-11        scales_1.3.0           
#>  [85] png_0.1-8               reformulas_0.4.0        spatstat.univar_3.1-1  
#>  [88] knitr_1.49              rjson_0.2.23            magic_1.6-1            
#>  [91] nloptr_2.1.1            nlme_3.1-167            curl_6.2.0             
#>  [94] cachem_1.1.0            rhdf5_2.51.2            RLRsim_3.1-8           
#>  [97] BiocVersion_3.21.1      KernSmooth_2.23-26      parallel_4.5.0         
#> [100] fda_6.2.0               vipor_0.4.7             AnnotationDbi_1.69.0   
#> [103] pillar_1.10.1           grid_4.5.0              vctrs_0.6.5            
#> [106] promises_1.3.2          dbplyr_2.5.0            xtable_1.8-4           
#> [109] cluster_2.1.8           beeswarm_0.4.0          evaluate_1.0.3         
#> [112] tinytex_0.54            magick_2.8.5            mvtnorm_1.3-3          
#> [115] cli_3.6.3               locfit_1.5-9.10         compiler_4.5.0         
#> [118] rlang_1.1.5             crayon_1.5.3            grpreg_3.5.0           
#> [121] labeling_0.4.3          mclust_6.1.1            fds_1.8                
#> [124] ggbeeswarm_0.7.2        stringi_1.8.4           rainbow_3.8            
#> [127] viridisLite_0.4.2       deldir_2.0-4            BiocParallel_1.41.0    
#> [130] nnls_1.6                hdrcde_3.4              munsell_0.5.1          
#> [133] Biostrings_2.75.3       tiff_0.1-12             spatstat.geom_3.3-5    
#> [136] Matrix_1.7-2            ExperimentHub_2.15.0    bit64_4.6.0-1          
#> [139] Rhdf5lib_1.29.0         KEGGREST_1.47.0         shiny_1.10.0           
#> [142] AnnotationHub_3.15.0    rbibutils_2.3           memoise_2.0.1          
#> [145] bslib_0.9.0             bit_4.5.0.1

Ali, Mayar, Merel Kuijs, Soroor Hediyeh-Zadeh, Tim Treis, Karin Hrovatin, Giovanni Palla, Anna C Schaar, and Fabian J Theis. 2024. “GraphCompass: Spatial Metrics for Differential Analyses of Cell Organization Across Conditions.” Bioinformatics.

Baddeley, Adrian, Ege Rubak, and Rolf Turner. 2015. Spatial Point Patterns: Methodology and Applications with R. CRC press.

Baddeley, Adrian, and Rolf Turner. 2005. “spatstat: An R Package for Analyzing Spatial Point Patterns.” Journal of Statistical Software 12 (6): 1–42.

Canete, Nicolas P, Sourish S Iyengar, John T Ormerod, Heeva Baharlou, Andrew N Harman, and Ellis Patrick. 2022. “SpicyR: Spatial Analysis of in Situ Cytometry Data in R.” Bioinformatics.

Damond, Nicolas, Stefanie Engler, Vito R. T. Zanotelli, Denis Schapiro, Clive H. Wasserfall, Irina Kusmartseva, Harry S. Nick, et al. 2019. “A Map of Human Type 1 Diabetes Progression by Imaging Mass Cytometry.” Cell Metabolism 29 (3): 755–768.e5.

Goldsmith, Jeff, Fabian Scheipl, Lei Huang, Julia Wrobel, Chongzhi Di, Jonathan Gellar, Jaroslaw Harezlak, et al. 2024. Refund: Regression with Functional Data.

Liebl, Dominik, and Matthew Reimherr. 2023. “Fast and Fair Simultaneous Confidence Bands for Functional Parameters.” Journal of the Royal Statistical Society Series B: Statistical Methodology 85 (3): 842–68.

Ramsay, James. 2024. Fda: Functional Data Analysis. https://CRAN.R-project.org/package=fda.

Ramsay, JO, and BW Silverman. 2005. “Principal Components Analysis for Functional Data.” Functional Data Analysis, 147–72.

Scheipl, Fabian, Jan Gertheiss, and Sonja Greven. 2016. “Generalized Functional Additive Mixed Models.” Electronic Journal of Statistics 10 (1): 1455–92.

Scheipl, Fabian, Ana-Maria Staicu, and Sonja Greven. 2015. “Functional Additive Mixed Models.” Journal of Computational and Graphical Statistics 24 (2): 477–501.

Sun, Ying, and Marc G Genton. 2011. “Functional Boxplots.” Journal of Computational and Graphical Statistics.

Wrobel, Julia, Alex C Soupir, Mitchell T Hayes, Lauren C Peres, Thao Vu, Andrew Leroux, and Brooke L Fridley. 2024. “Mxfda: A Comprehensive Toolkit for Functional Data Analysis of Single-Cell Spatial Data.” Bioinformatics Advances 4 (1): vbae155.