The weighted pair correlation function (wPCF)

The weighted pair correlation function (wPCF) is a spatial statistic that extends the standard pair correlation function to identify spatial relationships between two point populations, one of which is labelled with a continuous mark. For full details, see https://doi.org/10.1371/journal.pcbi.1010994 or https://doi.org/10.1017/S2633903X24000011.

In this tutorial, we briefly reproduce the example from Figure 3 of Quantification of spatial and phenotypic heterogeneity in an agent-based model of tumour-macrophage interactions (PLOS Computational Biology 19(3) https://doi.org/10.1371/journal.pcbi.1010994). We start by placing a line of points from population \(A\) across the line \(y = 0\), and a second population \(B\) placed randomly across the domain. Population \(B\) points are assigned a mark \(m_i = y_i^2\)

import numpy as np
import matplotlib.pyplot as plt
import muspan as ms

# Set the figure resolution
plt.rcParams['figure.dpi'] = 270

# Define the number of points in populations A and B
N_A = 100
N_B = 1000

# Create population A - N_A points across the domain center
points_A = np.asarray([(v / N_A, 0) for v in np.linspace(1, N_A, N_A)])

# Create population B - N_B points distributed via complete spatial randomness
points_B = np.random.rand(N_B, 2) * [1, 2] - [0, 1]

# Assign labels to population B points based on their y-coordinate squared
labels_B = points_B[:, 1] ** 2

# Create a domain for the wPCF example
domain = ms.domain('wPCF example')

# Add points to the domain
domain.add_points(points_B, 'B')
domain.add_points(points_A, 'A')

# Add labels to population B points
domain.add_labels('Mark', labels_B, add_labels_to='B')

# Visualize the domain with points colored by their mark
ms.visualise.visualise(domain, color_by='Mark', figure_kwargs={'figsize': (8, 6)})

(<Figure size 2160x1620 with 2 Axes>, <Axes: >)

# Query population A from the domain
popA = ms.query.query(domain, ('collection',), 'is', 'A')

# Query population B from the domain
popB = ms.query.query(domain, ('collection',), 'is', 'B')

# Calculate the weighted pair correlation function (wPCF)
# Parameters:
# - domain: the domain containing the points
# - popA: the first population of points (A)
# - popB: the second population of points (B)
# - mark_pop_B: the label/mark associated with population B
# - max_R: the maximum radius to consider for the wPCF
# - annulus_step: the step size for the annulus radii
# - annulus_width: the width of each annulus
radii, wPCF = ms.spatial_statistics.weighted_pair_correlation_function(
    domain, popA, popB, mark_pop_B='Mark', max_R=1, annulus_step=0.1, annulus_width=0.15
)

# Visualize the weighted pair correlation function (wPCF)
ms.visualise.visualise_wpcf(radii, wPCF, figure_kwargs={'figsize': (8, 6)})

# Plot the expected value of the target mark M (M = r^2) as a function of r
plt.gca().plot(np.arange(0, 1, 0.01), np.arange(0, 1, 0.01) ** 2, color=[0, 0, 0], linestyle=':', lw=5)

[<matplotlib.lines.Line2D at 0x320ab9a50>]

The black dashed line shows the expected value of the target mark \(M\) for which the wPCF is maximised (\(M = r^2\)), as a function of \(r\).