---
title: "Functional Data Analysis of Spatial Metrics"
author: 
  - name: "Martin Emons"
    affiliation:
      - &DMLS Department of Molecular Life Sciences, University of Zurich, Switzerland
      - &SIB SIB Swiss Institute of Bioinformatics, Lausanne, Switzerland
    email: "martin.emons@uzh.ch"
  - name: Fabian Scheipl
    affiliation:
      - &LMU Department of Statistics, Ludwig-Maximilians-Universität München, Germany
  - name: Mark D. Robinson
    affiliation:
      - *DMLS
      - *SIB
package: "`r BiocStyle::Biocpkg('spatialFDA')`"
output:
  BiocStyle::html_document
abstract: >
  This vignette shows how to perform a custom spatial analysis with the package `spatialFDA`. 
vignette: >
  %\VignetteIndexEntry{Functional Data Analysis of Spatial Metrics}
  %\VignetteEncoding{UTF-8}  
  %\VignetteEngine{knitr::rmarkdown}
bibliography: spatialFDA.bib
editor_options: 
  chunk_output_type: console
---

```{r v1, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    eval = TRUE
)
```

# Introduction

This vignette demonstrates how to use `r BiocStyle::Biocpkg('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 cell types given several samples/conditions. It does so by calculating spatial statistics metrics via the `r BiocStyle::CRANpkg('spatstat')` package and comparing differences using functional additive mixed models as implemented in the `r BiocStyle::CRANpkg('refund')` package [@spatstat2005; @refund2024].

The use case is a dataset from @damondMapHumanType2019 which contains images from 12 human donors. The raw data is published under a `CC-BY-4.0` License on [Mendeley](https://data.mendeley.com/datasets/cydmwsfztj/2). 

This package is similar to other packages in `python` and `R`. The following table shows the main differences in terms of functionality [@ali2024graphcompass; @canete2022spicyr; @wrobel2024mxfda]. 

|Package name                          | Foundation | Testing procedure   |
| ------------------------------------ | --------------- | ------------------- |
| `r BiocStyle::Biocpkg('spicyR')`     | $L$-function    | Scalar comparison   |
| `GraphCompass`                       | Graph-based     | Graph and  scalar comparison | 
| `r BiocStyle::CRANpkg('mxfda')`      | $K$-, $G$- and $L$- function | Function as input for survival modelling |
| `r BiocStyle::Biocpkg('spatialFDA')` | most `r BiocStyle::CRANpkg('spatstat')` functions   | Functional comparison over domain |

# Installation

`r BiocStyle::Biocpkg('spatialFDA')` can be installed and loaded from Bioconductor as follows

```{r installation, include = TRUE, eval = FALSE}
if (!requireNamespace("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install("spatialFDA")
```

```{r setup, warning = FALSE, message = FALSE}
library("spatialFDA")
library("dplyr")
library("ggplot2")
library("tidyr")
library("stringr")
library("dplyr")
library("patchwork")
library("SpatialExperiment")

set.seed(1234)
```

# Getting started

In this vignette we will analyse a diabetes dataset acquired by imaging mass cytometry (IMC) as acquired by @damondMapHumanType2019. 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 [@damondMapHumanType2019].

## Loading the data

The @damondMapHumanType2019 dataset is easily loaded from `ExperimentHub` via a small reader function `.loadExample()`. The entire dataset can be loaded by setting `full = TRUE`. For computational reasons, one can reduce to three patients as well by setting this flag to `FALSE`. We will subset the entire dataset to two samples per condition in order to have a multi-condition/multi-sample setting. The package offers multiple datatypes, we will use the `r BiocStyle::Biocpkg('SpatialExperiment')` (SPE) object [@righelli2022spatialexperiment].

```{r loading, warning = FALSE, message = FALSE}
# retrieve example data from Damond et al. (2019)
spe <- .loadExample(full = TRUE)

spe <- subset(spe, ,patient_id %in% c(6089,6180,6126,6134,6228,6414))
# set cell types as factors
colData(spe)$cell_type <- as.factor(colData(spe)$cell_type) 
```

## 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).

```{r plotting fovs, warning = FALSE, fig.width=8, fig.height=15}
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)
```

# Calculating Spatial Statistics Metrics

In a first step, we calculate a spatial statistic curve as implemented by `r BiocStyle::CRANpkg('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:

- Ripley's $K$ function and its variance stabilised form, Besag's $L$

- Pair correlation function $g$

- Nearest-neighbour function $G$

- Empty space function $F$

Note that all of these functions have different implementations to correct for inhomogeneity and for comparison between two marks (cross functions) [@baddeleySpatialPointPatterns].

All spatial statistics metrics are a function of radius $r$. As such, we have to determine the maximum radius sequence to consider. [@baddeleySpatialPointPatterns] recommend to take a maximum radius $r$ that is not larger than $1/4$ of the window length. As the window can be non-convex, we calculate the half of the minimum bounding radius of the observation window as heuristic.

```{r}
p <- rMaxHeuristic(spe,
subsetby = "image_number", marks = "cell_type"
)
p
```

We will consider a rMax of $50$ for analysing shorter-range interactions.

## 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 [@baddeleySpatialPointPatterns, pp. 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


```{r Lfunction, warning = FALSE, message = FALSE}
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)
```

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.

```{r plotLfunction, warning = FALSE, fig.width=8, fig.height=8}
# 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|6414",
                                                "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.

## 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 [@baddeleySpatialPointPatterns, pp. 255-266]:

- nearest-neighbor distance distribution $G$

- empty space function $F$

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

```{r Gfunction, warning = FALSE, message = FALSE}
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 plotGfunction, warning = FALSE, fig.width=8, fig.height=8}
# 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|6414",
                                                "Long-duration|6089",
                                                "Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "rs", 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. 

# 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 `r BiocStyle::CRANpkg('fda')` package [@sun2011functional; @ramsay2024fda].

```{r, funcBoxPlot, warning = FALSE, results='hide'}
# 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 %>% dplyr::group_by(ID) %>% dplyr::filter(sum(rs) >= 1)

collector <- plotFbPlot(metricRes, "r", "rs", "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.

```{r, variancetransform, warning = FALSE}
# can determine with a boxcox transformation what is the ideal parameter
# for the transformation use Fisher's variance-stabilising transformation
metricRes$rs <- asin(sqrt(metricRes$rs))

collector <- plotFbPlot(metricRes, 'r', 'rs', 'patient_stage')
```

# 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 [@ramsayPrincipalComponentsAnalysis2005]. We use the `r BiocStyle::CRANpkg('refund')` implementation of fPCA. 

```{r fPCA, warning = FALSE}
# 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", "rs", sample_id = "patient_id",
 image_id = "image_number", condition = "patient_stage")

dat$condition <- gsub("-","_", dat$patient_stage) %>% 
  as.factor() %>%
  relevel("Non_diabetic")

# drop rows with NA
dat <- dat |> drop_na()
# calculate the fPCA
pca <- functionalPCA(dat = dat, r = metricRes$r |> unique(), pve = 0.995)
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")
# print the eigenvalues
evalues
```

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. 

# 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 `r BiocStyle::CRANpkg('refund')` package [@scheiplFunctionalAdditiveMixed2015; @scheiplGeneralizedFunctionalAdditive2016; @refund2024].

The model implemented here is of the form:

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

With the following terms:

- $y_i(r)$: functional response, here the `r BiocStyle::CRANpkg('spatstat')` curves

- $g$: optional link function 

- $\alpha(r)$: global functional intercept varying over the domain $r$

- $\beta_{0,g(i)}(r)$: random functional intercept varying over the domain $r$ per grouping variable $g(i)$.

- $f_j(X_{ji},r)$: additive predictors

For the family we will use a scaled $t$-distribution with a logarithmic link function. This distribution has a positive support, thereby modelling the strictly positive response of the spatial statistics functions. 

In this context we need to specify a design matrix and contrasts. For the functional random intercepts we define a smooth function for each patient ID as implemented in `refund` [@refund2024].

```{r funcGamG, fig.height=10, warning = FALSE}
library('refund')
# create a design matrix
mm <- model.matrix(~condition, data = dat)
mm %>% head()

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 = mgcv::scat(link = "log"),
    algorithm = "bam"
)
summary(mdl)

plotLs <- lapply(colnames(mm), plotMdl, mdl = mdl,
                 shift = mdl$coefficients[["(Intercept)"]])
wrap_plots(plotLs, nrow = 3, axes = 'collect')
```


We note that there is a small non-significant 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 [@liebl2023fast]. 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.

## 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 model diagnostics.

```{r contour, warning = FALSE}
resid(mdl) |> cor() |> filled.contour(levels = seq(-1, 1, l = 40))
resid(mdl) |> cov() |> filled.contour()

qqnorm(resid(mdl), pch = 16)
qqline(resid(mdl))
```

In these model diagnostics, we note that there is still some variability in the residuals that is not considered by the model. The Q-Q plot indicates a good but not perfect 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 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.

```{r intercept, warning = FALSE, eval = TRUE}
# look at the smooth random intercepts per patient
data <- coef(mdl)$smterms$`character(0)`$coef

data <- data %>% left_join(dat %>% 
                             select(patient_id, condition) %>% unique)

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

p
```

We note that these random errors are not constrained to sum-to-zero over the domain $r$. This can lead to problems of identifiability between the global intercept and the functional random intercepts. 

```{r sessionInfo}
sessionInfo()
```
