# ANOVA for Data in Metric Spaces,
with Applications to Spatial Point Patterns

Raoul Müller111Work supported by DFG RTG 2088. 222raoul.mueller@uni-
goettingen.de 333Institute for Mathematical Stochastics, University of
Göttingen, 37077 Göttingen, Germany. Dominic Schuhmacher33footnotemark: 3
Jorge Mateu444Work partially funded by Ministry of Science and Innovation with
grant PID2019-107392RB-I00, and by grant UJI-B2021-37 from University Jaume I.
555Department of Mathematics, University Jaume I, 12071 Castellón, Spain.

###### Abstract

We give a review of recent ANOVA-like procedures for testing group differences
based on data in a metric space and present a new such procedure. Our
statistic is based on the classic Levene’s test for detecting differences in
dispersion. It uses only pairwise distances of data points and and can be
computed quickly and precisely in situations where the computation of
barycenters (“generalized means”) in the data space is slow, only by
approximation or even infeasible. We show the asymptotic normality of our test
statistic and present simulation studies for spatial point pattern data, in
which we compare the various procedures in a 1-way ANOVA setting. As an
application, we perform a 2-way ANOVA on a data set of bubbles in a mineral
flotation process.

##  1 Introduction

Real-world statistical data is often not Euclidean, involving components that
are most suitably analyzed in a more complicated space. Examples include
spaces of point patterns and more general subsets, trees and more general
graphs, functions and images.

In recent years a number of methods have been proposed for analyzing group
differences of such data by generalizing classical analysis of variance
(ANOVA) ideas to more complex data spaces. Examples include Cuevas et al.,
(2004) for functional data, Huckemann et al., (2009) for data on Riemannian
manifolds and Ramón et al., (2016) for point pattern data. A common feature of
the underlying spaces is that there is typically a more or less natural
concept of distance between data points available. In addition to the more
obvious choices of distances on function spaces and Riemannian manifolds,
suitable metrics for tree spaces, graph spaces and point pattern spaces can be
found in Billera et al., (2001), Ginestet et al., (2017) and Müller et al.,
(2020), respectively.

In the present paper we focus on generalized ANOVA-procedures for metric
spaces without using any more special structure of the space. There is a
number of preceding articles that work in similar generality.

Anderson, (2001) proposes to perform ANOVA based on pairwise dissimilarities
of observations rather than Euclidean distances between observations and their
group means, and introduces the name PERMANOVA for this procedure. While not
directly referring to any more abstract spaces than $\mathbb{R}^{d}$, that
article clearly discusses the abstract template of doing non-Euclidean ANOVA
without using a centroid object. We discuss this further in Subsection 3.1.
Anderson, (2006) proposes multidimensional scaling followed by a Levene’s test
(using the centroid object in the principal coordinate space) for detecting
differences of within-group dispersions (scatter, variability); this is
referred to as PERMDISP, see Anderson, (2017). Anderson et al., (2017) and
Hamidi et al., (2019) correct the PERMANOVA statistic for heteroscedasticity
in the unbalanced setting based on the variants of classical ANOVA by
Brown–Forsythe and Welch, respectively.

In an independent line of research, Dubey and Müller, (2019) design an ANOVA
procedure on metric spaces using Fréchet means as centroid objects. They
propose to use as statistic the sum of an ANOVA-term and a Levene-term. We
discuss this further in Subsection 3.2.

In the present paper we formulate Anderson’s PERMANOVA on general metric
spaces. We simply refer to the resulting method as Anderson ANOVA, because the
use of M (due to the use of $\mathbb{R}^{d}$ in Anderson’s work) seems
inappropriate in our context and the use of PER (referring to the fact that a
permutation test is performed) does not distinguish it from the other methods
used. Rather than pursuing the PERMDISP method mentioned above, we introduce a
new test for detecting differences of within-group dispersion based on
Levene’s procedure and refer to it as $L$-test. Our test statistic works
directly with the pairwise distances between observations without using any
kind of group centroid, neither in the original metric space nor in any
principal coordinate space. We show that it has an asymptotic
$\chi_{1}^{2}$-distribution, but we recommend using it with a permutation test
just as the other statistics.

We also study the two summands used by Dubey and Müller, (2019) as separate
test statistics for detecting differences in location and dispersion,
respectively. We refer to Table 1 for an overview of the methods discussed.

| location | dispersion
---|---|---
pairwise distances | Anderson, Subsec. 3.1 | New $L$-test, Section 4
Fréchet means | Dubey–Müller, Subsec. 3.2 | Dubey–Müller, Subsec. 3.2
Table 1:  Overview of the non-Euclidean ANOVA methods studied in this paper.
Procedures targeting _location_ are derived from the classic ANOVA statistic,
whereas those targeting _dispersion_ are derived from the classic Levene’s
statistic (ANOVA statistic for “deviations”). The rows distinguish whether
computationally a procedure is based on (simple arithmetics of) _pairwise
distances_ or on a centroid object (here a _Fréchet mean_) in the metric
space.

Although the methods described are applicable in general metric spaces, our
central goal in undertaking this research was to be able to perform ANOVA for
point pattern data, see also the discussion section of Müller et al., (2020).
Among all the metric spaces mentioned above, we therefore focus in the later
part of the present paper on the space of finite point patterns equipped with
the TT-metric from Müller et al., (2020). As in many other spaces, exact
Fréchet means can be computed within reasonable time only for (very) small
data sets and one typically has to resort to a heuristic algorithm that finds
only local minima of the Fréchet functional. We present simulation studies to
compare the powers of the four tests across various situations and to
understand the quality of approximation by the limiting
$\chi_{1}^{2}$-distribution from a practical point of view. We also present an
application of a 2-way ANOVA on a data set of bubbles in a mineral flotation
process.

The plan of the paper is as follows: Section 2 contains a brief reminder of
central aspects of classical ANOVA including Levene’s test. In Section 3 we
give a rather detailed presentation of Anderson ANOVA in metric spaces and the
two summands proposed by Dubey and Müller. In Section 4 we introduce our new
L-statistic, discuss its relation to the other methods and the original
Levene’s test, and show its asymptotic distribution. Section 5 is a short
overview of the metric space of point patterns. In Sections 6 and 7 we present
the simulation studies and the real-world data example. The paper ends with
some further conclusions in Section 8.

##  2 Classic ANOVA

For self-containedness and easy reference we briefly remind the reader of some
facts and formulae in the context of the classical ANOVA going back to Fisher,
(1925). Details can be found in Scheffé, (1967).

#### One-Way ANOVA

Given independent observations $x_{ij}\in\mathbb{R}$, $1\leq j\leq n_{i}$,
$1\leq i\leq k$, from $k$ potentially different distributions
$P_{1},\ldots,P_{k}$, we do the following sum-of-squares decomposition

$\mathrm{TSS}=\mathrm{MSS}+\mathrm{RSS},$

where

$\displaystyle\mathrm{TSS}$
$\displaystyle=\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}(x_{ij}-\bar{x}_{\cdot\cdot})^{2}$
(total sum of squares) $\displaystyle\mathrm{MSS}$
$\displaystyle=\sum_{i=1}^{k}n_{i}(\bar{x}_{i\cdot}-\bar{x}_{\cdot\cdot})^{2}$
(model sum of squares) $\displaystyle\mathrm{RSS}$
$\displaystyle=\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}(x_{ij}-\bar{x}_{i\cdot})^{2}$
(residual sum of squares).

Here $\bar{x}_{i\cdot}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}x_{ij}$ denotes the
$i$-th group mean and
$\bar{x}_{\cdot\cdot}=\frac{1}{n}\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}x_{ij}$
denotes the overall mean. We write $n=\sum_{i=1}^{k}n_{i}$ for the total
number of observations.

Assume for now that the group distributions $P_{i}$ are Gaussian with the same
variance. Under the null hypothesis that also the means are the same (hence
all data comes from the same normal distribution), it is well-known that the
ANOVA statistics

$F=\frac{n-k}{k-1}\frac{\mathrm{MSS}}{\mathrm{RSS}},$ (1)

describing the ratio between the variability explained by the model and the
total variability in the data, is $F$-distributed with $k-1$ and $n-k$ degrees
of freedom. Since $T\sim F(d_{1},d_{2})$ implies
$d_{1}T\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\chi^{2}_{d_{1}}$
as $d_{2}\to\infty$, we obtain

$(k-1)F\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\chi^{2}_{k-1}\quad\text{as
$n\to\infty$.}$ (2)

The _asymptotic_ result remains true even if $P_{1}=P_{2}=\ldots=P_{k}$ is
non-Gaussian, but has second moments and there are
$\lambda_{1},\ldots,\lambda_{k}>0$ such that the ratios of group sizes satisfy
$\frac{n_{i}}{n}\to\lambda_{i}$, see e.g. Wooldridge, (2010), Section 3.6.2.

######  Remark 1.

Strictly speaking ANOVA techniques are designed for inference within a linear
model of different group means plus errors. Based on an error distribution $P$
with mean zero, one considers the model equations

$x_{ij}=\mu_{i}+\varepsilon_{ij},\quad 1\leq j\leq n_{i},\,1\leq i\leq k,$

where $\mu_{i}\in\mathbb{R}$ are the different group means and
$\varepsilon_{ij}$ are i.i.d. $P$-distributed error terms. In terms of the
group distributions above this means that $P_{i}=P\ast\delta_{\mu_{i}}$, i.e.
$P_{i}$ is obtained by shifting $P$ by $\mu_{i}$. Note that the asymptotic
$\chi^{2}_{k-1}$-test does not need this assumption since in any case the null
hypothesis just correspond to having $k$ times the same distribution. At the
same time we cannot expect this test to achieve high power against all
alternatives that have substantially different group distributions (see also
the paragraph on Levene’s test below). We will take up this point when
discussing ANOVA on metric spaces, where typically “shifting the distribution”
is meaningless (but may have an intuitive counterpart).

#### Two-Way ANOVA

As soon as more than one grouping factor is involved, important design
decisions come into play, such as if factors are (partially) nested or if we
allow for interaction terms between several factors on the same level. ANOVA
has a long-standing history with many different designs. As an example which
is pursued further in later sections we remind the reader of the balanced two-
way ANOVA (two main factors, with interaction terms, same number $\tilde{n}$
of observations for each factor combination).

Given independent observations $x_{i_{1}i_{2}j}\in\mathbb{R}$, $1\leq
j\leq\tilde{n}$, $1\leq i_{1}\leq k_{1}$, $1\leq i_{2}\leq k_{2}$ from groups
obtained by crossing a Factor $a$ with $k_{1}$ levels and a Factor $b$ with
$k_{2}$ levels (with $n_{i_{1}i_{2}}:=\tilde{n}$ observations for each
combination), we can perform a finer sum-of-squares decomposition

$\mathrm{TSS}=\mathrm{SSa}+\mathrm{SSb}+\mathrm{SSi}+\mathrm{RSS},$

splitting up the model sum of squares into sums of squares for the individual
factors and an interaction sum of squares. In formulae:

$\displaystyle\mathrm{TSS}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\sum_{j=1}^{\tilde{n}}(x_{i_{1}i_{2}j}-\bar{x}_{\cdot\cdot\cdot})^{2}$
$\displaystyle\mathrm{RSS}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\sum_{j=1}^{\tilde{n}}(x_{i_{1}i_{2}j}-\bar{x}_{i_{1}i_{2}\cdot})^{2}$
$\displaystyle\mathrm{SSa}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}k_{2}\tilde{n}(\bar{x}_{i_{1}\cdot\cdot}-\bar{x}_{\cdot\cdot\cdot})^{2}$
$\displaystyle\mathrm{SSb}$
$\displaystyle=\sum_{i_{2}=1}^{k_{2}}k_{1}\tilde{n}(\bar{x}_{\cdot
i_{2}\cdot}-\bar{x}_{\cdot\cdot\cdot})^{2}$ $\displaystyle\mathrm{SSi}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\tilde{n}(\bar{x}_{i_{1}i_{2}\cdot}-\bar{x}_{i_{1}\cdot\cdot}-\bar{x}_{\cdot
i_{2}\cdot}+\bar{x}_{\cdot\cdot\cdot})^{2},$

where the various means are taken over the dot components while keeping the
given indices fixed. Set
$n=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}n_{i_{1}i_{2}}=k_{1}k_{2}\tilde{n}$.

In addition to performing an omnibus test for group differences as for one-way
ANOVA, we may then test for effects of Factor a and b separately, as well as
for an interaction effect. The corresponding statistics are

$\displaystyle
Fa=\frac{n-k_{1}k_{2}}{k_{1}-1}\frac{\mathrm{SSa}}{\mathrm{RSS}},\quad
Fb=\frac{n-k_{1}k_{2}}{k_{2}-1}\frac{\mathrm{SSb}}{\mathrm{RSS}},\quad
Fi=\frac{n-k_{1}k_{2}}{(k_{1}-1)(k_{2}-1)}\frac{\mathrm{SSi}}{\mathrm{RSS}}.$

If the observations come from Gaussian distributions with equal variances,
each of the three statistics is $F$-distributed again under the corresponding
null hypothesis that different levels of the factor or interaction to be
tested do not lead to different shifts in mean. The degrees of freedom can be
read from the denominator and the numerator, respectively, of the first ratio
in each statistic.

#### Levene’s Test

The test first proposed in Levene, (1960) was originally developed as a
preliminary test to check for equal variances _before_ applying the basic
ANOVA $F$-test in the Gaussian setting. This was important, as it was well-
known at the time that for the goal of inference about differences in the
means of the various groups (see Remark 1), the size of the $F$-test can
depart substantially from its nominal size if group variances are not equal.

Levene, (1960) proposed to use as test statistic the usual ANOVA statistic,
but to replace the observations $x_{ij}$ by the absolute differences from
their group means $z_{ij}=\lvert x_{ij}-\bar{x}_{i\cdot}\rvert$, i.e.

$\displaystyle\widetilde{F}=\frac{n-k}{k-1}\cdot\frac{\sum_{i=1}^{k}n_{i}(\bar{z}_{i\cdot}-\bar{z}_{\cdot\cdot})^{2}}{\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}(z_{ij}-\bar{z}_{i\cdot})^{2}}.$
(3)

If the observations are independently sampled from the same Gaussian
distributions, it is plausible that $\widetilde{F}$ is still approximately
$F$-distributed, because the dependence between the $z_{ij}$ is small even at
moderate group sizes. This was confirmed by simulation in Levene, (1960).
Brown and Forsythe, 1974a  present a larger simulation experiment suggesting
that replacing the $\bar{x}_{i\cdot}$ in the definition of $z_{ij}$ by a
trimmed mean or median leads to a more robust test for non-Gaussian data.

Current best practice suggests to perform a Welch-modified ANOVA directly if
the assumption of equal variance is unclear as it results only in a small loss
of power in the case where the variances are indeed equal. We refer to
Gastwirth et al., (2009) for a comprehensive presentation on Levene’s test
including this question and many further developments.

Levene’s test and its variants remain highly important today as differences in
variances (or some other measure of dispersion) are often in the center of
attention in their own rights. In the rest of the paper we present tests on
differences in “location” of groups and differences in “dispersion” of groups,
both based on inter-point distances in a metric space. Our goal is to combine
one of either kind in order to detect group differences in some universality.

##  3 Non-Euclidean ANOVA

In this and the next sections we assume that our data lies in a general metric
space $(\mathcal{X},d)$. We present existing methods of testing for group
differences based on ANOVA-like ideas. For the presentation we focus on
generalizations of 1-way ANOVA, but provide further information on which
methods can easily be extended to more complex designs. We always assume
having $n=\sum_{i=1}^{k}n_{i}$ independent observations
$x_{ij}\in\mathcal{X}$, $1\leq j\leq n_{i}$, $1\leq i\leq k$ from $k$
potentially different distributions $P_{1},\ldots,P_{k}$ on $\mathcal{X}$
(with Borel $\sigma$-algebra).

###  3.1 Anderson ANOVA

Anderson, (2001) argues, in the context of data sets in ecology, that
traditional multivariate analogues of ANOVA are too stringent in their
assumptions. These are typically based on similar statistics as (1), but with
absolute values replaced by Euclidean norms, see e.g. Mardia et al., (1979)
Section 12.3. We may avoid the use of means of observations by writing
$\mathrm{TSS}-\mathrm{RSS}$ instead of $\mathrm{MSS}$ and replacing the sums
of squared deviations from the mean with the help of the formula

$\displaystyle\sum_{j=1}^{m}\lVert
y_{j}-\overline{y}\rVert^{2}=\frac{1}{2m}\sum_{j_{1},j_{2}=1}^{m}\lVert
y_{j_{1}}-y_{j_{2}}\rVert^{2}=\frac{1}{m}\sum_{j_{1},j_{2}=1}^{m,\,<}\lVert
y_{j_{1}}-y_{j_{2}}\rVert^{2},$

where we indicate by “$<$” in the summation bound that the sum is to be taken
over strictly ordered summands only, here $j_{1}<j_{2}$. Anderson proposes to
replace the pairwise Euclidean distances by more general dissimilarities
between observations and performs a permutation test. In our context we simply
use the pairwise distances in the metric space. Thus

$\displaystyle\mathrm{TSS}$
$\displaystyle=\frac{1}{n}\biggl{(}\sum_{i_{1},i_{2}=1}^{k,\,<}\sum_{j_{1}=1}^{n_{i_{1}}}\sum_{j_{2}=1}^{n_{i_{2}}}d^{2}(x_{i_{1}j_{1}},x_{i_{2}j_{2}})+\sum_{i=1}^{k}\sum_{j_{1},j_{2}=1}^{n_{i},\,<}d^{2}(x_{ij_{1}},x_{ij_{2}})\biggr{)}$
$\displaystyle\mathrm{RSS}$
$\displaystyle=\sum_{i=1}^{k}\frac{1}{n_{i}}\sum_{j_{1},j_{2}=1}^{n_{i},\,<}d^{2}(x_{ij_{1}},x_{ij_{2}})$
$\displaystyle\mathrm{MSS}$ $\displaystyle=\mathrm{TSS}-\mathrm{RSS}$

and the final Anderson ANOVA statistic becomes

$F_{\text{A}}=\frac{n-k}{k-1}\frac{\mathrm{MSS}}{\mathrm{RSS}}.$

It has been noted in various places that this statistic may suffer from type I
error inflation (in terms of a null hypothesis of equal _means_ in Euclidean
space) and substantial loss of power in the unbalanced setting if the groups
are heteroscedastic; see e.g.Alekseyenko, (2016). Anderson et al., (2017) and
Hamidi et al., (2019) propose improvements based on the classical ANOVA
variants by Brown and Forsythe, 1974b  and Welch, (1951), respectively. In the
former, the $F$-statistic is replaced by

$F_{\text{BF}}=\frac{\mathrm{MSS}}{\sum_{i=1}^{k}(1-\frac{n_{i}}{n})\frac{1}{n_{i}(n_{i}-1)}\sum_{j_{1},j_{2}=1}^{n_{i},\,<}d^{2}(x_{ij_{1}},x_{ij_{2}})}.$

For the simulation studies in Section 6 we concentrate on the balanced
setting, for which Anderson $F_{A}$ performs typically well even in presence
of heteroscedacity. We therefore do not discuss these improvements further,
which in the balanced setting do not change the statistic.

###  3.2 Fréchet ANOVA

Dubey and Müller, (2019) introduce ANOVA-like terms that use distances in the
metric $d$ to Fréchet means rather than absolute differences to averages. For
observation $y_{1},\ldots y_{m}\in\mathcal{X}$ the Fréchet mean is defined as

$\displaystyle\bar{y}=\operatorname*{argmin\,}_{z\in\mathcal{X}}\sum_{i=1}^{m}d^{2}(y_{i},z).$
(4)

One of the assumptions in Dubey and Müller, (2019) is that all Fréchet means
considered exist and are unique. For our usual set of observations we denote
by $\bar{x}_{i\cdot}$ the Fréchet mean of $x_{i1},\ldots,x_{in_{i}}$,
$i=1,\ldots,k$ and by $\bar{x}_{\cdot\cdot}$ the Fréchet mean of all
observations. Following the notation in Dubey and Müller, (2019), we write the
Fréchet variance for the $i$-th group and the total Fréchet variance as

$\displaystyle\hat{V}_{i}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}d^{2}(x_{ij},\bar{x}_{i\cdot})\quad\text{and}\quad\hat{V}_{p}=\frac{1}{n}\sum_{i=1}^{k}\sum_{j=1}^{n_{i}}d^{2}(x_{ij},\bar{x}_{\cdot\cdot}),$

respectively. While $\hat{V}_{i}$ is the mean of
$d^{2}(x_{ij},\bar{x}_{i\cdot})$, $j=1,\ldots,n_{i}$, we also require the
corresponding variance

$\displaystyle\hat{\sigma}_{i}^{2}=\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}d^{4}(x_{ij},\bar{x}_{i\cdot})-\hat{V}_{i}^{2}.$

Setting $\lambda_{i}=\frac{n_{i}}{n}$ one finally obtains

$\displaystyle U_{n}$
$\displaystyle=\sum_{i_{1},i_{2}=1}^{k,\,<}\frac{\lambda_{i_{1}}\lambda_{i_{2}}}{\hat{\sigma}_{i_{1}}^{2}\hat{\sigma}_{i_{2}}^{2}}(\hat{V}_{i_{1}}-\hat{V}_{i_{2}})^{2}$
$\displaystyle F_{n}$
$\displaystyle=\hat{V}_{p}-\sum_{i=1}^{k}\lambda_{i}\hat{V}_{i}$
$\displaystyle T$
$\displaystyle=\frac{nU_{n}}{\sum_{i=1}^{k}\frac{\lambda_{i}}{\hat{\sigma}_{i}^{2}}}+\frac{nF_{n}^{2}}{\sum_{i=1}^{k}\lambda_{i}^{2}\hat{\sigma}_{i}^{2}}=:T_{L}+T_{F}.$

In the Euclidean setting of Section 2 the term $F_{n}$ is equal to
$\frac{1}{n}(\mathrm{TSS}-\mathrm{RSS})$ and the denominator of $T_{F}$ is
then an estimator for the variance of $\frac{1}{n}\mathrm{RSS}$, so that
$T_{F}$ has close ties to the ANOVA F-statistic. The unweighted summands
$(\hat{V}_{i_{1}}-\hat{V}_{i_{2}})^{2}$ of $U_{n}$ are similar in spirit to
the terms $(\bar{z}_{i\cdot}-\bar{z}_{\cdot\cdot})^{2}$ from the definition of
Levene’s statistic, and in fact it appears that in the Euclidean case $T_{L}$
corresponds exactly to a simpler variant of Welch’s ANOVA applied to
$d^{2}(x_{ij},\bar{x}_{i\cdot})$, $j=1,\ldots,n_{i}$, $i=1,\ldots,k$; see the
computation in Formulae (8)–(16) in Hamidi et al., (2019). Thus $T_{L}$ has
close ties to Levene’s statistic.

Dubey and Müller show under a list of conditions pertaining to existence and
uniqueness of theoretical and empirical Fréchet means and the complexity of
the metric space (in terms of entropy integrals) that

$\frac{nU_{n}}{\sum_{i=1}^{k}\frac{\lambda_{i}}{\hat{\sigma}_{i}^{2}}}\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\chi^{2}_{k-1}\quad\text{and}\quad\frac{nF_{n}^{2}}{\sum_{i=1}^{k}\lambda_{i}^{2}\hat{\sigma}_{i}^{2}}\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}0\quad\text{as
$n\to\infty$}.$

The authors advocate the simple addition of the two terms in order to obtain a
single test statistic $T$, maybe with weights if there is prior information
available whether to rather look out for inequality of Fréchet means or of
Fréchet variances. However, due to the unbalanced convergence of the two terms
and the fact that the reason for the concrete normalization (especially) of
$T_{F}$ remains a bit inscrutable to us, we prefer to analyze the two summands
separately in Section 6.

##  4 A New Non-Euclidean Method of Levene Type

What appears to be missing is a test for detecting differences of within-group
dispersion that is based directly on pairwise distances between observations
in the metric space. The idea of the PERMDISP-test mentioned in the
introduction, i.e. performing multidimensional scaling and applying Levene’s
test in the principal coordinate space, is to some extent applicable here.
However, it is rather an indirect method and it is methodologically not on the
same level as the Anderson $F_{A}$. Indeed multidimensional scaling can be
applied in combination with _any_ Euclidean procedure, so the PERMDISP-method
should be rather paired up with the analog method of multidimensional scaling
plus applying Euclidean (M)ANOVA. What is more, it contains an unwelcome
tuning parameter, the number of principal coordinates, which is not easy to
choose, but may be crucial. Instead we propose the following test of Levene
type for data in a metric space.

###  4.1 Form and Properties

We assume the same setup as in the previous section, i.e. there are
$n=\sum_{i=1}^{k}n_{i}$ independent observations $x_{ij}\in\mathcal{X}$,
$1\leq j\leq n_{i}$, $1\leq i\leq k$ from $k$ potentially different
distributions $P_{1},\ldots,P_{k}$ on $\mathcal{X}$. Set
$N_{i}=\binom{n_{i}}{2}$ and $N=\sum_{i=1}^{k}N_{i}$. As a surrogate for the
individual deviation terms $z_{ij}$ from Levene’s statistic (3), which in a
general metric space would require the use of a Fréchet or similar mean, we
use $d_{i,\\{j_{1},j_{2}\\}}:=\frac{1}{2}d(x_{ij_{1}},x_{ij_{2}})$. To
simplify the notation, we enumerate the two-element subsets of
$\\{1,\ldots,n_{i}\\}$ by $j=1,\ldots,N_{i}$ and use $d_{ij}$ rather than
$d_{i,\\{j_{1},j_{2}\\}}$ for the $j$-th half-distance in the $i$-th group.

In a first step we assume that $n_{1}=\ldots=n_{k}$ (balanced case) and
emulate the statistics (3) by setting

$\displaystyle L$
$\displaystyle:=\frac{N-k}{k-1}\frac{\sum_{i=1}^{k}n_{i}(\bar{d}_{i\cdot}-\bar{d}_{\cdot\cdot})^{2}}{\sum_{i=1}^{k}\sum_{j=1}^{N_{i}}(d_{ij}-\bar{d}_{i\cdot})^{2}}$
(5)

where

$\displaystyle\bar{d}_{i\cdot}=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}d_{ij}\quad\text{and}\quad\bar{d}_{\cdot\cdot}=\frac{1}{N}\sum_{i=1}^{k}\sum_{j=1}^{N_{i}}d_{ij}$

denote the $i$-th group mean and the overall mean over pairwise distances,
respectively.

Typographically the main fractions of Equations (5) and (3) are very similar,
but the way they use the data $x_{ij}$ is quite different in that we replace
$z_{ij}=\lvert x_{ij}-\bar{x}_{i\cdot}\rvert$, $1\leq j\leq n_{i}$ by
$d_{i,\\{j_{1},j_{2}\\}}=\frac{1}{2}d(x_{ij_{1}},x_{ij_{2}})$, $1\leq
j_{1}<j_{2}\leq n_{i}$. Note that we keep $n_{i}$ in the numerator rather than
replacing it by $N_{i}$, which might have seemed more natural at first glance.
The reason is the substantial dependence of the random variables
$d_{i,\\{j_{1},j_{2}\\}}$ (as opposed to the less substantial dependence
between the $z_{ij}$) for each $i$, which implies that $n_{i}$, not $N_{i}$,
is the correct scaling factor; see Subsection 4.2. Note further that, for the
same reason, the main denominator is not the most natural choice here, but it
is convenient since it keeps the statistic similar to the original Levene
statistic, is fast to compute and empirically performs no worse than the more
natural choice discussed in Subsection 4.2.

There are various ways how one might generalize (5) to general group sizes. We
propose using

$\displaystyle L$
$\displaystyle:=\frac{N-k}{k-1}\frac{\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(\bar{d}_{i\cdot}-\bar{d}_{j\cdot})^{2}}{\sum_{i=1}^{k}\sum_{j=1}^{N_{i}}(d_{ij}-\bar{d}_{i\cdot})^{2}}.$
(6)

Direct computation shows that Equations (6) and (5) agree in the balanced
case, but not in general; see Remark 8. The statistic (6) performs well in
several respects: it allows for an asymptotic distribution ($\chi_{k-1}^{2}$
up to a deterministic factor, see Corollary 3), is still fast to compute and
shows a reasonable performance for unequal group sizes, though it may well be
that a more judicious scaling that takes more proper care of different group
sizes would be superior.

We briefly come back to this last point in Section 6, but do not go much
deeper in the present paper because based on additional considerations, both
theoretical and from simulation studies, we do not see any clear improvements
when choosing different normalizations.

In spite of the limit distribution which we compute in the next section, we
recommend performing a permutation test as for the other statistics
considered. For this we permute the observations, not only their distances,
i.e. new permutations use distances that are potentially different from the
pairwise within-group distances of the original data. As a consequence not
only the RSS changes with permutations, but also the TSS.

It is easy enough to generalize the construction of the above test statistic
to more complex experimental designs. As an example we take up the balanced
two-way ANOVA from Section 2 and form the corresponding Levene-type statistics
for $(\mathcal{X},d)$. For the specific statistics see Section 7.1.

###  4.2 Limit Distribution

In this subsection we derive asymptotic distributions for the statistic $L$
from (6) and for the related statistic

$\displaystyle\widetilde{L}:=\frac{N^{*}-k}{k-1}\frac{\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(\bar{d}_{i\cdot}-\bar{d}_{j\cdot})^{2}}{4\,T_{n}},$
(7)

where $N^{*}=\sum_{i=1}^{k}n_{i}(n_{i}-1)^{2}$ and

$\displaystyle
T_{n}=\sum_{i=1}^{k}\sum_{\begin{subarray}{c}j_{1},j_{2},j_{3}=1\\\
j_{1}\not\in\\{j_{2},j_{3}\\}\end{subarray}}^{n_{i}}\bigl{(}d_{i,\\{j_{1},j_{2}\\}}-\bar{d}_{i\cdot}\bigr{)}\bigl{(}d_{i,\\{j_{1},j_{3}\\}}-\bar{d}_{i\cdot}\bigr{)}.$
(8)

The previous formula makes it necessary to use the more complicated notation
$d_{i,\\{j_{1},j_{2}\\}}=\frac{1}{2}d(x_{ij_{1}},x_{ij_{2}})$ from the
beginning of Subsection 4.1. Note that $\frac{1}{N^{*}-k}T_{n}$ is a natural
group based estimator of
$\mathrm{Cov}\bigl{(}\frac{1}{2}d(X_{1},X_{2}),\frac{1}{2}d(X_{1},X_{3})\bigr{)}$,
where $X_{1},X_{2},X_{3}$ are three independent random variables sampled from
the distribution of the group. The normalization by ${N^{*}-k}$ rather than
${N^{*}}$ is simply modeled after the bias correcting term for independent
data points.

In spite of the ANOVA-like construction, we cannot use the asymptotic theory
for ANOVA directly, because the distances $d_{i,\\{j_{1},j_{2}\\}}$, our
“data”, stem from dependent random variables for each $i$. This dependence is
taken into account by using the factor $\frac{n_{i}n_{j}}{n}$ rather than
$N_{i}$ or $N_{j}$ in the numerator and by normalizing with
$\frac{1}{N^{*}-k}4\,T_{n}$ in (7), which then still allows to obtain the
asymptotic $\chi^{2}_{k-1}$-distribution for $(k-1)\widetilde{L}$. In contrast
$(k-1)L$ converges “only” towards a multiple of $\chi^{2}_{k-1}$ that depends
on parameters of the group distribution.

######  Theorem 2.

Assume that the Borel $\sigma$-algebra for $(\mathcal{X},d)$ is countably
generated. In the usual 1-way setup of Subsection 4.1 assume that
$P_{1}=\ldots=P_{k}=P$ for a distribution $P$ that is not a Dirac distribution
and satisfies
$\int_{\mathcal{X}}\int_{\mathcal{X}}d^{2}(x,y)\,P(dx)\,P(dy)<\infty$. Suppose
that there are $\lambda_{i}>0$ such that $n_{i}/n\to\lambda_{i}$ for every $i$
as $n\to\infty$. Then we have

$(k-1)\,\widetilde{L}\,\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\,\chi^{2}_{k-1}\quad\text{as
$n\to\infty$.}$

######  Corollary 3.

Under the conditions of Theorem 2, we obtain

$(k-1)\,L\,\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\,\frac{4\gamma^{2}}{\sigma^{2}}\chi^{2}_{k-1}\quad\text{as
$n\to\infty$},$

where with independent $X,Y,Z\sim P$ we have

$\displaystyle\gamma^{2}$ $\displaystyle=\mathrm{Cov}(d(X,Y),d(X,Z));$
$\displaystyle\sigma^{2}$ $\displaystyle=\mathrm{Var}(d(X,Y)).$

###### Proof of Theorem 2.

Under the null hypothesis our data is generated by independent
$\mathcal{X}$-valued random elements $X_{ij}\sim P$, $1\leq j\leq n_{i}$,
$1\leq i\leq k$, and the distances $d_{i,\\{j_{1},j_{2}\\}}$ are realizations
of the random variables $\frac{1}{2}d(X_{ij_{1}},X_{ij_{2}})$, $1\leq
j_{1}<j_{2}\leq n_{i}$, $1\leq i\leq k$. Under the conditions on $P$ we have
asymptotic normality of the $U$-statistics

$\displaystyle
U_{i}=U_{i}^{(n)}={\binom{n_{i}}{2}}^{-1}\sum_{j_{1},j_{2}=1}^{n_{i},<}\tfrac{1}{2}d(X_{ij_{1}},X_{ij_{2}}),\quad
i=1,\ldots,k$ (9)

by a straightforward generalization of Hoeffding’s theorem to random elements
in $\mathcal{X}$, see Theorem 5 in the appendix. More precisely, we have with
$X,Y,Z\sim P$ independent that

$\sqrt{n_{i}}\bigl{(}U_{i}-\tfrac{1}{2}\mathbb{E}d(X,Y)\bigr{)}\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\mathcal{N}(0,\gamma^{2})\quad\text{as
$n_{i}\to\infty$},$ (10)

where
$\gamma^{2}=\mathrm{Cov}(d(X,Y),d(X,Z))=\mathrm{Var}(\mathbb{E}(d(X,Y)\,|\,X))=4\gamma_{h}^{2}$
in the notation of the appendix with $h=\frac{1}{2}d$. In view of the 1-way
ANOVA construction, on which $L$ is based, we define the “design matrix”
$D=D_{n}\in\mathbb{R}^{n\times k}$ by

$\displaystyle
D^{\prime}:=\begin{pmatrix}1&\ldots&1&0&\cdots&0&\cdots&0&\cdots&0\\\
0&\cdots&0&1&\ldots&1&\cdots&0&\cdots&0\\\ &\vdots&&&&&\ddots&&\vdots&\\\
0&\cdots&0&0&\cdots&0&\cdots&1&\ldots&1\end{pmatrix}\in\mathbb{R}^{k\times
n},$ (11)

where the $i$-th row has exactly $n_{i}$ ones, and the “contrast matrix”

$\displaystyle C:=\begin{pmatrix}1&0&\cdots&0&-1\\\ 0&1&&0&-1\\\
\vdots&&\ddots&&\vdots\\\ 0&0&&1&-1\end{pmatrix}\in\mathbb{R}^{(k-1)\times
k}.$ (12)

Setting
$\Delta=\lim_{n\to\infty}\frac{1}{n}D_{n}^{\prime}D_{n}=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{n})$,
we obtain with $U=U^{(n)}=(U_{1},\ldots,U_{k})^{\prime}$ by independence of
the components and $n_{i}\to\infty$ as $n\to\infty$ (since $\lambda_{i}>0$)
that

$Z_{n}:=\gamma^{-1}\sqrt{n}\,\Delta^{1/2}\bigl{(}U-\mathbb{E}U\bigr{)}\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\mathcal{N}_{k}(0,I_{k})\quad\text{as
$n\to\infty$}.$ (13)

Setting $\nu=(n_{1},\ldots,n_{k})^{\prime}$, we may further compute
$C^{\prime}(C(D^{\prime}D)^{-1}C^{\prime})^{-1}C=D^{\prime}D-\frac{1}{n}\nu\nu^{\prime}$
(see Lemma 7 in the Appendix for the calculation) and therefore

$\widetilde{L}=\frac{N^{*}-k}{k-1}\frac{U^{\prime}C^{\prime}(C(D_{n}^{\prime}D_{n})^{-1}C^{\prime})^{-1}CU}{4\,T_{n}}.$
(14)

Since
$\mathbb{E}U=\frac{1}{2}\mathbb{E}d(X,Y)\cdot\mathds{1}\in\mathbb{R}^{k}$ and
$C\cdot\mathds{1}=0$, we obtain

$(k-1)\widetilde{L}=\gamma^{2}\frac{Z_{n}^{\prime}\bigl{(}\tfrac{1}{n}W_{n}\bigr{)}Z_{n}}{\frac{4}{N^{*}-k}T_{n}},$

where
$W_{n}:=\Delta^{-1/2}C^{\prime}(C(D_{n}^{\prime}D_{n})^{-1}C^{\prime})^{-1}C\Delta^{-1/2}$.
Note that

$W:=\lim_{n\to\infty}\frac{1}{n}W_{n}=\Delta^{-1/2}C^{\prime}(C\Delta^{-1}C^{\prime})^{-1}C\Delta^{-1/2}$

is a symmetric and idempotent matrix of rank $k-1$, and therefore
$Z^{\prime}WZ\sim\chi^{2}_{k-1}$ for $Z\sim\mathcal{N}_{k}(0,I_{k})$ by Lemma
9 from the Appendix. Using (13) it is straightforward to show with the help of
the continuous mapping theorem that

$Z_{n}^{\prime}\bigl{(}\tfrac{1}{n}W_{n}\bigr{)}Z_{n}\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\chi^{2}_{k-1}.$

So it suffices to show that $\frac{1}{N^{*}-k}T_{n}\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\gamma_{d/2}^{2}$. For this we note that the normalized
inner sum of (8) satisfies

$\displaystyle\frac{1}{n_{i}(n_{i}-1)^{2}}\sum_{\begin{subarray}{c}j_{1},j_{2},j_{3}=1\\\
j_{1}\not\in\\{j_{2},j_{3}\\}\end{subarray}}^{n_{i}}\bigl{(}d_{i,\\{j_{1},j_{2}\\}}-\bar{d}_{i\cdot}\bigr{)}\bigl{(}d_{i,\\{j_{1},j_{3}\\}}-\bar{d}_{i\cdot}\bigr{)}$
$\displaystyle\hskip
17.07164pt=\underbrace{\frac{n_{i}(n_{i}-1)(n_{i}-2)}{n_{i}(n_{i}-1)^{2}}}_{\longrightarrow\,1}\underbrace{\frac{1}{n_{i}(n_{i}-1)(n_{i}-2)}\sum_{j_{1},j_{2},j_{3}=1}^{n_{i},\,\neq}\bigl{(}d_{i,\\{j_{1},j_{2}\\}}-\bar{d}_{i\cdot}\bigr{)}\bigl{(}d_{i,\\{j_{1},j_{3}\\}}-\bar{d}_{i\cdot}\bigr{)}}_{\longrightarrow\,\mathrm{Cov}(\frac{1}{2}d(X_{1},X_{2}),\frac{1}{2}d(X_{1},X_{3}))\,=\,\gamma^{2}_{d/2}}$
(15) $\displaystyle\hskip
34.1433pt{}+\underbrace{\frac{1}{n_{i}-1}}_{\longrightarrow\,0}\underbrace{\frac{1}{n_{i}(n_{i}-1)}\sum_{j_{1},j_{2}}^{n_{i},\,\neq}\bigl{(}d_{i,\\{j_{1},j_{2}\\}}-\bar{d}_{i\cdot}\bigr{)}^{2}}_{\longrightarrow\,\mathrm{Var}(\frac{1}{2}d(X_{1},X_{2}))=\sigma^{2}_{d/2}},$

where convergence of the averages is almost surely and follows after expansion
of the products by the strong law of large numbers for $U$-statistics using
the prerequisite $\mathbb{E}(d(X_{1},X_{2})^{2})<\infty$; see Hoeffding,
(1961).

Thus for the total term

$\displaystyle\frac{1}{N^{*}-k}T_{n}$
$\displaystyle=\frac{1}{N^{*}-k}\sum_{i=1}^{k}n_{i}(n_{i}-1)^{2}\cdot\frac{1}{n_{i}(n_{i}-1)^{2}}\sum_{\begin{subarray}{c}j_{1},j_{2},j_{3}=1\\\
j_{1}\not\in\\{j_{2},j_{3}\\}\end{subarray}}^{n_{i}}\bigl{(}d_{i,\\{j_{1},j_{2}\\}}-\bar{d}_{i\cdot}\bigr{)}\bigl{(}d_{i,\\{j_{1},j_{3}\\}}-\bar{d}_{i\cdot}\bigr{)}\longrightarrow\gamma^{2}_{d/2}.$

∎

###### Proof of Corollary 3.

This follows from Theorem 2 because

$L=4\frac{\frac{1}{N^{*}-k}T_{n}}{\frac{1}{N-k}\sum_{i=1}^{k}\sum_{j=1}^{N_{i}}(d_{ij}-\bar{d}_{i\cdot})^{2}}\widetilde{L},$

where the numerator is a consistent estimator of $\gamma^{2}/4$ and the
denominator is a consistent estimator of $\sigma^{2}/4$, see (15). ∎

##  5 Metric Space of Finite Point Patterns

In Sections 6 and 7 we apply the four statistics from Table 1 for the space of
finite point patterns equipped with the metric introduced in Müller et al.,
(2020). For self-containedness we give a short summary of the relevant
concepts and results, referring to the paper as MSM20.

For $n\in\mathbb{Z}_{+}$ write $[n]=\\{1,2,\ldots,n\\}$ (including
$[0]=\emptyset$). Denote by $\mathfrak{N}_{\mathrm{fin}}$ the space of finite
multisets on a complete separable metric space $(\mathcal{R},\varrho)$. We
refer to the elements
$\xi=\\{x_{1},x_{2},\ldots,x_{n}\\}\in\mathfrak{N}_{\mathrm{fin}}$ as point
patterns, where $n\in\mathbb{Z}_{+}=\\{0,1,2,\ldots\\}$ and
$x_{i}\in\mathcal{X}$, $i\in[n]$. Note that $x_{i}=x_{j}$ for $i\neq j$ is
allowed and that the point patterns can be identified with the counting
measure $\sum_{i=1}^{n}\delta_{x_{i}}$, which is often helpful for theoretical
considerations. We write $\lvert\xi\rvert$ to denote the total number of
points in the pattern $\xi$.

######  Definition 4 (Definition 1 of MSM20).

Let $C>0$ and $p\geq 1$ be two parameters, referred to as _penalty_ and
_order_ , respectively.

For
$\xi=\\{x_{1},\ldots,x_{m}\\},\eta=\\{y_{1},\ldots,y_{n}\\}\in\mathfrak{N}_{\mathrm{fin}}$
define the _transport-transform (TT) metric_ by

$\begin{split}d_{\mathrm{TT}}(\xi,\eta)=d_{\mathrm{TT}}^{(C,p)}(\xi,\eta)=\biggl{(}\min\biggl{(}(m+n-2l)C^{p}+\sum_{r=1}^{l}\varrho(x_{i_{r}},y_{j_{r}})^{p}\biggr{)}\biggr{)}^{1/p},\end{split}$
(16)

where the minimum is taken over equal numbers of pairwise different indices
$i_{1},\ldots,i_{l}$ in $[m]$ and $j_{1},\ldots,j_{l}$ in $[n]$, i.e. over the
set

$\begin{split}S(m,n)=\bigl{\\{}(i_{1},&\ldots,i_{l};\hskip
1.5ptj_{1},\ldots,j_{l})\,;\;l\in\\{0,1,\ldots,\min\\{m,n\\}\\},\\\
&i_{1},\ldots,i_{l}\in[m]\text{ pairwise
different},\,j_{1},\ldots,j_{l}\in[n]\text{ pairwise
different}\bigr{\\}}.\end{split}$

The distance $d_{\mathrm{TT}}(\xi,\eta)$ can be computed by filling up the
smaller point pattern with dummy points located at distance $C$ until it has
the same cardinality $n$ as the larger point pattern and then solving a
standard assignment problem with cost $\min\\{d(x,y),2^{1/p}C\\}$ between
points $x,y$ (MSM20, Theorem 1). The classical worst-case complexity of this
is $O(n^{3})$ (MSM20, Remark 1), which can be somewhat improved to order
$n^{2.5}$ up to polylogarithmic factors (Lee and Sidford,, 2014). Practical
computation times for well over $n=1000$ points are less than one second (R
package ttbary, Müller and Schuhmacher,, 2021, using the auction algorithm
from Bertsekas,, 1988).

The TT-metric can be interpreted as an unbalanced Wasserstein metric (Remark
3). Computing Fréchet means in Wasserstein spaces is a topic of active
research; see e.g. Borgwardt and Patterson, (2020), Borgwardt and Patterson,
(2021), Heinemann et al., (2021) and references therein for recent
developments the space of discrete measures. In our context an additional
increase in difficulty comes from the constraint that the result must be a
discrete measure with integer cardinality. In MSM20 we therefore apply an
alternating heuristics to obtain local minima of the Fréchet functional in
(4). The resulting “pseudo-barycenters” are obtained much faster and appear to
be of good quality (consistent objective function values and results conform
with intuition), but are by no means perfect and still require considerable
computation time for hundreds of patterns with hundred of points (Table 1–4 in
MSM20).

A related metric that we take up in Section 7 is the relative TT-metric
defined as

$d_{\mathrm{RTT}}(\xi,\eta)=d_{\mathrm{RTT}}^{(C,p)}(\xi,\eta)=\frac{1}{\max\\{\lvert\xi\rvert,\lvert\eta\rvert\\}^{1/p}}d_{\mathrm{TT}}^{(C,p)}(\xi,\eta).$
(17)

This metric is in a sense more robust to individual outliers if there are many
points. In particular note that $d_{\mathrm{RTT}}(\xi_{N},\xi_{N}\cup\zeta)\to
0$ as $N\to\infty$ if $\lvert\xi_{N}\rvert\to\infty$ and $\zeta$ is a fixed
point pattern.

In view of the conditions for Theorem 2, completeness and separability are
inherited from $(\mathcal{X},\varrho)$ to
$(\mathfrak{N}_{\mathrm{fin}},d_{\mathrm{TT}})$ and
$(\mathfrak{N}_{\mathrm{fin}},d_{\mathrm{RTT}})$. This is straightforward to
see after checking that $d_{\mathrm{RTT}}(\xi_{N},\xi)\to 0$ iff
$d_{\mathrm{TT}}(\xi_{N},\xi)\to 0$ iff
$\lvert\xi_{N}\rvert\to\lvert\xi\rvert$ and each point $x$ of $\xi$ is
approximated by exactly one point of $\xi_{N}$ (if $x$ is a multipoint of
cardinality $k$ this means that there is a total of exactly $k$ points in
$\xi_{N}$, possibly forming multipoints of their own, that converge towards
$x$). The condition
$\int_{\mathcal{X}}\int_{\mathcal{X}}d(x,y)\,P(dx)\,P(dy)<\infty$ is always
satisfied for $d_{\mathrm{RTT}}$ because it is bounded by $C$. Since
$d_{\mathrm{TT}}(\xi,\eta)\leq
C\max\\{\lvert\xi\rvert,\lvert\eta\rvert\\}^{1/p}$ it is satisfied for
$d_{\mathrm{TT}}$ if $\Xi\sim P$ satisfies
$\mathbb{E}\lvert\Xi\rvert^{2/p}<\infty$, which is for example the case for
all point process distributions considered in Section 6.

For the simulation study in the next section it is helpful to understand some
basic probability measures on $\mathfrak{N}_{\mathrm{fin}}$. Suppose that
$\mathcal{R}\subset\mathbb{R}^{d}$ is compact (in the next section we only use
a unit square in $\mathbb{R}^{2}$). A random element in the metric space
$(\mathfrak{N}_{\mathrm{fin}},d_{\mathrm{RTT}})$, equipped with its Borel
$\sigma$-algebra is called a _point process_ , i.e. a point process is a
measurable map from a probability space to $\mathfrak{N}_{\mathrm{fin}}$. The
Borel $\sigma$-algebra coincides with the smallest $\sigma$-algebra that makes
$\xi\mapsto\xi(A)$ measurable for every measurable $A\subset\mathcal{R}$,
which is the usual $\sigma$-algebra considered on
$\mathfrak{N}_{\mathrm{fin}}$; see Proposition 9.1.IV in Daley and Vere-Jones,
(2008).

We say a point process $\Xi$ satisfies _complete spatial randomness (CSR)_ if
it is a Poisson process with intensity measure
$\nu=\lambda\,\mathrm{Leb}^{d}$, where $\lambda\geq 0$ and $\mathrm{Leb}^{d}$
is Lebesgue measure (on $\mathcal{R}$). This means that
$\Xi(A)\sim\mathrm{Po}(\nu(A))$ for all measurable $A\subset\mathcal{R}$ and
that $\Xi(A_{1}),\ldots,\Xi(A_{l})$ are independent for all $l\in\mathbb{N}$
and all measurable $A_{1},\ldots,A_{l}\subset\mathcal{R}$ that are pairwise
disjoint. See e.g. Section 2.4 in Daley and Vere-Jones, (2003) for more
details on the Poisson process.

##  6 Simulation Study

We tested the different statistics from Table 1 for various point process
distributions and present the results in what follows. First we investigate
the practical use of our asymptotics in Subsection 6.1. In spatial statistics
there are usually two fundamentally different ways how distributions can
deviate from CSR. One is spatial inhomogeneity of points, i.e. points may be
more or less likely to occur in different regions of the space. The ability of
tests to detect deviations from CSR against various spatially inhomogeneous
alternatives is studied in Subsection 6.2. The other way is interaction of
points, i.e. presence of points in one region may excite or inhibit the
presence of other points nearby. In Subsection 6.3 we study how well the
statistics discern between various interaction strengths in homogeneous
Strauss processes.

For the evaluations in Subsections 6.2 and 6.3 we perform permutation tests.
These are based on generating $M$ independent uniform permutations of the
indices of the data points resulting in alternative split-ups of the data into
$k$ groups of sizes $n_{i}$, $1\leq i\leq k$. We then determine the rank $r$
of the statistic-value for the original split-up within the statistic-values
of the alternative split up (from $r=1$ for the highest value to $r=M+1$ for
the lowest value). It is easily checked (and well-known) that
$p=\frac{r}{M+1}$ is an honest p-value (i.e.
$\mathbb{P}(p\leq\alpha)\leq\alpha$ for every $\alpha\in(0,1)$). We reject the
null if $p\leq 0.05$.

In Subsections 6.2 and 6.3 we have $k=2$, $n_{i}=\tilde{n}=20$ and use $M=999$
permutations if no barycenter computation is needed and $M=99$ permutations if
barycenter computation is needed. In view of the $\binom{40}{20}\approx
1.4\cdot 10^{11}$ possible split-ups, this means that there is a high degree
of randomization in each individual test. The small $M=99$ was necessary due
to the large computational burden of computing pseudo-barycenters in point
pattern space (see Section 5). For statistics that do not require barycenter
computation, choosing $M=999$ typically results in much faster computation
time than the choice of $M=99$ for statistics that do require barycenter
computation. For reproducibility of individual test results, a higher $M$ or
(where possible) comparing within all possible split-ups into groups would be
desirable in both cases.

Preferring exact permutation tests over tests based on the limit
$\chi^{2}$-distribution is in agreement with the recommendations from previous
papers and corresponds to our own experience. However, the
$\chi^{2}$-approximation of our $L$-statistic is quite fast as we can see in
Subsection 6.1, where we compare the finite sample distributions of the new
$L$\- and the Dubey–Müller statistics.

In all tests we use as the underlying space
$\mathcal{R}=[0,1]^{2}\subset\mathbb{R}^{2}$ with the Euclidean metric. The
significance level is always $\alpha=0.05$. Furthermore we choose as order
$p=2$ and as penalty $C=0.25$, which means that $\sqrt{2}\cdot 0.25\approx
0.35$ is the maximal contribution that a single matched point pair makes to
the TT-distance, i.e. the actual Euclidean distances are cut off at this
value. In applications the choice of $C$ is often based on the physical
reality of the data and possibly the goal of the analysis. For the present
simulation study we tried not to restrict a substantial proportion of matching
distances while keeping the contribution of additional points reasonably low.
Table 2 gives an overview of how many pairs are matched above and below the
cutoff distance for various values of $C$ based on pairwise comparisons of
1000 point patterns simulated according to CSR with intensity $\lambda=35$.
For $C=0.25$ we have for every matching above the cutoff distance
$1/0.038\approx 26$ matchings below the cutoff distance.

|  mean $d_{\mathrm{TT}}$ |  below cutoff |  above cutoff |  unpaired
---|---|---|---|---
$C=0.1$ |  0.309 |  11123388 |  4469722 |  3309250
relative |  |  1 |  0.424 |  0.311
$C=0.15$ |  0.393 |  13456877 |  2136233 |  3309250
relative |  |  1 |  0.167 |  0.257
$C=0.2$ |  0.457 |  14514420 |  1078690 |  3309250
relative |  |  1 |  0.078 |  0.24
$C=0.25$ |  0.512 |  15054688 |  538422 |  3309250
relative |  |  1 |  0.038 |  0.233
$C=0.3$ |  0.561 |  15347529 |  245581 |  3309250
relative |  |  1 |  0.017 |  0.23
$C=0.35$ |  0.609 |  15498175 |  94935 |  3309250
relative |  |  1 |  0.006 |  0.229
Table 2:  Pairwise comparison within $1000$ patterns simulated from CSR on
$[0,1]^{2}$ with intensity $\lambda=35$ for various penalties $C$. The columns
give the mean $d_{\mathrm{TT}}$-distance, and the (absolute and relative)
number of matchings below cutoff and above cutoff, as well as the number of
unpaired points for these $\binom{1000}{2}$ pairwise comparisons. Note that
the relative numbers are with respect to the number of matchings below cutoff.

###  6.1 Asymptotics

In the present subsection we numerically assess the speed of convergence of
our new $\widetilde{L}$ statistic under the null hypothesis of equal group
distributions towards the $\chi^{2}_{k-1}$ distribution as presented in
Subsection 4.2. For comparison we also consider the Fréchet $T_{L}$ and $T$
statistics, which were shown in Dubey and Müller, (2019) to have a limiting
$\chi^{2}_{k-1}$ distribution as well.

Our experiments are based on $k=2$ groups, both simulated from the same
distribution, which is either $\mathrm{CSR}(35)$ or the Strauss hard core
distribution with $\lambda=35$. These are the extreme distributions having
either no interaction or very strong interaction in Subsection 6.3. As group
size we consider $\tilde{n}=5,20,50,200$. The computation of the Fréchet $T$
and $T_{L}$ depend on the calculation of a barycenter. For this we used the
heuristic algorithm presented in Müller et al., (2020). The calculation of an
exact barycenter is computationally infeasible for this kind of data. To
compensate that we do not get the optimal solution, we did $5$ restarts in
every barycenter calculation and used the best of the $5$ solutions as the
barycenter.

Figure 1 shows QQ-plots for the empirical distributions of our new Levene
statistic $\widetilde{L}$, the Fréchet statistic $T_{L}$ and the Fréchet
statistic $T$ on the $y$-axis and the theoretical $\chi_{1}^{2}$ distribution
on the $x$-axis. The data are the $\mathrm{CSR}(35)$ point patterns. In the
first column the groups consist of $\tilde{n}=5$ patterns, in the second
column of $\tilde{n}=20$ patterns and so on. For the two Levene statistics
$\widetilde{L}$ and $T_{L}$ we can see the computed quantiles approach the
theoretical quantiles as the group size $\tilde{n}$ gets larger. Even for a
medium group size $\tilde{n}=50$ the computed quantiles are very close to the
theoretical quantiles of a $\chi_{1}^{2}$ distribution.

Similarly Figure 2 shows QQ-plots for hardcore Strauss distributed point
patterns. Again the four columns correspond to the four group sizes
$\tilde{n}=5,20,50,200$ and the three rows correspond to the three statistics.
For this data the computed quantiles are already very close to the theoretical
quantiles of a $\chi_{1}^{2}$ distribution for $\tilde{n}=20$ for the two
Levene statistics.

In both cases the third row, the combined Fréchet statistic $T$, yields
quantiles that are far from the theoretical quantiles. This is solely due to
the summand $T_{F}$ that is not considered in the second row.

In spite of the asymptotic results we use permutation based tests in what
follows. This is in the tradition of previous methods, see e.g. Anderson,
(2017), Dubey and Müller, (2019). Comparison of different statistics for
different data sets are presented in the following two subsections.

Figure 1: QQ-plots of the percentiles based on 500 statistics values (on the
$y$-axis) versus $\chi^{2}_{1}$-percentiles. Based on $k=2$ groups of
$\tilde{n}=5,20,50,200$ patterns from $\mathrm{CSR}(35)$. The first row is our
new $\widetilde{L}$ statistic (7), the second and third rows are the Fréchet
$T_{L}$ statistic from Section 3.2 and the Fréchet $T$ statistic,
respectively.  Figure 2: QQ-plots of the percentiles based on 500 statistics
values (on the $y$-axis) versus $\chi^{2}_{1}$-percentiles. Based on $k=2$
groups of $\tilde{n}=5,20,50,200$ patterns from a Strauss hard core
distribution with $\lambda=35$. The first row is our new $\widetilde{L}$
statistic (7), the second and third rows are the Fréchet $T_{L}$ statistic
from Section 3.2 and the Fréchet $T$ statistic, respectively.

###  6.2 Inhomogeneity

Here we compare $k=2$ groups of $\tilde{n}=20$ point patterns. Patterns in
Group 2 are simulated from CSR with $\lambda=35$. In Group 1, they are
simulated from various inhomogeneous scenarios, i.e. from Poisson process
distributions where the intensity function (the density of the measure $\nu$
with respect to Lebesgue measure) deviates more or less from a constant but
still integrates up to $35$ over the whole window $\mathcal{R}=[0,1]^{2}$.

In Scenarios 1–3 the intensity is obtained by adding a number of rotation-
invariant Gaussian distributions with different means but the same covariance
matrix $\sigma^{2}I$ and scaling to total mass $35$. For simplicity we do not
restrict the intensity to $\mathcal{R}$, but as can be seen from Figure 3 only
very few points outside $\mathcal{R}$ occur. Scenarios 4–6 use as intensity an
exponential function that is constant in the $y$-coordinate and induces a
certain tendency for points to lie in the left part of the window rather than
in the right part.

Table 3 provides more information about the chosen parameters. Figure 3 shows
five example point patterns for each scenario. In addition we add a Scenario
0, which corresponds to simulating the first group also from CSR with
$\lambda=35$.

Scenario |  $\lambda(x,y)$ proportional to
---|---
$1$ | $\sum_{i=1}^{3}\varphi_{\mu_{i},0.075}(x,y)$
$2$ | $\sum_{i=1}^{3}\varphi_{\mu_{i},0.1}(x,y)$
$3$ | $\sum_{i=1}^{4}\varphi_{\mu_{i},0.1}(x,y)$
$4$ | $\exp(-2x)$
$5$ | $\exp(-1x)$
$6$ | $\exp(-0.02x)$
Table 3:  Overview of the Poisson process intensities for the six scenarios.
The proportionality constant is chosen such that the expected number of points
in each scenario is $35$. By $\varphi_{\mu,\sigma^{2}}$ we denote the density
of the bivariate normal distribution with mean $\mu\in\mathbb{R}^{2}$ and
covariance matrix $\sigma^{2}I$. The different $\mu_{i}$ used are seen in
Figure 3.

Figure 3:  Five example patterns for each of the six scenarios together with
the window $[0,1]^{2}$ (dotted line) in which the homogeneous Poissons
processes of the second group are sampled.

Table 4 gives the results in terms of numbers of rejections (out of 100) of
the null hypothesis of equal distribution in both groups.

Scenario |  $1$ |  $2$ |  $3$ |  $4$ |  $5$ |  $6$ |  $0$
---|---|---|---|---|---|---|---
Anderson $F_{A}$ |  100 |  100 |  100 |  100 |  99 |  39 |  2
new $L$ |  93 |  76 |  77 |  14 |  7 |  9 |  3
Fréchet $T_{F}$ |  100 |  100 |  100 |  99 |  11 |  0 |  4
Fréchet $T_{L}$ |  59 |  24 |  47 |  13 |  9 |  12 |  4
Table 4:  Numbers of rejections of the null hypothesis “equal distribution in
both groups” based on 100 data sets per column. In each data set the first
group is sampled from the scenario indicated in the column and the second
group is sampled from Scenario 0.

We observe that the direct ANOVA procedures perform much better than the
Levene (or indirect ANOVA) procedures. This is not so surprising, because the
inhomogeneity experiment considers two groups of distributions that are
different in terms of their location in the point pattern space. To see this
intuitively, think about the distributions in Scenarios 1–6 (and 0 as a
boundary case) in terms of producing locally perturbed versions of a typical
point pattern, which is more or less any one of the example point patterns in
Figure 3 (more appropriately one would rather think of an idealized version of
these patterns, such as the Fréchet mean). Among the direct ANOVA methods,
Anderson $F_{A}$ performs substantially better than Fréchet $T_{F}$ and has
still a reasonable chance to detect the faint differences between Scenarios 6
and 0 when presented with the 20 patterns from each group. Our new L-test
performs somewhat better than the Fréchet L-test, but both tests are only able
to detect the inhomogeneity (with reasonable probability) when it is very
obvious (Scenarios 1–3).

###  6.3 Interaction between Points

Again we compare $k=2$ groups of $\tilde{n}=20$ point patterns. This time the
group distributions differ in the degree of point interaction. For this we
consider the distribution of the homogeneous Strauss process on the unit
square $\mathcal{R}=[0,1]^{2}$, which is obtained by specifying the density
$f\colon\mathfrak{N}_{\mathrm{fin}}\to\mathbb{R}_{+}$,

$f(\xi):=c\cdot\beta^{\lvert\xi\rvert}\cdot\gamma^{s_{R}(\xi)},$

with respect to CSR with intensity 1 on $\mathcal{R}$, where

$s_{R}(\xi)=\sum_{\\{x,y\\}\subset\xi}\mathds{1}\\{\lVert x-y\rVert\leq R\\}$

is the number of pairs of points at distance $\leq R$ from one another. Here
$R>0$ is the range of the interaction, $\beta>0$ is the so-called activity
(which controls the intensity of the process via an increasing function, that
is however only accessible numerically) and $\gamma\in[0,1]$ is the strength
of the interaction. The constant $c$ normalizes the density to an overall
integral of $1$ and is also not available in closed form. We write
$\mathrm{Strauss}(\beta,\gamma;R)$ for this point process distribution.
Intuitively a $\mathrm{Strauss}(\beta,\gamma;R)$-process is obtained from a
$\mathrm{CSR}(\beta)$ process by penalizing each outcome according to a factor
$\gamma$ per $R$-close point pair. Correspondingly we have
$\mathrm{Strauss}(\beta,1;R)=\mathrm{CSR}(\beta)$ (regardless of $R$). At the
other end of the spectrum $\mathrm{Strauss}(\beta,0;R)$ is the distribution of
a hard core process with no points allowed within distance $R$ of other
points.

For the simulation we set $R=0.1$ and consider scenarios based on the six
different values $\gamma=0,0.2,0.4,0.6,0.8,1$. The activity $\beta$ is adapted
so that each time $\lambda=35$. Figure 4 shows one realization for each of the
six scenarios.

Figure 4: Simulations from $\mathrm{Strauss}(\beta,\gamma;0.1)$-distributions,
where rowwise from left to right $\gamma=0,0.2,0.4,0.6,0.8,1$ and $\beta$ is
adjusted such that $\lambda=35$. For $\gamma=0$ we have a realization of a
hard core process, for $\gamma=1$ a realization from CSR.

We perform two different experiments here. In the first one the patterns in
Group 1 are sampled from $\mathrm{CSR}(35)$ corresponding to $\gamma=1$, in
the second one they are sampled from the mentioned hard core process with
$\lambda=35$ corresponding to $\gamma=0$. The patterns in Group 2 are sampled
in both experiment from each of the six $\gamma$-values in turn. The results
are listed in Table 5.

$\gamma=1$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Anderson $F_{A}$ |  $98$ |  $41$ |  $13$ |  $8$ |  $4$ |  $5$
new $L$ |  $100$ |  $100$ |  $95$ |  $67$ |  $20$ |  $3$
Fréchet $T_{F}$ |  $100$ |  $98$ |  $78$ |  $45$ |  $20$ |  $4$
Fréchet $T_{L}$ |  $100$ |  $99$ |  $88$ |  $45$ |  $20$ |  $4$

$\gamma=0$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Anderson $F_{A}$ |  $3$ |  $55$ |  $90$ |  $98$ |  $97$ |  $99$
new $L$ |  $11$ |  $60$ |  $96$ |  $100$ |  $100$ |  $100$
Fréchet $T_{F}$ |  $6$ |  $57$ |  $91$ |  $97$ |  $100$ |  $100$
Fréchet $T_{L}$ |  $9$ |  $33$ |  $82$ |  $95$ |  $99$ |  $100$

Table 5: Numbers of rejections of the null hypothesis “equal distribution in
both groups” based on 100 data sets per column. In each data the point
patterns in both groups are sampled from a Strauss distribution with
$\lambda=35$ and $R=0.1$. The first group is sampled using $\gamma=1$ or
$\gamma=0$ as indicated on the top left of the table and the second group uses
$\gamma$ as indicated in the column.

In contrast to the situation in the previous subsection (different
inhomogeneity), we now observe that the _indirect_ ANOVA procedures, i.e. the
Levene-type tests perform considerably better than the direct ANOVA
procedures. Again this is intuitively understandable because a small $\gamma$
in the Strauss process leads to less dispersion, both in terms of a smaller
variance for the total number of points and also with respect to typical
distances of points from one another: for small $\gamma$ the points are quite
regularly placed, whereas for larger $\gamma$ there are erratic patches that
are free of points leading typically to some points that have to be matched
over longer distances, which in the squared Euclidean metric has quite some
influence. A small $\gamma$ will also lead to smaller average distances than a
larger $\gamma$ (either between point patterns or relative to a barycenter),
which may explain why the difference in the performance of the indirect and
direct ANOVA tests is somewhat less pronounced than in the inhomogeneity
experiment.
Note again that the powers of the tests based on pairwise distances are
slightly better than those of the tests based on barycenters.

##  7 Applications

In this section we apply our Levene’s test to a real data example. We
investigate the location of bubbles in a mineral flotation experiment. The
structure of the data calls for a two factor design. We establish a distance
based two-way Levene’s test and compare its performance to existing methods.
The classical two-way ANOVA design can be found in Section 2.

###  7.1 Balanced Two-Way Levene’s Test

As mentioned in Subsection 4.1 it is easy to generalize statistic (5) to a
two-way design, that will further be useful for the bubble data analyzed in
the next Subsection.

Suppose we have independent observations $x_{i_{1}i_{2}j}\in\mathcal{X}$,
$1\leq j\leq\tilde{n}$, $1\leq i_{1}\leq k_{1}$, $1\leq i_{2}\leq k_{2}$ from
groups obtained by crossing a Factor $a$ with $k_{1}$ levels and a Factor $b$
with $k_{2}$ levels with $\tilde{n}$ observations for each combination. In a
similar way as above we denote by $d_{i_{1}i_{2}j}$ the $j$-th half-distance
in the group $(i_{1},i_{2})$, where
$j=1,\ldots,\tilde{N}:=\binom{\tilde{n}}{2}$. Set then

$\displaystyle\mathrm{RSS}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\sum_{j=1}^{\tilde{N}}(d_{i_{1}i_{2}j}-\bar{d}_{i_{1}i_{2}\cdot})^{2}$
$\displaystyle\mathrm{MSS}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\tilde{n}(\bar{d}_{i_{1}i_{2}\cdot}-\bar{d}_{\cdot\cdot\cdot})^{2}$
$\displaystyle\mathrm{SSa}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}k_{2}\tilde{n}(\bar{d}_{i_{1}\cdot\cdot}-\bar{d}_{\cdot\cdot\cdot})^{2}$
$\displaystyle\mathrm{SSb}$
$\displaystyle=\sum_{i_{2}=1}^{k_{2}}k_{1}\tilde{n}(\bar{d}_{\cdot
i_{2}\cdot}-\bar{d}_{\cdot\cdot\cdot})^{2}$ $\displaystyle\mathrm{SSi}$
$\displaystyle=\sum_{i_{1}=1}^{k_{1}}\sum_{i_{2}=1}^{k_{2}}\tilde{n}(\bar{d}_{i_{1}i_{2}\cdot}-\bar{d}_{i_{1}\cdot\cdot}-\bar{d}_{\cdot
i_{2}\cdot}+\bar{d}_{\cdot\cdot\cdot})^{2},$

where the various means are taken over the dot components in the usual way.
Note that we never use any distances between observations of different factor
combinations.

In addition to the omnibus test for group differences as in one-way ANOVA, we
may then perform Levene-type tests for effects of Factor a and b separately,
as well as for an interaction effect. The corresponding statistics are

$\displaystyle
L=\frac{N-k_{1}k_{2}}{(k_{1}k_{2}-1)}\frac{\mathrm{MSS}}{\mathrm{RSS}},\hskip
5.69054ptLa=\frac{N-k_{1}k_{2}}{k_{1}-1}\frac{\mathrm{SSa}}{\mathrm{RSS}},\hskip
5.69054ptLb=\frac{N-k_{1}k_{2}}{k_{2}-1}\frac{\mathrm{SSb}}{\mathrm{RSS}},\hskip
5.69054ptLi=\frac{N-k_{1}k_{2}}{(k_{1}-1)(k_{2}-1)}\frac{\mathrm{SSi}}{\mathrm{RSS}}.$

###  7.2 Bubble Data

Figure 5: Arrangement of floating bubbles data. Rows represent the three
frother concentration levels and columns the three volumetric air flowrate
levels (treatments). Each cell contains six spatial point patterns
(responses).

We consider the data from González et al., (2021) which provides locations of
bubbles in a mineral flotation experiment, where the interest is analysing if
the spatial distribution might be affected by frother concentrations and
volumetric airflow rates. Indeed, the data set consists of 54 images
containing a total of 8385 floating bubbles. The images of bubbles can be
regarded as spatial point patterns where the centroids of the bubbles
correspond to the points. In addition, we have three frother concentration
levels (5 ppm, 10 ppm, 15 ppm) as well as three volumetric airflow rate levels
(5 l/min, 8 l/min, 10 l/min), and we have six replicates of point patterns at
each combination of levels of such factors. The treatment combinations of the
experiment, as well as the observed bubble point patterns, are represented in
Figure 5.

We used the two-way design of Levene’s statistic from Section 7.1 to test for
influence of the different factors, interaction and differences between the
groups. For comparison we also used the two factor statistics from Anderson,
(2001), we performed a two factor ANOVA on the number of points per pattern,
and finally complemented our analysis with a two factor ANOVA with
$K$-functions, so as to link our analysis with that of González et al.,
(2021). We did a permutation test with $999$ permutations.

In Section 6, the cutoff was always fixed to $C=0.25$. This was a reasonable
value for point patterns with expected $35$ points in the unit square. In the
bubble data, the number of points per observed pattern ranges from $21$ to
$353$. With such a great variability in the number of points we suggest
adjusting the cutoff to prevent that distances between two patterns are
dominated by their different numbers of points. For the results presented in
this section we computed the mean number of points of the tested patterns
$\bar{n}$ and used the cutoff $\bar{C}=0.25\cdot 35/\bar{n}$ for the
computations of $d_{TT}$. For more details to the cutoff see (16).

$p$-values | FC | VA | Interaction | Overall
---|---|---|---|---
Anderson $F_{A}$ | 0.003 | 0.001 | 0.001 | 0.001
new $L$ | 0.001 | 0.001 | 0.001 | 0.001
Number of points | 0.001 | 0.001 | 0.001 | 0.001
$K$-functions | 0.005 | 0.001 | 0.001 | 0.001∗∗

∗∗ this is the p-value for the sum of both factors, not the overall ANOVA
statistic.

Table 6: Results for the different tests for the bubble data. Quantiles are
obtained by a permutation test with $999$ permutations. The cutoff is
$C=0.0564$, the maximal radius for the $K$-functions is $r=0.15$.

$p$-values | FC | VA | Interaction | Overall
---|---|---|---|---
Anderson $F_{A}$ | 0.043 | 0.001 | 0.019 | 0.001
new $L$ | 0.001 | 0.001 | 0.001 | 0.001
Number of points | 0.001 | 0.002 | 0.001 | 0.001
$K$-functions | 0.002 | 0.022 | 0.002 | 0.006∗∗

∗∗ this is the p-value for the sum of both factors, not the overall ANOVA
statistic.

Table 7: Results for the different tests for the bubble data, leaving out the
frother concentration of $15$ppm. Quantiles are obtained by a permutation test
with $999$ permutations. The cutoff is $C=0.0636$, the maximal radius for the
$K$-functions is $r=0.15$.

The p-values of the permutation tests are shown in Tables 6 and 7. In
particular, Table 6 shows results for the whole data set, while Table 7
depicts results for only part of the data, leaving out the third column, i.e.
any patterns from frother concentration of $15$ ppm. In both cases, our new
Levene, Anderson $F_{A}$, the ANOVA on the number of points per pattern, and
the ANOVA for $K$-functions detect significant influence of each of the two
factors and the interaction. We already recommended to always perform both,
the tests for differences in variability and the test for differences of
means. In the second test scenario, both Levene’s test and Anderson $F_{A}$
detect significance for the frother concentration and the interaction of both
parameters for our usual significance level of $5\%$. But the relative
difference between the $p$-values of the two tests is very large. For the
smaller significance level of $1\%$, our Levene’s test still detects
significance where Anderson $F_{A}$ does not. So the test for differences of
means might not be enough in a practical application. This is particularly
important in cases where, as it is the case for the bubble data, the number of
points plays a crucial role in the behavior and structure of the point
patterns.

We see that for this data apparently the numbers of points per pattern contain
enough information to detect significant influence of the factors. This is not
very surprising since the number of points per pattern is similar in the $6$
patterns of a single cell, but the differences between cells are large.
This observation is reinforced by a classical multidimensional scaling (mds).
Based on the TT-distances between the point patterns, we translated every
point pattern into a single point in $\mathbb{R}^{2}$. The mds was applied
first for the whole bubble data set, see Figure 6, and then for a subset of
the data consisting of the first and second columns, leaving out the data with
a frother concentration of $15$ ppm, see Figure 7. This is the same data that
we used for our analyses in Tables 6 and 7. The three levels of the air flow
are encoded by the colors ‘red’, ‘green’ and ‘blue’, same color means same air
flow rate, and the three levels of the frother concentration are encoded by
the symbols ‘circle’, ‘triangle’ and ‘cross’. When we compare these plots to
the images of the point patterns in Figure 5 we can see that the
multidimensional scaling sorts the point patterns from left to right in
ascending order by their number of points per pattern. In Figure 7 we can see
that the points that correspond to the data with a frother concentration of
$5$ ppm, i.e. the circles, and the data with a frother concentration of $10$
ppm, i.e. the triangles, are scattered differently. The (coordinate-wise)
means of the triangles and circles are similar, but we can see that the
circles are more scattered along both axes. We conjecture that it is this
difference in scatter that our Levene’s test is able to detect in Table 7,
whereas the Anderson $F_{A}$ only barely detects a slight difference in means.

Figure 6: The bubble data after a multidimensional scaling into two dimensions
based on the distance matrix w.r.t the TT-metric. Figure 7: The bubble data
without the third column, i.e. without data where FC=$15$ppm, after a
multidimensional scaling into two dimensions based on the distance matrix
w.r.t the TT-metric.

##  8 Conclusions and Discussion

In this paper we gave an overview of some ANOVA procedures that can be used
for data in general metric spaces. We introduced a new method that is similar
to Levene’s test and compared it to existing methods with regard to point
pattern data in Section 6. In the studies, see Tables 4 and 5 for the results,
we compared the distance-based ANOVA from Anderson, (2001), the distance-based
Levene’s test, (5), introduced in this paper and the tests based on the ANOVA
statistic $T_{F}$ and the Levene statistic $T_{L}$ from Dubey and Müller,
(2019).
The latter proposed in their paper the combined statistic $T=T_{L}+T_{F}$. In
our simulations we wanted to put a focus on the two fundamentally different
ways of “location” and “dispersion” in which group distributions can differ,
even in an abstract metric space. We therefore did our tests with the
statistics $T_{L}$ and $T_{F}$ seperately, which allows us to have a direct
comparison between the distance-based statistics and the statistics of Dubey
and Müller, (2019). We also did the tests with the proposed combined statistic
$T$, and for completeness we give the results in Tables 8 and 9. In the tests
for differences in interaction, comparing Tables 9 and 5, we can see that the
performance of the statistic $T$ is “between” the performance of $T_{L}$ and
$T_{F}$. For the tests of inhomogeneity, see Tables 8 and 4, the combined
statistic $T$ performs almost as good as the better statistic of $T_{L}$ and
$T_{F}$, which is in this case $T_{F}$.
For the presented scenarios and the chosen parameters all the statistics
worked well for their designated purpose, in particular the ANOVA statistics
for detecting inhomogeneity and the Levene’s statistics for detecting
differences in point interaction. But for different scenarios this might not
be the case. For a cutoff of $C=0.1$ instead of the proposed $C=0.25$, we
observed that the statistic $T$ and the Anderson $F_{A}$ do not work as well
anymore in detecting differences in interaction, see Tables 10 and 11. The
performance of $T$ is considerably worse than the performance of $T_{L}$ and
$T_{F}$ is working even more poorly as well. The distance based ANOVA
statistic of Anderson also performs very poorly compared to the cutoff of
$C=0.25$, while our statistic $L$ is working even better with the smaller
cutoff.

For future research it would be interesting to take a closer look at the
(co)variance estimator $\gamma$ of our $\widetilde{L}$ statistic. This
estimator is not unbiased and it remains open if a statistic with an unbiased
estimator works even better, in particular for the asymptotics.
There are also more complex designs for the $2$-factor ANOVA, with different
sums of squares or designs that allow for different group sizes. Our $L$
statistic could be generalized to more complex $2$-factor designs, or even
$k$-factor designs.
Additionally it would be interesting to test our statistics with more kinds of
data. On the one hand there are different kinds of point pattern data, e.g.
marked point patterns, but also data from other metric spaces, e.g. image data
or graph data.

In Section 6 we already mentioned that the computation of the Fréchet $T$
statistic is more time consuming than the distance based tests, because of the
barycenter calculation. If we take for example the data from Section 6.3, the
distance based Anderson $F_{A}$ and our new $L$ take about $8$ seconds and $2$
seconds, respectively, for $100$ permutation tests with $999$ permutations
each. The calculation of $100$ permutation tests of Fréchet $T$ with $99$
permutations each takes about $45$ minutes. The workhorse computation of all
tests is done in C++. However, parts of the overhead for Anderson $F_{A}$ and
Fréchet $T$ are programmed in R and a complete implementation in C++ might
improve the runtime a little bit. These numbers are merely meant to give a
general impression of the order of magnitude of the runtimes. The distance
based tests take only a few seconds, while for the Fréchet tests the
computation of the barycenters takes many minutes, even with the fast
heuristics instead of the exact solution and with fewer permutations.

Overall we find that the new $L$ in combination with the Anderson $F_{A}$ has
a similar performance and allows for considerably faster computation than the
other methods in settings where the computation of barycenters is costly.

Scenario |  $1$ |  $2$ |  $3$ |  $4$ |  $5$ |  $6$ |  $0$
---|---|---|---|---|---|---|---
Fréchet $T$ |  100 |  100 |  100 |  98 |  11 |  2 |  6
Table 8: Performance of the sum statistic Fréchet $T$. Numbers of rejections
of the null hypothesis “equal distribution in both groups” based on 100 data
sets per column for the $7$ scenarios of inhomogeneity.

$\gamma=1$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Fréchet $T$ |  $100$ |  $98$ |  $77$ |  $48$ |  $16$ |  $2$

$\gamma=0$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Fréchet $T$ |  $5$ |  $48$ |  $91$ |  $98$ |  $100$ |  $100$

Table 9: Performance of the sum statistic Fréchet $T$. Numbers of rejections of the null hypothesis “equal distribution in both groups” based on 100 data sets per column for the different scenarios of interaction between points. $\gamma=1$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Anderson $F_{A}$ |  37 |  16 |  8 |  5 |  12 |  3
new $L$ |  100 |  100 |  96 |  69 |  15 |  6
Fréchet $T_{F}$ |  49 |  31 |  10 |  4 |  5 |  9
Fréchet $T_{L}$ |  100 |  99 |  87 |  47 |  13 |  7
Fréchet $T$ |  99 |  84 |  38 |  12 |  7 |  9
Table 10: C=0.1, The first scenario: Groups of Poisson-distributed point patterns vs groups of Strauss-distributed point patterns with $6$ different gammas: $\gamma=0,0.2,0.4,0.6,0.8,1$, $R=0.1$, $\lambda=35$, $\alpha=0.05$, $20$ patterns per group. Numbers indicate how many times the hypothesis ”equal distributions in both groups” is rejected out of $100$ times. The tests should see no difference between groups of Poisson-patterns and Strauss-patterns with $\gamma=1$. $\gamma=0$ vs. |  $\gamma=0$ |    $\gamma=0.2$ |    $\gamma=0.4$ |    $\gamma=0.6$ |    $\gamma=0.8$ |  $\gamma=1$
---|---|---|---|---|---|---
Anderson $F_{A}$ |  5 |  11 |  35 |  41 |  40 |  50
new $L$ |  5 |  99 |  100 |  100 |  100 |  100
Fréchet $T_{F}$ |  8 |  14 |  26 |  41 |  40 |  53
Fréchet $T_{L}$ |  7 |  87 |  100 |  100 |  100 |  100
Fréchet $T$ |  7 |  31 |  73 |  95 |  100 |  100
Table 11: C=0.1, The second scenario: Groups of Strauss-distributed point
patterns with a fixed $\gamma=0$ vs groups of Strauss-distributed point
patterns with $6$ different gammas: $\gamma=0,0.2,0.4,0.6,0.8,1$, $R=0.1$,
$\lambda=35$, $\alpha=0.05$, $20$ patterns per group. Numbers indicate how
many times the hypothesis ”equal distributions in both groups” is rejected out
of $100$ times. The tests should see no difference for $\gamma=0$.

## References

  * Alekseyenko, (2016) Alekseyenko, A. V. (2016).  Multivariate Welch t-test on distances.  Bioinformatics, 32(23):3552–3558.
  * Anderson, (2001) Anderson, M. J. (2001).  A new method for non-parametric multivariate analysis of variance.  Austral ecology, 26(1):32–46.
  * Anderson, (2006) Anderson, M. J. (2006).  Distance-based tests for homogeneity of multivariate dispersions.  Biometrics, 62(1):245–253.
  * Anderson, (2017) Anderson, M. J. (2017).  Permutational multivariate analysis of variance (PERMANOVA).  Wiley statsref: statistics reference online, pages 1–15.
  * Anderson et al., (2017) Anderson, M. J., Walsh, D. C., Robert Clarke, K., Gorley, R. N., and Guerra-Castro, E. (2017).  Some solutions to the multivariate Behrens–Fisher problem for dissimilarity-based analyses.  Australian & New Zealand Journal of Statistics, 59(1):57–79.
  * Bertsekas, (1988) Bertsekas, D. P. (1988).  The auction algorithm: A distributed relaxation method for the assignment problem.  Annals of operations research, 14(1):105–123.
  * Billera et al., (2001) Billera, L. J., Holmes, S. P., and Vogtmann, K. (2001).  Geometry of the space of phylogenetic trees.  Advances in Applied Mathematics, 27(4):733–767.
  * Borgwardt and Patterson, (2020) Borgwardt, S. and Patterson, S. (2020).  Improved linear programs for discrete barycenters.  Informs Journal on Optimization, 2(1):14–33.
  * Borgwardt and Patterson, (2021) Borgwardt, S. and Patterson, S. (2021).  On the computational complexity of finding a sparse Wasserstein barycenter.  Journal of Combinatorial Optimization, 41(3):736–761.
  * (10) Brown, M. B. and Forsythe, A. B. (1974a).  Robust tests for the equality of variances.  Journal of the American Statistical Association, 69(346):364–367.
  * (11) Brown, M. B. and Forsythe, A. B. (1974b).  The small sample behavior of some statistics which test the equality of several means.  Technometrics, 16(1):129–132.
  * Cuevas et al., (2004) Cuevas, A., Febrero, M., and Fraiman, R. (2004).  An ANOVA test for functional data.  Computational statistics & data analysis, 47(1):111–122.
  * Daley and Vere-Jones, (2003) Daley, D. and Vere-Jones, D. (2003).  An introduction to the theory of point processes. Vol. I.  Springer, 2nd edition.  Elementary theory and methods.
  * Daley and Vere-Jones, (2008) Daley, D. and Vere-Jones, D. (2008).  An introduction to the theory of point processes. Vol. II.  Springer, 2nd edition.  General theory and structure.
  * Denker and Keller, (1983) Denker, M. and Keller, G. (1983).  On U-statistics and v. Mises’ statistics for weakly dependent processes.  Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete, 64(4):505–522.
  * Dubey and Müller, (2019) Dubey, P. and Müller, H.-G. (2019).  Fréchet analysis of variance for random objects.  Biometrika, 106(4):803–821.
  * Fisher, (1925) Fisher, R. (1925).  Statistical Methods for Research Workers.  Oliver & Boyd.
  * Gastwirth et al., (2009) Gastwirth, J. L., Gel, Y. R., and Miao, W. (2009).  The impact of Levene’s test of equality of variances on statistical theory and practice.  Statistical Science, 24(3):343–360.
  * Ginestet et al., (2017) Ginestet, C. E., Li, J., Balachandran, P., Rosenberg, S., and Kolaczyk, E. D. (2017).  Hypothesis testing for network data in functional neuroimaging.  The Annals of Applied Statistics, pages 725–750.
  * González et al., (2021) González, J. A., Lagos-Álvarez, B. M., and Mateu, J. (2021).  Two-way layout factorial experiments of spatial point pattern responses in mineral flotation.  TEST, pages 1–30.
  * Hamidi et al., (2019) Hamidi, B., Wallace, K., Vasu, C., and Alekseyenko, A. V. (2019).  $W_{d}^{*}$-test: robust distance-based multivariate analysis of variance.  Microbiome, 7(1):1–9.
  * Heinemann et al., (2021) Heinemann, F., Munk, A., and Zemel, Y. (2021).  Randomised Wasserstein barycenter computation: Resampling with statistical guarantees.  Preprint.  Available at https://arxiv.org/abs/2012.06397.
  * Hoeffding, (1948) Hoeffding, W. (1948).  A class of statistics with asymptotically normal distribution.  The Annals of Mathematical Statistics, 19(3):293–325.
  * Hoeffding, (1961) Hoeffding, W. (1961).  The strong law of large numbers for U-statistics.  Technical Report, Mimeograph Series 302, Department of Statistics, University of North Carolina.
  * Huckemann et al., (2009) Huckemann, S., Hotz, T., and Munk, A. (2009).  Intrinsic MANOVA for Riemannian manifolds with an application to Kendall’s space of planar shapes.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(4):593–603.
  * Lee and Sidford, (2014) Lee, Y. T. and Sidford, A. (2014).  Path finding methods for linear programming: Solving linear programs in Õ(vrank) iterations and faster algorithms for maximum flow.  In 2014 IEEE 55th Annual Symposium on Foundations of Computer Science, pages 424–433. IEEE.
  * Levene, (1960) Levene, H. (1960).  Robust tests for equality of variances.  Contributions to Probability and Statistics, pages 278–292.
  * Mardia et al., (1979) Mardia, K., Kent, J., and Bibby, J. (1979).  Multivariate Analysis.  Academic Press, London.
  * Müller and Schuhmacher, (2021) Müller, R. and Schuhmacher, D. (2021).  ttbary: Barycenter Methods for Spatial Point Patterns.  R package version 0.2-0. https://CRAN.R-project.org/package=ttbary.
  * Müller et al., (2020) Müller, R., Schuhmacher, D., and Mateu, J. (2020).  Metrics and barycenters for point pattern data.  Statistics and Computing, 30:953–972.
  * Ramón et al., (2016) Ramón, P., de la Cruz, M., Chacón-Labella, J., and Escudero, A. (2016).  A new non-parametric method for analyzing replicated point patterns in ecology.  Ecography, 39(11):1109–1117.
  * Scheffé, (1967) Scheffé, H. (1967).  The analysis of variance.  John Wiley & Sons, 5th printing, 1st edition.
  * Welch, (1951) Welch, B. L. (1951).  On the comparison of several mean values: An alternative approach.  Biometrika, 38(3/4):330–336.
  * Wooldridge, (2010) Wooldridge, J. M. (2010).  Econometric Analysis of Cross Section and Panel Data.  MIT press.

##  Appendix A Auxiliary Results Used for the Proof of Theorem 2

For completeness and self-containedness we state here (consequences of)
results from the literature as well as some additional calculations needed for
the proof of Theorem 2.

Firstly we formulate a straightforward generalization of Hoeffding’s theorem
for the asymptotic normality of $U$-statistics (univariate version of Theorem
7.1 in Hoeffding,, 1948) for random elements in the general metric space
$\mathcal{X}$ with countably generated Borel $\sigma$-algebra. See also
Theorem 1(b) of Denker and Keller, (1983), where this result is further
generalized to (weakly) dependent sequences of random elements.

######  Theorem 5.

Let $(X_{n})_{n\in\mathbb{N}}$ be an i.i.d. sequence of $\mathcal{X}$-valued
random elements. Let $h\colon\mathcal{X}^{m}\to\mathbb{R}$ be symmetric and
non-degenerate in the sense that there are $x_{2},\ldots,x_{m}\in\mathcal{X}$
such that

$\mathbb{E}h(X_{1},x_{2},\ldots,x_{m})\neq 0.$

Suppose further that
$\mathbb{E}\bigl{(}h(X_{1},\ldots,X_{m})^{2}\bigr{)}<\infty$. We write

$U_{n}={\binom{n}{m}}^{-1}\sum_{\begin{subarray}{c}i_{1},\ldots,i_{m}=1\\\
i_{1}<\ldots<i_{m}\end{subarray}}^{n}h(X_{i_{1}},\ldots,X_{i_{m}}).$

for the $U$-statistic with kernel $h$. Then

$\sqrt{n}(U_{n}-\mathbb{E}(U_{n}))\stackrel{{\scriptstyle\mathcal{D}}}{{\longrightarrow}}\mathcal{N}(0,m^{2}\gamma_{h}^{2}),$

where for an independent copy $(\tilde{X}_{2},\ldots,\tilde{X}_{m})$ of
$(X_{2},\ldots,X_{m})$

$\gamma_{h}^{2}=\mathrm{Cov}\bigl{(}h(X_{1},X_{2},\ldots,X_{m}),h(X_{1},\tilde{X}_{2},\ldots,\tilde{X}_{m})\bigr{)}=\mathrm{Var}\bigl{(}\mathbb{E}(h(X_{1},\ldots,X_{m})\,|\,X_{1})\bigr{)}.$

######  Remark 6.

In the setting of Theorem 5 above, Theorem 5.2 of Hoeffding, (1948) yields

$m^{2}\gamma_{h}^{2}\leq n\mathrm{Var}(U_{n})\leq
m\mathrm{Var}(h(X_{1},\ldots,X_{m}))$

for all $n\geq m$. The right hand bound is sharp for $n=m$ and
$n\mathrm{Var}(U_{n})\searrow m^{2}\gamma_{h}^{2}$ as $n\to\infty$.

The above inequality means in particular that for finite $n$ the expression
$\frac{m^{2}}{n}\gamma_{h}^{2}$ can only underestimate $\mathrm{Var}(U_{n})$.
The exact formula for $m=2$ is

$n\mathrm{Var}(U_{n})=\frac{n-2}{n-1}\cdot 4\gamma_{h}^{2}+\frac{1}{n-1}\cdot
2\mathrm{Var}(h(X_{1},X_{2})).$

The next result is similar to classical ANOVA. For completeness we give its
proof.

######  Lemma 7.

Let $C\in\mathbb{R}^{(k-1)\times k}$ as in (12), $D\in\mathbb{R}^{n\times k}$
as in (11), $U=(u_{1},\ldots,u_{k})^{\prime}$ and
$\nu=(n_{1},\ldots,n_{k})^{\prime}$. We have

$C^{\prime}(C(D^{\prime}D)^{-1}C^{\prime})^{-1}C=D^{\prime}D-\frac{1}{n}\nu\nu^{\prime}$

and

$U^{\prime}(D^{\prime}D-\frac{1}{n}\nu\nu^{\prime})U=\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(u_{i}-u_{j})^{2}$

###### Proof.

Define

$\displaystyle\nu_{(i)}$ $\displaystyle:=(n_{1},\ldots,n_{i})^{\prime}\text{ ,
}\;\Lambda_{(i)}:=\mathrm{diag}(\nu_{(i)})\in\mathbb{R}^{i\times i}\;\text{
and }\;\mathds{1}_{(i)}:=(1,\ldots,1)^{\prime}\in\mathbb{R}^{i}.$

Then $\mathds{1}_{(i)}\mathds{1}_{(i)}^{\prime}$ is the $i\times i$ matrix of
$1$’s. We build up the equality step by step. Since
$D^{\prime}D=\Lambda_{(k)}$ and therefore

$\displaystyle(D^{\prime}D)^{-1}=(\Lambda_{(k)})^{-1}=\mathrm{diag}(1/n_{1},\ldots,1/n_{k}),$

We obtain

$\displaystyle C(D^{\prime}D)^{-1}C^{\prime}$
$\displaystyle=(\Lambda_{(k-1)})^{-1}+\frac{1}{n_{k}}\cdot\mathds{1}_{(k-1)}\mathds{1}_{(k-1)}^{\prime}$

and

$\displaystyle(C(D^{\prime}D)^{-1}C^{\prime})^{-1}$
$\displaystyle=\Lambda_{(k-1)}-\frac{1}{n}\cdot\nu_{(k-1)}\nu_{(k-1)}^{\prime}$

and finally

$\displaystyle
C^{\prime}(C(D^{\prime}D)^{-1}C^{\prime})^{-1}C=\Lambda_{(k)}-\frac{1}{n}\cdot\nu\nu^{\prime}$

When we multiply the vector $U$ from left and right, the $ij$-th entry in the
matrix is the coefficient of $u_{i}u_{j}$. This leads to

$\displaystyle U^{\prime}(D^{\prime}D$
$\displaystyle-\frac{1}{n}\nu\nu^{\prime})U$
$\displaystyle=\sum_{i=1}^{k}n_{i}u_{i}^{2}-\frac{1}{n}\sum_{i=1}^{k}n_{i}^{2}u_{i}^{2}-\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}2n_{i}n_{j}u_{i}u_{j}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{k}n_{i}u_{i}^{2}\sum_{j=1}^{k}n_{j}-\frac{1}{2n}\sum_{i=1}^{k}n_{i}^{2}(u_{i}^{2}+u_{i}^{2})-\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}2n_{i}n_{j}u_{i}u_{j}$
$\displaystyle=\frac{1}{2n}\sum_{i=1}^{k}\sum_{j=1}^{k}n_{i}n_{j}(u_{i}^{2}+u_{j}^{2})-\frac{1}{2n}\sum_{i=1}^{k}n_{i}^{2}(u_{i}^{2}+u_{i}^{2})-\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}2n_{i}n_{j}u_{i}u_{j}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(u_{i}^{2}+u_{j}^{2})-\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}2n_{i}n_{j}u_{i}u_{j}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(u_{i}-u_{j})^{2}.$

∎

######  Remark 8.

Let $\bar{u}=\frac{1}{k}\sum_{i=1}^{k}u_{i}$. An equivalent expression for
$U^{\prime}(D^{\prime}D-\frac{1}{n}\nu\nu^{\prime})U$ in Lemma 7 can be
computed as

$\displaystyle\frac{1}{n}$
$\displaystyle\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}n_{i}n_{j}(u_{i}-u_{j})^{2}$
$\displaystyle=\frac{1}{2n}\sum_{i=1}^{k}\sum_{j=1}^{k}n_{i}n_{j}((u_{i}-\bar{u})+(\bar{u}-u_{j}))^{2}$
$\displaystyle=\frac{1}{2n}\sum_{i=1}^{k}\sum_{j=1}^{k}n_{i}n_{j}(u_{i}-\bar{u})^{2}+\frac{1}{2n}\sum_{i=1}^{k}\sum_{j=1}^{k}n_{i}n_{j}(\bar{u}-u_{j})^{2}+\frac{1}{2n}\sum_{i=1}^{k}\sum_{j=1}^{k}n_{i}n_{j}(u_{i}-\bar{u})(\bar{u}-u_{j})$
$\displaystyle=\sum_{i=1}^{k}n_{i}(u_{i}-\bar{u})^{2}+\frac{1}{2n}\sum_{i=1}^{k}n_{i}(u_{i}-\bar{u})\sum_{j=1}^{k}n_{j}(\bar{u}-u_{j})$
$\displaystyle=\sum_{i=1}^{k}n_{i}(u_{i}-\bar{u})^{2}+\frac{1}{2n}\left(\sum_{i=1}^{k}n_{i}(u_{i}-\bar{u})\right)^{2}$

If $n_{1}=\ldots=n_{k}=\tilde{n}$, we see directly from the right-hand side
that

$\displaystyle\frac{1}{n}\sum_{i=1}^{k-1}\sum_{j=i+1}^{k}$ $\displaystyle
n_{i}n_{j}(u_{i}-u_{j})^{2}=\sum_{i=1}^{k}n_{i}(u_{i}-\bar{u})^{2}=\tilde{n}\sum_{i=1}^{k}(u_{i}-\bar{u})^{2}.$
(18)

The following lemma is well known. It follows by spectral decomposition, see
e.g. Kent, Mardia and Bibby (1979), Theorem 3.4.4(b), setting $p=1$ and
$\Sigma=I$.

######  Lemma 9.

Let $Z\sim\mathcal{N}_{n}(0,I)$ and let $C\in\mathbb{R}^{n\times n}$ be
symmetric and idempotent. Then $Z^{\prime}CZ\sim\chi^{2}_{r}$, where
$r=\mathrm{trace}(C)=\mathrm{rank}(C)$.