Contents

1 Introduction

multipointR is a package to fit point process models from spatstat.model on SpatialExperiment and SpatialFeatureExperiment objects. Cells are approximated as points and their distribution in space is modelled as the log-linear combination of spatial (or non-spatial) covariates. multipointR has a mode to do inferences on the level of a single image with fitModel and across multiple images with fitModelAcrossImages.

2 Installation

multipointR can be installed and loaded from Bioconductor as follows

if (!requireNamespace("BiocManager")) {
    install.packages("BiocManager")
}
BiocManager::install("multipointR")

3 Setup

library("multipointR")
library("SpatialExperiment")
library("SpatialFeatureExperiment")
library("dplyr")
library("ggplot2")
library("spatstat.model")
library("patchwork")

4 Dataset

For this document, we will use the well established MIBI TOF dataset from Keren et al. (2018).

spe <- SpatialDatasets::spe_Keren_2018()
spe
#> class: SpatialExperiment 
#> dim: 48 197678 
#> metadata(0):
#> assays(1): intensities
#> rownames(48): Na Si ... Ta Au
#> rowData names(0):
#> colnames(197678): 1 2 ... 197677 197678
#> colData names(40): CellID imageID ... Censored sample_id
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x y
#> imgData names(0):

The preprocessing and visualisation of this dataset is taken from the Statial vignette.

# Code source: https://www.bioconductor.org/packages/release/bioc/vignettes/
# Statial/inst/doc/Statial.html
# #kontextual-identifying-discrete-changes-in-cell-state

# Examine all cell types in image
unique(spe$cellType)
#>  [1] "Keratin_Tumour" "dn_T_CD3"       "B_cell"         "CD4_T_cell"    
#>  [5] "DC_or_Mono"     "Unidentified"   "Macrophages"    "CD8_T_cell"    
#>  [9] "Other_Immune"   "Endothelial"    "Mono_or_Neu"    "Mesenchymal"   
#> [13] "Neutrophils"    "NK"             "Tumour"         "DC"            
#> [17] "Tregs"

# Set up cell populations
tumour <- c("Keratin_Tumour", "Tumour")

bcells <- c("B_cell")
tcells <- c("dn_T_CD3", "CD4_T_cell", "CD8_T_cell", "Tregs")
myeloid <- c("DC_or_Mono", "DC", "Mono_or_Neu", "Macrophages", "Neutrophils")

endothelial <- c("Endothelial")
mesenchymal <- c("Mesenchymal")

tissue <- c(endothelial, mesenchymal)
immune <- c(bcells, tcells, myeloid, "NK", "other immune")

all <- c(tumour, tissue, immune, "Unidentified")

# Lets define a new cell type vector
spe$cellTypeNew <- spe$cellType

# Select for all cells that express higher than baseline level of p53
p53Pos <- assay(spe)["p53", ] > -0.300460

# Find p53+ tumour cells
spe$cellTypeNew[spe$cellType %in% tumour] <- "Tumour"
spe$cellTypeNew[p53Pos & spe$cellType %in% tumour] <- "p53_Tumour"

# Group all immune cells under the name "Immune"

spe$cellTypeNew[spe$cellType %in% immune] <- "Immune"

spe$cellTypeNew <- as.factor(spe$cellTypeNew)
speSub <- subset(spe, , cellTypeNew %in% c("Immune", "Tumour", "p53_Tumour"))
speSub$cellTypeNew <- factor(speSub$cellTypeNew,
    levels = c("Immune", "Tumour", "p53_Tumour")
)

5 Inference with parametric point process models

5.1 The inhomogeneous Poisson point process

Parametric point process models (PPMs), as implemented in the spatstat package, model the distribution of points in space as an (in)homogeneous Poisson distribution (Baddeley, Turner, and others 2014). Thus, the expected number of points falling in a region \(B\) follows an inhomogeneous Poisson process with rate \(\lambda(u)\), the local intensity (Baddeley et al. 2016).

\[ \mathbb{E}[n(\mathbf{X} \cap B)] = \int_B \lambda(u)du \]

The PPM then models this inhomogeneous intensity \(\lambda(u)\) as a function of spatial and non-spatial covariates

\[ \lambda_{\Theta}(u) = \exp(\beta_0(u) + \Theta^TZ(u)) \]

where \(\beta_0(u)\) is the baseline intensity, \(Z(u)\) are spatially varying covariate functions and \(\Theta\) is the parameter vector to be estimated. In particular, parameters within \(\Theta\) would be relevant for making inferences in order to make claims about spatial associations.

5.2 The Gibbs point process

In the Poisson point process, points can be arbitrarily close together, which in a tissue of cells is not feasible. However, a constraint on cell co-localization can be imposed by a Gibbs point process. This class of models allows inhibitory or repulsive interactions between points to be modeled. In our case, we model the interaction between cells as a so-called “hard core” process, whereby the resulting conditional intensity at a location \(u\) is given by

\[ \lambda(u \mid \mathbf{x}) = \cases{ \phi(u) & \text{if $u$ is permissible}\\ 0 & \text{if $u$ is not permissible}\\ } \]

where \(\phi(u)\) is the intensity function of an inhomogeneous Poisson process (note the change in notation to before). The hardcore permission is given if the distance of \(u\) to any other point \(x_i\) is greater than the hard core diameter \(r\). This diameter \(r\) can be either provided by the user (e.g. prior knowledge on the average cell size) or estimated from the data (Baddeley et al. 2016).

The estimation from a dataset then becomes the following:

\[ \lambda_{\Theta}(u \mid \mathbf{x}) = \exp(\beta_0(u) + \Theta^TZ(u)) \cdot h(u,r,\mathbf{x}) \]

where \(\Theta\) are the first order terms (similar to the Poisson point process), whereas \(h(u,r,\mathbf{x})\), the second order terms, define the hard core interaction with interaction radius \(r\), subject to the constraint:

\[ h(u,r,\mathbf{x}) = \cases{ 0 & if $\lVert u-v \rVert \leq r, v \in \mathbf{x}\setminus\{u\}$\\ 1 & if $\lVert u-v \rVert > r, v \in \mathbf{x}\setminus\{u\}$\\ } \]

(Baddeley et al. 2016; Baddeley and Turner 2000).

5.3 Single image inference

For the single image case, we will analyse the relationship between the distribution of p53+ tumour cells and immune cells in image \(6\).

# Code source: https://www.bioconductor.org/packages/release/bioc/vignettes/
# Statial/inst/doc/Statial.html
# #kontextual-identifying-discrete-changes-in-cell-state

# Plot image 6

df <- spe |>
    colData() |>
    cbind(spatialCoords(spe)) |>
    as.data.frame() |>
    dplyr::filter(imageID == "6") |>
    dplyr::filter(cellTypeNew %in%
        c("Immune", "Tumour", "p53_Tumour"))

df$cellTypeNew <- factor(df$cellTypeNew,
    levels = c("Immune", "Tumour", "p53_Tumour")
)
p1 <- df |>
    arrange(cellTypeNew) |>
    ggplot(aes(x = x, y = y, color = cellTypeNew)) +
    geom_point(size = 1) +
    scale_colour_manual(
        values = c("#505050", "#D6D6D6", "#64BC46"),
        labels = c("Immune", "Tumour", "p53+ Tumour")
    ) +
    guides(colour = guide_legend(
        title = "Cell types",
        override.aes = list(size = 5)
    )) +
    coord_equal() +
    theme_light()

p1

We notice qualitatively a clearly separated tumour with p53+ cells and immune cells at the tumour border.

5.3.1 Model of homogeneous intensity

First, we will estimate the distribution of p53+ tumour cells as a homogeneous intensity in space, i.e., we fit a model of the form

\[ \lambda(u) = \exp(\beta_0) \cdot h(u,r,\mathbf{x}) \]

with a constant intercept \(\beta_0\)

speSub <- subset(spe, , imageID == "6")

m0 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    interaction = "Hardcore",
    formula = as.formula("p53_Tumour ~ 1")
)

The model is very basic and has not many parameters, let’s look at the spatial trend:

plot(m0)

The resulting spatial trend is a flat surface.

This is not a very sensible model since we see clear a inhomogeneous distribution of points in space.

Looking at the conditional intensity function we obtain an estimate of the intensity of the process given the point pattern.

plot(m0, type = "cif")

We note that in the vicinity of cells no other point can be placed, a result of the hardcore process.

This leads to a conceptual question on how to parametrise the Gibbs model best for biological tissue. An alternative to the hard core interaction model is the Strauss interaction model that parametrises not a hard cut-off to zero but rather an interaction probability \(\gamma\). A middle ground is the Fiksel process or a Strauss-Hardcore process that combines both a hard core effect for effects \(r_h\) where no cells are found and a Strauss/Fiksel process where cells are less likely to be found. This can either be a fixed probability \(\gamma\) (Strauss-Hard) or follow a double exponential decay (Fiksel). We will use a Fiksel process for this tutorial as it is more flexible in modelling interactions (Takacs and Fiksel 1986).

\[ h(u,r_h,r_f,\mathbf{x}) = \cases{ 0 & if $\lVert u-v \rVert \leq r_h, v \in \mathbf{x}\setminus\{u\}$\\ \exp(a\exp(-\kappa d)) & if $r_h < \lVert u-v \rVert \leq r_f, v \in \mathbf{x}\setminus\{u\}$\\ 1 & if $\lVert u-v \rVert > r_f, v \in \mathbf{x}\setminus\{u\}$\\ } \]

where \(r_h\) is the hard core radius (estimated from the data) and \(r_f\) is the Fiksel interaction radius (user provided) (Baddeley et al. 2016).

5.3.2 Model of inhomogeneous intensity

Next, we will formulate a null model of the distribution of p53 positive cells. For this model, we will specify an inhomogeneous intensity varying with a bivariate spline model of \(x\) and \(y\) (\(s(x,y)\)). The model we fit is of the form

\[ \lambda(u) = \exp(\beta_0(u)) \cdot Z_0(u) \cdot h(u,r_h,r_f,\mathbf{x}) \]

with a spatially-varying intercept \(\beta_0(u)\), which we specify as an power image of the intensity \(Z_0(u)\) and an inhomogeneous intercept parameterised as a tensor product spline.

m1 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ log(lambda) + 
        splines::ns(x, df = 3)*splines::ns(y, df = 3)"),
    interaction = "Hardcore",
    improve.type = "enet",
    improve.args = list(alpha = 1),
    relaxed = TRUE
)
#> Warning in .resolve_control(control = control, nvars = nvars, deprecated =
#> list(thresh = if (!missing(thresh)) thresh, : Passing 'thresh' to glmnet() is
#> deprecated. Use control = list(thresh = ...) instead.
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit

We can look at the influence of the parameters with a likelihood ratio test. This is a bit less verbose than the individual \(z\)-tests for each basis function that summary(m1) would provide.

anova(m1, test = "LRT")
#> Warning: anova.ppm now computes the *adjusted* deviances when the models are
#> not Poisson processes.
#> Analysis of Deviance Table
#> 
#> Terms added sequentially (first to last)
#> 
#> Model 1: ~splines::ns(x, df = 3) + splines::ns(y, df = 3)     Hardcore
#> Model 2: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3)   Hardcore
#> Model 3: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3) + splines::ns(x, df = 3):splines::ns(y, df = 3)   Hardcore
#>   Npar Df AdjDeviance  Pr(>Chi)    
#> 1    7                             
#> 2    8  1     1706.18 < 2.2e-16 ***
#> 3   17  9       37.84  1.86e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Given the degree to which the log-likelihood improves by adding the spline intensity function (i.e., the LRT), we can conclude that model of inhomogeneity is more suitable than the homogeneous model for this dataset.

plot(m1)

This is also clear from the spatial trend, highlighting that most regions have a very low probability of p53+ cells.

5.3.3 Model of distance to immune cells

From the literature, we know that p53+ tumour cells have immune-modulatory function. Our hypothesis is therefore that p53+ tumour cells will interact spatially with immune cells. In order to test this, we formulate a Gibbs model of distance to immune cells while accounting for an underlying inhomogeneous distribution of p53+ tumour cells.

The model is then:

\[ \lambda(u) = \exp(\beta_0 + \beta_{\text{dist}}Z_{\text{dist}}(u)) \cdot Z_0(u) \cdot h(u,r_h,r_f,\mathbf{x}) \]

m2 <- fitModel(
    spe = speSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ log(lambda) + 
        splines::ns(x, df = 3)*splines::ns(y, df = 3) + 
        spatstat.geom::distfun(Immune)"),
    interaction = "Hardcore",
    improve.type = "enet",
    improve.args = list(alpha = 1),
    relaxed = TRUE
)
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit

We can again look at the contribution of the model parameters with a LRT

anova(m1, m2, test = "LRT")
#> Analysis of Deviance Table
#> 
#> Model 1: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3) + splines::ns(x, df = 3):splines::ns(y, df = 3)   Hardcore
#> Model 2: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3) + spatstat.geom.distfun.Immune. + splines::ns(x, df = 3):splines::ns(y, df = 3)   Hardcore
#>   Npar Df AdjDeviance Pr(>Chi)
#> 1   17                        
#> 2   18  1     0.74223   0.3889

We see that the model including the distance to immune cells does not lead to an improvement in model fit for this specific image.

plot(m2)

Also the spatial trend looks pretty similar as in the model above.

An important step in fitting complex spatial models is to check goodness-of-fit

diagnose.ppm(m2)

#> Model diagnostics (raw residuals)
#> Diagnostics available:
#>  four-panel plot
#>  mark plot 
#>  smoothed residual field
#>  x cumulative residuals
#>  y cumulative residuals
#>  sum of all residuals
#> sum of raw residuals in clipped window = -8.225e-10
#> area of clipped window = 3864000
#> quadrature area = 3931000
#> range of smoothed field =  [-2.374e-05, 1.555e-05]

We note that there is still some unexplained variance of the residuals along both the \(x\) and \(y\) coordinates but overlaying both in a 2D density the deviations are only minor and centered around zero.

We can also look at the residuals of the model which are in this case not completely normally distributed, indicating some residual unexplained variance. For computational reasons we will not run this

p <- qqplot.ppm(m2, nsim = 50)

Looking at a simulated Q-Q plot of the residuals, we see slightly heavier tails than expected.

Given the complexity of the dataset, the model fit is acceptable.

5.3.4 Model of segmented objects

Another option that users have is to define a segmented polygon as spatial covariate in their ppm model.

To do this, we first perform a segmentation on our image. For simplicity, we will use the Bioconductor package sosta. Of course, one can use any other segmentation method, the only requirement being that the polygon is stored as an sf object.

segmentedTumour <- sosta::reconstructShapeDensityImage(
    speSub,
    marks = "cellTypeNew",
    markSelect = c("Tumour", "p53_Tumour"),
    thres = 5e-4
)

p1 + geom_sf(
    data = segmentedTumour, inherit.aes = FALSE,
    fill = NA, color = "red", linewidth = 1
)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.

In order to save this polygon we will convert our SpatialExperiment object into a SpatialFeatureExperiment object and store it in the annotGeometries.

sfeSub <- toSpatialFeatureExperiment(speSub)
annotGeometry(sfeSub, "tumour_mask") <- segmentedTumour

After this, we can specify our ppm with the segmented tumour as spatial covariate. The result will be a logical covariate, indicating the intensity of the response within the tumour (tumour_mask == TRUE) as compared to outside.

m3 <- fitModel(
    spe = sfeSub,
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ log(lambda) + 
        splines::ns(x, df = 3)*splines::ns(y, df = 3) + 
        tumour_mask + 
        spatstat.geom::distfun(Immune)"),
    interaction = "Hardcore",
    improve.type = "enet",
    improve.args = list(alpha = 1),
    relaxed = TRUE
)
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit

We will check again for the importance of this covariate with a LRT

anova(m2, m3, test = "LRT")
#> Analysis of Deviance Table
#> 
#> Model 1: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3) + spatstat.geom.distfun.Immune. + splines::ns(x, df = 3):splines::ns(y, df = 3)   Hardcore
#> Model 2: ~log(lambda) + splines::ns(x, df = 3) + splines::ns(y, df = 3) + tumour_mask + spatstat.geom.distfun.Immune. + splines::ns(x, df = 3):splines::ns(y, df = 3)     Hardcore
#>   Npar Df AdjDeviance  Pr(>Chi)    
#> 1   18                             
#> 2   19  1      25.378 4.712e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We can interpret from this test that the covariate tumour_mask is important for describing the distribution of p53+ cells in space. This makes sense, since p53+ tumour cells are a subset of tumour cells and can only be found within the tumour.

plot(m3) + geom_sf(
    data = segmentedTumour, inherit.aes = FALSE,
    fill = NA, color = "darkred", linewidth = 0.6
)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.

The model looks very similar to the spatial trend above in the model with just the distance to immune cells defined. We can again look at the improvement of the model diagnostics

diagnose.ppm(m3)

#> Model diagnostics (raw residuals)
#> Diagnostics available:
#>  four-panel plot
#>  mark plot 
#>  smoothed residual field
#>  x cumulative residuals
#>  y cumulative residuals
#>  sum of all residuals
#> sum of raw residuals in clipped window = -9.887e-10
#> area of clipped window = 3864000
#> quadrature area = 3931000
#> range of smoothed field =  [-2.176e-05, 1.401e-05]

The model diagnostics look very similar to the model above. For computational reasons, we did not include the Q-Q plot here, but the result is also similar to the Q-Q plot above.

5.4 Multi-image inference

5.4.1 Fitting a model per image

As mentioned above, this dataset contains many images. We may want to fit separate models like those above to each image and look at the resulting distribution of the coefficients relating the distance to immune cells \(\beta_{\text{dist}}\)

mdlLs <- fitModelAcrossImages(
    spe = spe,
    imageId = "imageID",
    marks = "cellTypeNew",
    formula = as.formula("p53_Tumour ~ log(lambda) + 
        splines::ns(x, df = 3)*splines::ns(y, df = 3) + 
        spatstat.geom::distfun(Immune)"),
    interaction = "Hardcore",
    improve.type = "enet",
    improve.args = list(alpha = 1),
    relaxed = TRUE,
    threshold = 10
)
#> Fitting ppm to image 1
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 2
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 3
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 4
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 5
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 6
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 7
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 8
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 9
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 10
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 11
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 12
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 13
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 14
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 15
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 16
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 17
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 18
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 19
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 20
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 21
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 22
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 23
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 24
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 25
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 26
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 27
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: algorithm did not converge
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 28
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 29
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 31
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 32
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 33
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 34
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 35
#> There were less than 10 points to compute an intensity on
#> Fitting ppm to image 36
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 37
#> Warning: glm.fit: algorithm did not converge
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 38
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 39
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: algorithm did not converge
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 40
#> Model was fit with relaxed enet. The number of coefficients can
#>         therefore be different than an unregularised fit
#> Fitting ppm to image 41
#> There were less than 10 points to compute an intensity on

The dataset describes three tumour subtypes, cold, compartmentalised and mixed. We want to investigate the differences in the distance to immune cells coefficient \(\beta_{\text{dist}}\) across these subtypes. First, we need to extract the coefficients from the model and convert this to a data.frame.

mdlDf <- mdlToDf(
    mdlLs = mdlLs,
    imageCovariates = c(
        "imageID",
        "tumour_type"
    )
)
mdlDfSub <- mdlDf %>% filter(covariate %in% c("spatstat.geom.distfun.Immune."))

ggplot(mdlDfSub, aes(x = covariate, y = Estimate, label = imageID)) +
    geom_boxplot(outlier.shape = NA, alpha = 0.3) +
    geom_jitter(aes(color = log10(S.E.)),
        position = position_jitter(seed = 123)
    ) +
    geom_text(aes(color = log10(S.E.)),
        hjust = 0,
        vjust = 0, position = position_jitter(seed = 123)
    ) +
    theme_light() +
    facet_wrap(~tumour_type)

We see a difference between the tumour types: cold tumours show no effect and a generally positive effect of immune cell distance on p53+ cells in compartmentalised and mixed tumours despite not being significant.

There are two of these models that are clear outliers, let’s look at them in more detail

# model 12
plot(mdlLs[[12]])

diagnose.ppm(mdlLs[[12]])

#> Model diagnostics (raw residuals)
#> Diagnostics available:
#>  four-panel plot
#>  mark plot 
#>  smoothed residual field
#>  x cumulative residuals
#>  y cumulative residuals
#>  sum of all residuals
#> sum of raw residuals in clipped window = -1.783e-06
#> area of clipped window = 3882000
#> quadrature area = 3934000
#> range of smoothed field =  [-4.974e-05, 6.233e-05]
# model 39
plot(mdlLs[[39]])

diagnose.ppm(mdlLs[[39]])

#> Model diagnostics (raw residuals)
#> Diagnostics available:
#>  four-panel plot
#>  mark plot 
#>  smoothed residual field
#>  x cumulative residuals
#>  y cumulative residuals
#>  sum of all residuals
#> sum of raw residuals in clipped window = 3.564e-11
#> area of clipped window = 3845000
#> quadrature area = 3935000
#> range of smoothed field =  [-6.873e-06, 5.516e-06]

The first model seems reasonable in terms of diagnostics, however it contains many cells and thus shows a stronger effect than the other images. The second image shows larger variation in the residuals at both higher \(x\) and \(y\) values, suggesting that the model is misspecified for this image.

5.4.2 Second-stage model of coefficients

We can also test the effect of tumour type (cold, compartmentalised, mixed) on the distance to immune cells coefficient \(\beta_{\text{dist}}\) with a linear model. We weight each observation by the inverse of the standard error, as done by Gerber et al. (2026).

The model is then parametrised as follows:

\[ {\hat{\beta}}_{\text{dist}} = \gamma_0 + \gamma Z_{\text{stage}} + \epsilon \quad \epsilon \sim \mathcal{N}\!\left(0,\, \sigma^2 \text{SE}_{\text{dist}}\right) \]

We can fit this model

mdl <- lm(Estimate ~ tumour_type,
    data = mdlDf,
    weights = 1 / (mdlDf$S.E.),
    subset = covariate == "spatstat.geom.distfun.Immune."
)

print(summary(mdl))
#> 
#> Call:
#> lm(formula = Estimate ~ tumour_type, data = mdlDf, subset = covariate == 
#>     "spatstat.geom.distfun.Immune.", weights = 1/(mdlDf$S.E.))
#> 
#> Weighted Residuals:
#>      Min       1Q   Median       3Q      Max 
#> -0.18076 -0.06589  0.00171  0.04652  0.56867 
#> 
#> Coefficients:
#>                                Estimate Std. Error t value Pr(>|t|)  
#> (Intercept)                  -0.0002697  0.0014541  -0.186   0.8544  
#> tumour_typecompartmentalised  0.0005215  0.0027152   0.192   0.8493  
#> tumour_typemixed              0.0045231  0.0024613   1.838   0.0785 .
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.1541 on 24 degrees of freedom
#>   (2 observations deleted due to missingness)
#> Multiple R-squared:  0.1287, Adjusted R-squared:  0.05611 
#> F-statistic: 1.773 on 2 and 24 DF,  p-value: 0.1914

In agreement with the findings above, there are positive but non-significant effects for both compartmentalised and mixed tumours in comparison to cold tumours.

5.4.3 Fitting a shared model across images

We can also fit a shared model across all images. In that case, we parametrise an interaction effect between the tumour type and the distance to immune cells. Furthermore, we specify a random intercept per sample ID.

In comparison to the ppm models from spatstat.model mppm offers less flexibility in terms of model parametrisation. If you need more flexibility, we recommend you to consider the previous section.

For instance, gam fits are not available in mppm, therefore we account for inhomogeneity with a kernel density estimation of the point distribution. As offsets are also not supported by mppm we can add the KDE as a parameter to our estimation via a power law.

mdl <- fitModelAcrossImages(
    spe = spe,
    imageId = "imageID",
    imageLs = seq(15, 18),
    marks = "cellTypeNew",
    interaction = "Fiksel",
    formula = as.formula("p53_Tumour ~  log(lambda) +
        splines::bs(x, df = 3):splines::bs(y, df = 3) +
        tumour_type:spatstat.geom::distfun(Immune) + 
        (1|DONOR_NO)"),
    improve.type = "enet",
    improve.args = list(alpha = 1),
    sharedModel = TRUE,
    threshold = 10
)
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: algorithm did not converge
#> Warning: glm.fit: algorithm did not converge
#> Warning in formula.character(object, env = baseenv()): Using formula(x) is deprecated when x is a character vector of length > 1.
#>   Consider formula(paste(x, collapse = " ")) instead.
#> (Fixing subset index..)
#> iteration 1 
#> iteration 2 
#> iteration 3
mdl |> summary()
#> Point process model fitted to 4 point patterns
#> Call: spatstat.model::mppm(formula = feFormula, random = reFormula,
#> Log trend formula: ~log(lambda) + splines::bs(x, df = 3):splines::bs(y, df = 3) 
#> + tumour_type:spatstat.geom.distfun.Immune.
#> 
#> Fixed effects:
#>                                                (Intercept) 
#>                                                8.373577457 
#>                                                log(lambda) 
#>                                                1.950965820 
#>              tumour_typecold:spatstat.geom.distfun.Immune. 
#>                                               -0.000577975 
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. 
#>                                               -0.003209734 
#>             tumour_typemixed:spatstat.geom.distfun.Immune. 
#>                                                0.005086109 
#> Random effects:
#>       (Intercept)
#> 30734  -0.2455181
#> 30740   1.8399851
#> 30742  -0.1420693
#> 30744  -1.4523976
#> 
#> Fitted coefficients of interpoint interaction:
#>                                               x 
#>                                    8.946268e-05 
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1 
#>                                   -9.591993e-01 
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1 
#>                                    1.220953e-01 
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1 
#>                                   -1.181115e+00 
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2 
#>                                    1.146289e+00 
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2 
#>                                    5.049595e-01 
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2 
#>                                   -1.385607e-01 
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3 
#>                                   -7.066477e-01 
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3 
#>                                    1.074888e-01 
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3 
#>                                    7.432865e-01 
#> All fitted coefficients:
#>                                                (Intercept) 
#>                                               8.373577e+00 
#>                                                log(lambda) 
#>                                               1.950966e+00 
#>                                                          x 
#>                                               8.946268e-05 
#>            splines::bs(x, df = 3)1:splines::bs(y, df = 3)1 
#>                                              -9.591993e-01 
#>            splines::bs(x, df = 3)2:splines::bs(y, df = 3)1 
#>                                               1.220953e-01 
#>            splines::bs(x, df = 3)3:splines::bs(y, df = 3)1 
#>                                              -1.181115e+00 
#>            splines::bs(x, df = 3)1:splines::bs(y, df = 3)2 
#>                                               1.146289e+00 
#>            splines::bs(x, df = 3)2:splines::bs(y, df = 3)2 
#>                                               5.049595e-01 
#>            splines::bs(x, df = 3)3:splines::bs(y, df = 3)2 
#>                                              -1.385607e-01 
#>            splines::bs(x, df = 3)1:splines::bs(y, df = 3)3 
#>                                              -7.066477e-01 
#>            splines::bs(x, df = 3)2:splines::bs(y, df = 3)3 
#>                                               1.074888e-01 
#>            splines::bs(x, df = 3)3:splines::bs(y, df = 3)3 
#>                                               7.432865e-01 
#>              tumour_typecold:spatstat.geom.distfun.Immune. 
#>                                              -5.779750e-04 
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. 
#>                                              -3.209734e-03 
#>             tumour_typemixed:spatstat.geom.distfun.Immune. 
#>                                               5.086109e-03 
#> 
#> Random effects summary:
#> Random effects:
#>  Formula: ~1 | DONOR_NO
#>         (Intercept) Residual
#> StdDev:    1.381674        1
#> 
#> Variance function:
#>  Structure: fixed weights
#>  Formula: ~invwt 
#> 
#> Interaction for each pattern:    Fiksel process
#> 
#> Interaction on row 1:
#> Fiksel process
#> Interaction distance:    24.55344
#> Hard core distance:  21.95344
#> Rate parameter:  0.5
#> Fitted interaction strength a:   NaN
#> 
#> Interaction on row 2:
#> Fiksel process
#> Interaction distance:    10.02976
#> Hard core distance:  7.929762
#> Rate parameter:  0.5
#> Fitted interaction strength a:   NaN
#> 
#> Interaction on row 3:
#> Fiksel process
#> Interaction distance:    14.60179
#> Hard core distance:  10.00179
#> Rate parameter:  0.5
#> Fitted interaction strength a:   NaN
#> 
#> Interaction on row 4:
#> Fiksel process
#> Interaction distance:    17.57489
#> Hard core distance:  12.97489
#> Rate parameter:  0.5
#> Fitted interaction strength a:   NaN
#> 
#> --- Gory details: ---
#> Combined data frame has 12132 rows
#> Linear mixed-effects model fit by maximum likelihood
#>   Data: moadf 
#>   Subset: glmmsubset 
#>   AIC BIC    logLik
#>    NA  NA -52195.91
#> 
#> Random effects:
#>  Formula: ~1 | DONOR_NO
#>         (Intercept)  Residual
#> StdDev:    1.182787 0.8560531
#> 
#> Variance function:
#>  Structure: fixed weights
#>  Formula: ~invwt 
#> Fixed effects:  .mpl.Y ~ log(lambda) + splines::bs(x, df = 3):splines::bs(y,      df = 3) + tumour_type:spatstat.geom.distfun.Immune. + x 
#>                                                                Value Std.Error
#> (Intercept)                                                 8.373577 0.7079399
#> log(lambda)                                                 1.950966 0.0433769
#> x                                                           0.000089 0.0000789
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1            -0.959199 0.6991599
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1             0.122095 0.8842056
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1            -1.181115 0.4997039
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2             1.146289 0.9633558
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2             0.504959 1.2042996
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            -0.138561 0.6774567
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            -0.706648 0.5703537
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.107489 0.8162226
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.743287 0.5665547
#> tumour_typecold:spatstat.geom.distfun.Immune.              -0.000578 0.0020791
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.003210 0.0006516
#> tumour_typemixed:spatstat.geom.distfun.Immune.              0.005086 0.0014542
#>                                                               DF  t-value
#> (Intercept)                                                11570 11.82809
#> log(lambda)                                                11570 44.97712
#> x                                                          11570  1.13326
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1            11570 -1.37193
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1            11570  0.13808
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1            11570 -2.36363
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2            11570  1.18989
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2            11570  0.41930
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            11570 -0.20453
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            11570 -1.23896
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3            11570  0.13169
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3            11570  1.31194
#> tumour_typecold:spatstat.geom.distfun.Immune.              11570 -0.27799
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. 11570 -4.92570
#> tumour_typemixed:spatstat.geom.distfun.Immune.             11570  3.49752
#>                                                            p-value
#> (Intercept)                                                 0.0000
#> log(lambda)                                                 0.0000
#> x                                                           0.2571
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1             0.1701
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1             0.8902
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1             0.0181
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2             0.2341
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2             0.6750
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2             0.8379
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3             0.2154
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.8952
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.1896
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.7810
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.0000
#> tumour_typemixed:spatstat.geom.distfun.Immune.              0.0005
#>  Correlation: 
#>                                                            (Intr) lg(lm) x     
#> log(lambda)                                                 0.530              
#> x                                                          -0.036  0.084       
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1            -0.070 -0.070  0.216
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1             0.035  0.030 -0.345
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1            -0.049 -0.129 -0.533
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2             0.004  0.020  0.010
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2             0.017  0.032  0.057
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            -0.015 -0.017  0.087
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            -0.030 -0.015  0.137
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.005 -0.019 -0.175
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.037  0.056 -0.266
#> tumour_typecold:spatstat.geom.distfun.Immune.              -0.061  0.021 -0.035
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.219 -0.389 -0.041
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.027  0.044  0.127
#>                                                            s::(,d=3)1:::(,d=3)1
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1            -0.741              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1             0.313              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2            -0.719              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2             0.608              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            -0.328              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3             0.541              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3            -0.407              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.162              
#> tumour_typecold:spatstat.geom.distfun.Immune.              -0.014              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.008              
#> tumour_typemixed:spatstat.geom.distfun.Immune.              0.025              
#>                                                            s::(,d=3)2:::(,d=3)1
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1            -0.301              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2             0.587              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2            -0.771              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2             0.465              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            -0.373              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.435              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3            -0.100              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.041              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.062              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.093              
#>                                                            s::(,d=3)3:::(,d=3)1
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2            -0.302              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2             0.441              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            -0.692              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3             0.191              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3            -0.183              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.540              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.015              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.109              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.096              
#>                                                            s::(,d=3)1:::(,d=3)2
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2            -0.812              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2             0.468              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            -0.688              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.593              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3            -0.284              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.001              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.096              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.021              
#>                                                            s::(,d=3)2:::(,d=3)2
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2            -0.706              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3             0.605              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3            -0.670              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.387              
#> tumour_typecold:spatstat.geom.distfun.Immune.              -0.023              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.024              
#> tumour_typemixed:spatstat.geom.distfun.Immune.              0.060              
#>                                                            s::(,d=3)3:::(,d=3)2
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3            -0.369              
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3             0.485              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3            -0.695              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.004              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.009              
#> tumour_typemixed:spatstat.geom.distfun.Immune.              0.009              
#>                                                            s::(,d=3)1:::(,d=3)3
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3            -0.886              
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3             0.462              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.001              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.082              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.015              
#>                                                            s::(,d=3)2:::(,d=3)3
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3            -0.592              
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.003              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.018              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.006              
#>                                                            s::(,d=3)3:::(,d=3)3
#> log(lambda)                                                                    
#> x                                                                              
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                                
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3                                
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3                                
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3                                
#> tumour_typecold:spatstat.geom.distfun.Immune.               0.022              
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.  0.029              
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.060              
#>                                                            tmr_typcl:...I.
#> log(lambda)                                                               
#> x                                                                         
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3                           
#> tumour_typecold:spatstat.geom.distfun.Immune.                             
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune. -0.010         
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.004         
#>                                                            tmr_typcm:...I.
#> log(lambda)                                                               
#> x                                                                         
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)1                           
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)2                           
#> splines::bs(x, df = 3)1:splines::bs(y, df = 3)3                           
#> splines::bs(x, df = 3)2:splines::bs(y, df = 3)3                           
#> splines::bs(x, df = 3)3:splines::bs(y, df = 3)3                           
#> tumour_typecold:spatstat.geom.distfun.Immune.                             
#> tumour_typecompartmentalised:spatstat.geom.distfun.Immune.                
#> tumour_typemixed:spatstat.geom.distfun.Immune.             -0.014         
#> 
#> Standardized Within-Group Residuals:
#>           Min            Q1           Med            Q3           Max 
#> -1.883041e+00 -4.146539e-01 -5.915583e-02 -7.511065e-06  1.733827e+01 
#> 
#> Number of Observations: 11588
#> Number of Groups: 4

We notice again a significant effect of the interaction of compartmentalised tumour with distance to immune cells as well as the mixed tumour types. The effect sizes are in both cases very subtle.

In order to assess the model fit we can calculate the sum of the residuals per image and plot them as a function of the condition. In expectation these should be zero.

res <- residuals(mdl, type = "raw")
resInt <- sapply(res, integral.msr)
print(resInt)
#>           1           2           3           4 
#>   -1.363447 -176.501711   -4.309934  -11.789741

These residuals are not perfectly centered around zero for some of the models, implying that the model fit from mppm is not ideal in this scenario.

6 Session Information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.24-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] SpatialDatasets_1.11.0          ExperimentHub_3.3.2            
#>  [3] AnnotationHub_4.3.2             BiocFileCache_3.3.0            
#>  [5] dbplyr_2.6.0                    patchwork_1.3.2                
#>  [7] spatstat.model_3.7-2            rpart_4.1.27                   
#>  [9] spatstat.explore_3.8-2          nlme_3.1-170                   
#> [11] spatstat.random_3.5-1           spatstat.geom_3.8-2            
#> [13] spatstat.univar_3.2-0           spatstat.data_3.1-9            
#> [15] ggplot2_4.0.3                   dplyr_1.2.1                    
#> [17] SpatialFeatureExperiment_1.15.0 SpatialExperiment_1.23.0       
#> [19] SingleCellExperiment_1.35.2     SummarizedExperiment_1.43.0    
#> [21] Biobase_2.73.2                  GenomicRanges_1.65.1           
#> [23] Seqinfo_1.3.1                   IRanges_2.47.3                 
#> [25] S4Vectors_0.51.9                BiocGenerics_0.59.12           
#> [27] generics_0.1.4                  MatrixGenerics_1.25.0          
#> [29] matrixStats_1.5.0               multipointR_0.99.6             
#> [31] BiocStyle_2.41.0               
#> 
#> loaded via a namespace (and not attached):
#>   [1] splines_4.6.1             bitops_1.1-0             
#>   [3] filelock_1.0.3            tibble_3.3.1             
#>   [5] R.oo_1.27.1               polyclip_1.10-7          
#>   [7] lifecycle_1.0.5           httr2_1.3.0              
#>   [9] Rdpack_2.6.6              formula.tools_1.7.1      
#>  [11] sf_1.1-2                  edgeR_4.11.6             
#>  [13] lattice_0.23-1            MASS_7.3-66              
#>  [15] backports_1.5.1           magrittr_2.0.5           
#>  [17] limma_3.69.4              sass_0.4.10              
#>  [19] rmarkdown_2.31            jquerylib_0.1.4          
#>  [21] yaml_2.3.12               otel_0.2.0               
#>  [23] sp_2.2-3                  spatstat.sparse_3.2-0    
#>  [25] DBI_1.3.0                 RColorBrewer_1.1-3       
#>  [27] multcomp_1.4-32           abind_1.4-8              
#>  [29] spatialreg_1.4-3          purrr_1.2.2              
#>  [31] R.utils_2.13.0            RCurl_1.98-1.20          
#>  [33] TH.data_1.1-5             rappdirs_0.3.4           
#>  [35] sandwich_3.1-3            spatstat.utils_3.2-4     
#>  [37] terra_1.9-46              units_1.0-1              
#>  [39] goftest_1.2-3             dqrng_0.4.1              
#>  [41] DelayedMatrixStats_1.35.0 codetools_0.2-20         
#>  [43] DropletUtils_1.33.0       DelayedArray_0.39.6      
#>  [45] scuttle_1.23.1            shape_1.4.6.1            
#>  [47] tidyselect_1.2.1          farver_2.1.2             
#>  [49] jsonlite_2.0.0            BiocNeighbors_2.7.3      
#>  [51] e1071_1.7-17              iterators_1.0.14         
#>  [53] survival_3.8-11           sosta_1.5.1              
#>  [55] smoothr_1.3.0             foreach_1.5.2            
#>  [57] tools_4.6.1               Rcpp_1.1.2               
#>  [59] glue_1.8.1                BiocBaseUtils_1.15.1     
#>  [61] SparseArray_1.13.2        xfun_0.60                
#>  [63] mgcv_1.9-4                EBImage_4.55.2           
#>  [65] HDF5Array_1.41.3          withr_3.0.3              
#>  [67] BiocManager_1.30.27       fastmap_1.2.0            
#>  [69] boot_1.3-32               rhdf5filters_1.25.4      
#>  [71] spData_2.3.5              digest_0.6.39            
#>  [73] R6_2.6.1                  wk_0.9.5                 
#>  [75] LearnBayes_2.15.2         tensor_1.5.1             
#>  [77] jpeg_0.1-11               dichromat_2.0-1          
#>  [79] RSQLite_3.53.3            R.methodsS3_1.8.2        
#>  [81] h5mread_1.5.2             data.table_1.18.6.1      
#>  [83] class_7.3-24              httr_1.4.8               
#>  [85] htmlwidgets_1.6.4         S4Arrays_1.13.0          
#>  [87] spdep_1.4-2               pkgconfig_2.0.3          
#>  [89] gtable_0.3.6              blob_1.3.0               
#>  [91] S7_0.2.2                  XVector_0.53.0           
#>  [93] htmltools_0.5.9           bookdown_0.47            
#>  [95] fftwtools_0.9-11          scales_1.4.0             
#>  [97] png_0.1-9                 reformulas_0.4.4         
#>  [99] knitr_1.51                rjson_0.2.23             
#> [101] coda_0.19-4.1             curl_8.0.0               
#> [103] proxy_0.4-29              cachem_1.1.0             
#> [105] zoo_1.9-0                 rhdf5_2.57.12            
#> [107] operator.tools_1.6.3.1    BiocVersion_3.24.0       
#> [109] KernSmooth_2.23-27        parallel_4.6.1           
#> [111] AnnotationDbi_1.75.2      s2_1.1.11                
#> [113] pillar_1.11.1             grid_4.6.1               
#> [115] vctrs_0.7.3               beachmat_2.29.1          
#> [117] sfheaders_0.4.5           evaluate_1.0.5           
#> [119] tinytex_0.60              zeallot_0.2.0            
#> [121] magick_2.9.1              mvtnorm_1.4-2            
#> [123] cli_3.6.6                 locfit_1.5-9.12          
#> [125] compiler_4.6.1            crayon_1.5.3             
#> [127] rlang_1.3.0               labeling_0.4.3           
#> [129] classInt_0.4-11           viridisLite_0.4.3        
#> [131] deldir_2.0-4              BiocParallel_1.47.0      
#> [133] Biostrings_2.81.6         tiff_0.1-12              
#> [135] marginaleffects_0.32.0    glmnet_5.0               
#> [137] Matrix_1.7-6              sparseMatrixStats_1.25.0 
#> [139] bit64_4.8.4               Rhdf5lib_2.1.0           
#> [141] KEGGREST_1.53.6           statmod_1.5.2            
#> [143] rbibutils_2.4.1           memoise_2.0.1            
#> [145] bslib_0.12.0              bit_4.6.0

References

Baddeley, Adrian, Ege Rubak, Rolf Turner, and others. 2016. Spatial Point Patterns: Methodology and Applications with R. Vol. 1. CRC press Boca Raton.

Baddeley, Adrian, and Rolf Turner. 2000. “Practical Maximum Pseudolikelihood for Spatial Point Patterns: (With Discussion).” Australian & New Zealand Journal of Statistics 42 (3): 283–322.

Baddeley, Adrian, Rolf Turner, and others. 2014. “Package ‘Spatstat’.” The Comprehensive R Archive Network () 146.

Gerber, Reto, Jake Griner, Silvia Guglietta, Carsten Krieg, and Mark D Robinson. 2026. “MIMIC: A Flexible Pipeline to Register and Summarize Imc-Msi Experiments.” Communications Biology.

Keren, Leeat, Marc Bosse, Diana Marquez, Roshan Angoshtari, Samir Jain, Sushama Varma, Soo-Ryum Yang, et al. 2018. “A Structured Tumor-Immune Microenvironment in Triple Negative Breast Cancer Revealed by Multiplexed Ion Beam Imaging.” Cell 174 (6): 1373–87.

Takacs, Roland, and Th Fiksel. 1986. “Interaction Pair-Potentials for a System of Ant’s Nests.” Biometrical Journal 28 (8): 1007–13.