# Linguistics from a topological viewpoint

Rui Dong<EMAIL_ADDRESS>Centre for Linguistics, Leiden University,
Leiden, the Netherlands

###### Abstract

Typological databases in linguistics are usually categorical-valued. As a
result, it is difficult to have a clear visualization of the data. In this
paper, we describe a workflow to analyze the topological shapes of South
American languages by applying multiple correspondence analysis technique and
topological data analysis methods.

Multiple correspondence analysis, Topological data analysis, Quechuan,
Nuclear-Macro-Jê, South American languages

###### Contents

  1. I Introduction
  2. II Related work
  3. III Multiple correspondence analysis
  4. IV Topological data analysis
    1. IV.1 A toy example of pretzel
    2. IV.2 The TDA workflow
  5. V Data analysis procedure
    1. V.1 Data preprocessing
    2. V.2 MCA processing
    3. V.3 The shapes of languages
  6. VI Applications of TDA
    1. VI.1 Nuclear-Macro-Jê
    2. VI.2 Quechuan
  7. VII Discussion
  8. A Multiple correspondence analysis
  9. B Background about topological data analysis
    1. B.1 Simplices and complexes
    2. B.2 Chains and boundary operators
    3. B.3 Basics about persistent homology
    4. B.4 Permutation test

##  I Introduction

Typological data collected by linguists is usually categorical-valued. For
instance, in the recently released database Grambank
grambank_dataset_zenodo_v1  $189$ out of the $195$ features are binary, and
the rest $6$ features are ternary. In the Grambank dataset, each language is
represented by multiple categorical values, as a result, it is difficult to
measure the difference between two languages. For instance, in the case of the
commonly used Gower distance, although it is simple and concise by definition,
the Gower distance does not reflect the frequency of each value in each
feature, which is a crucial characteristic held by a categorical-valued
datasets.

In contrast to this, in this paper we first perform dimensional reduction via
the multiple correspondence analysis (MCA) method, the advantage of the MCA
method is that the position of a feature value encodes the information of its
frequency inside the dataset, so that for each language we can form a sub-
point cloud consisting of the values of all features. We then apply the
framework of topological data analysis (TDA) to analyze the topological
structure of the distributions of each point cloud. Technically speaking we
can use TDA to analyze the topological invariants in any general dimension,
however in practice the meaning of higher dimensional topological invariants
is usually vague. Hence our analysis focuses only on the $1$-dimensional
topological invariants in each sub-point cloud, that is, the circular
structures.

In this paper we restrict ourselves to South American languages. Specifically,
we apply the topological methods to analyze the Nuclear-Macro-Jê (NMJ) family
and the Quechuan family. In sub-Section VI.1, we show that, within NMJ, there
is a significant distinction between languages of the Jê-proper and the non-
Jê-proper languages, and in sub-Section VI.2 we show that, within the Quechuan
family, there is a significant distinction between northern and southern
Quechuan languages.

##  II Related work

The topological method has already been applied into linguistics, as is
presented in MR4368966 ; MR3767893 , where the PCA method is applied to the
Syntactic Structures of the World’s Languages (SSWL) database
111https://ling.yale.edu/syntactic-structures-worlds-language-cross-
linguistic-database. In these previous works, each language is represented as
a point inside a Euclidean space, and the TDA method is applied to analyze the
topological structures of different language families.

Our method is inspired by these approaches. However, in contrast to these
previous works, in our paper we apply MCA method to display a point cloud for
each single language, so that it is possible to form a visualization for each
single language and to compare the shapes of different languages.

##  III Multiple correspondence analysis

The Multiple correspondence analysis (MCA)MCA  is a dimension reduction
technique that is analogous to principal component analysis (PCA) for dealing
with categorical-valued features instead of quantitative ones. One main
difference between MCA and PCA is that in the case of MCA, one can obtain
either a set of points representing all the samples or a set of points
representing all the categorical values. In our context, we apply MCA to
project all the categorical values into a Euclidean space $\mathbb{R}^{d}$.
The Appendix A contains the details of how to implement the MCA in practice.

There are two noteworthy characteristics in MCA:

  * •

The distance of a point away from the original point relies on the frequency
of the corresponding categorical values. The lower frequency categorical
values appear further away from the original point in $\mathbb{R}^{d}$.

  * •

The direction of a point depends on the intrinsic linguistical meaning. In the
case of binary features, the two categorical values of the same feature lie on
a line through the original point of $\mathbb{R}^{d}$, and these two points
are in the opposite directions apart from the original point.

Hence in the point cloud of all categorical values, if a point is far away
from the original point, it implies that the representative value is special
in the dataset, the languages that share this specific value should be more
relevant. On the other hand, if a point is close to the original point, it
means that the corresponding categorical value is commonly appeared in the
dataset, hence in some way we hope to put less weight on it.

This phenomenon suggests that, if we extract the corresponding values of a
specific language out of this point cloud, we can then obtain a sub-point
cloud, the “shape” of which should exhibit some linguistical meaning of this
chosen language. As we have discussed above, for a chosen language, if the
value of a feature is common, it will be close to the original point, and if
the value is special, it will be far away from the original point, the
direction of which relies on the intrinsic meaning from linguistical
viewpoint.

More than obtaining only an intuitive visualization of a language, we would
like to be able to investigate the topological structures of the corresponding
sub-point clouds. This is exactly what the framework of topological data
analysis is meant to do.

##  IV Topological data analysis

The framework of topological data analysis (TDA) normally starts from a point
cloud $P$ inside a Euclidean space $\mathbb{R}^{d}$. We can regard the TDA as
a higher dimensional generalization of clustering analysis. That means, more
than detecting cluster structures of a point cloud, the TDA method is able to
detect higher dimensional topological structures (e.g., holes, voids, etc.) of
the point cloud $P$. The Appendix B contains the basic technical definitions
and results necessary for this context. To get a more intuitive understanding
of the main idea of TDA, we start from a simple example of a pretzel-shaped
point cloud.

###  IV.1 A toy example of pretzel

Consider the following pretzel-shaped point cloud in Figure 1. Intuitively,
this point cloud forms three loops.

Figure 1: A pretzel-shaped point cloud

However, since the points are discrete, the topology of it is trivial, we need
some more technical method to describe the shape of this point cloud.

One way to deal with this issue is to cover each point with a disk of radius
$r$, as the radius $r$ ranges from $0$ to infinity, the union of all disks
form a connected topological space.

We display the persistent behaviors of circular structures in Figure 2, from
which we can observe that as the radius increases, three circular structures
are born (b), then two of them disappear (c), finally all disks get
intersected together and the third circular structure is dead (d). This is the
core idea of the TDA framework, as the radius of disks increases, we can
detect the persistent topological structures (such as circular structures) of
the point cloud via investigating the union of all the disks.

However, the union of all the disks (the green objects in Figure 2) is not
easy to deal with, to get around this problem, we construct a “skeleton” (the
purple objects in Figure 2) out of the union of disks. In terms of TDA, this
“skeleton” is referred to as “Vietoris-Rips simplicial complex”MR2572029 . One
fundamental result in TDA is that the union of all the disks and the
“skeleton” have the same circular structures (this is referred to as the
“Nerve theorem”MR2572029  in the framework of TDA.). Hence, we can work on the
“skeleton” instead of the union of all disks.

###  IV.2 The TDA workflow

Before describing the framework of TDA, we first need to recall a topological
object called “$n$-simplex”. Intuitively speaking, an $n$-simplex is a
generalized “$n$-dimensional triangle”. For instance, a $0$-simplex is a
point, a $1$-simplex is an edge, a $2$-simplex is simply a triangle, a
$3$-simplex is a tetrahedron. From these examples, we observe that one point
determines a $0$-simplex, two points determine a $1$-simplex, three points
determine a $2$-simplex, four points determine a $3$-simplex, and in a similar
way, a general $k$-simplex can be defined as an $k$-dimensional object
determined by $k+1$ points. We denote the $k$-simplex determined by the set of
points $\\{p_{0},p_{1},...,p_{k}\\}$ as $\langle
p_{0},p_{1},...,p_{k}\rangle$.

We can then construct multiple simplices out of a point cloud $P$ by the
pairwise distances of the points. In more detail, we fix a threshold $r>0$, if
the pairwise distances of any $k+1$ points $\\{p_{0},p_{1},...,p_{k}\\}\in P$
are all less than the threshold $r$, then these $k+1$ points determine a
$k$-simplex. The collection of all such simplices determined by the threshold
$r$ is referred to as the Vietoris-Rips simplicial complex with respect to the
point cloud $P$, which is denoted as $\operatorname{\mathrm{Rips}}_{r}(P)$.

Figure 2: From (a) to (d): the radius ranges from $0$, the circular structures
were born, remained, and dead.

The object of $\operatorname{\mathrm{Rips}}_{r}(P)$ is still too abstract to
deal with, however, we can construct a collection of vector spaces out of
$\operatorname{\mathrm{Rips}}_{r}(P)$ over the field
$\mathbf{F_{2}}=\\{0,1\\}$, we refer to these vector spaces as spaces of
$k$-chains, and denote them as $C_{k}$ for $k\geq 0$. In more detail, a
$k$-chain $\gamma$ is a formal finite sum that looks like:
$\gamma=c_{1}\sigma_{1}+c_{2}\sigma_{2}+...+c_{n}\sigma_{n}$, with
$c_{i}=0\textrm{ or }1$ and $\sigma_{i}$ being a $k$-simplex, and the vector
space $C_{k}$ contains all such objects.

It is obvious that the simplices are not independent to each other, for
example, if $\operatorname{\mathrm{Rips}}_{r}(P)$ contains a $2$-simplex
$\langle p_{0},p_{1},p_{2}\rangle$ as is displayed in Figure 3, it indicates
that all the pairwise distances $d(p_{0},p_{1})$, $d(p_{1},p_{2})$ and
$d(p_{2},p_{0})$ are less than or equal to the threshold $r$, hence
$\operatorname{\mathrm{Rips}}_{r}(P)$ definitely contains the $1$-simplices
$\langle p_{0},p_{1}\rangle$, $\langle p_{1},p_{2}\rangle$ and $\langle
p_{2},p_{0}\rangle$, and the $0$-simplices $\langle p_{0}\rangle$, $\langle
p_{1}\rangle$, $\langle p_{2}\rangle$ as well.

The relationship between $k$-simplices and $(k-1)$-simplices is described by
the so-called boundary operator $\partial_{k}$ which extracts the boundaries
of all $k$-simplices. In the case of $2$-simplex $\langle
p_{0},p_{1},p_{2}\rangle$ as displayed in Figure 3, the boundary contains
three $1$-simplices $\langle p_{0},p_{1}\rangle$, $\langle p_{1},p_{2}\rangle$
and $\langle p_{2},p_{0}\rangle$, we can describe it with the following
equation:

$\partial_{2}(\langle p_{0},p_{1},p_{2}\rangle)=\langle
p_{0},p_{1}\rangle+\langle p_{1},p_{2}\rangle+\langle p_{2},p_{0}\rangle.$ (1)

Figure 3: A $2$-simplex $\langle p_{0},p_{1},p_{2}\rangle$

One significant property of the boundary operators is that if you take the
boundary of a simplex twice, you always get $0$. In other words,

$\partial_{k-1}\circ\partial_{k}=0.$ (2)

To explain it, we still take the $2$-simplex $\langle
p_{0},p_{1},p_{2}\rangle$ as an example. We claim that
$\partial_{1}\circ\partial_{2}(\langle p_{0},p_{1},p_{2}\rangle)=0$. In fact,
for a $1$-simplex $\langle p_{i},p_{j}\rangle$,

$\partial_{1}\langle p_{i},p_{j}\rangle=\langle p_{i}\rangle+\langle
p_{j}\rangle,$ (3)

combining Eqn (1) and Eqn (3) together, we obtain that

$\partial_{1}\circ\partial_{2}(\langle p_{0},p_{1},p_{2}\rangle)=2\langle
p_{0}\rangle+2\langle p_{1}\rangle+2\langle p_{2}\rangle.$

Remind that in our computation of boundary operators, all the coefficients are
over the field $\mathbf{F}_{2}$, that is, $2=0$, therefore,
$\partial_{1}\circ\partial_{2}(\langle p_{0},p_{1},p_{2}\rangle)=0$.

According to the Eqn (2), all the images of the boundary operator
$\partial_{k-1}$ lie inside the kernel of the boundary operator
$\partial_{k}$, i.e.,
$\operatorname{\mathrm{im}}(\partial_{k-1})\subset\ker(\partial_{k})$. We then
define the Betti number $\beta_{k}$ as the difference of dimensions of
$\ker(\partial_{k})$ and $\operatorname{\mathrm{im}}(\partial_{k-1})$, that
is,
$\beta_{k}=\dim(\ker(\partial_{k}))-\dim(\operatorname{\mathrm{im}}(\partial_{k-1}))$.
Intuitively speaking, the Betti number $\beta_{k}$ gives the counting of
$k$-dimensional structures. More precisely, $\beta_{0}$ counts the connected
components, $\beta_{1}$ counts the circular structures, $\beta_{2}$ counts the
void structures, and $\beta_{k}$ counts the $n$-dimensional topological
structures.

######  Remark.

In the context of algebraic topology, the quotient space
$\ker(\partial_{k})/\operatorname{\mathrm{im}}(\partial_{k-1})$ is called the
$k$-th homology, and $\beta_{k}$ is just the dimension of the $k$-th homology.

In our context, we focus mainly on the first Betti numbers $\beta_{1}$, i.e.,
the counting of circular structures. To explain $\beta_{1}$ more explicitly,
we show how to compute $\beta_{1}$ of the boundary of the $2$-simplex in the
following example.

Figure 4: The boundary of the $2$-simplex $\langle p_{0},p_{1},p_{2}\rangle$.

######  Example 1.

Consider the boundary of the $2$-simplex $\langle p_{0},p_{1},p_{2}\rangle$ as
is displayed in Figure 4. According to the definition of simplex, it contains
three $1$-simplices: $\langle p_{0},p_{1}\rangle$, $\langle
p_{1},p_{2}\rangle$, and $\langle p_{2},p_{0}\rangle$, thus the space of
$1$-chains $C_{1}$ is a $3$-dimensional vector space over $\mathbf{F}_{2}$. We
then need to find the kernel of $\partial_{1}$, that is, to find out all the
solutions of the equation:

$\partial_{1}(c_{1}\langle p_{0},p_{1}\rangle+c_{2}\langle
p_{1},p_{2}\rangle+c_{3}\langle p_{2},p_{0}\rangle)=0.$

It is easy to see that the only solution is $c_{1}=c_{2}=c_{3}=1$, i.e.,
$\ker(\partial_{1})=\\{\langle p_{0},p_{1}\rangle+\langle
p_{1},p_{2}\rangle+\langle p_{2},p_{0}\rangle\\}$ and therefore
$\dim(\ker(\partial_{1}))=1$. We then need to find out
$\dim(\operatorname{\mathrm{im}}(\partial_{2}))$, however, there is no
$2$-simplex in the boundary of $\langle p_{0},p_{1},p_{2}\rangle$, hence the
space of $2$-chains $C_{2}$ is a zero dimensional vector space, and hence
$\dim(\operatorname{\mathrm{im}}(\partial_{2}))=0$. Therefore the Betti number
$\beta_{1}=1-0=1$ and we show that the boundary of a $2$-simplex $\langle
p_{0},p_{1},p_{2}\rangle$ has one circular structure.

As the threshold $r$ increases from $0$, we can obtain a series of different
Vietoris-Rips simplicial complex
$\\{\operatorname{\mathrm{Rips}}_{r}(P)\\}_{r\geq 0}$. Roughly speaking, the
procedure of TDA is to study the persistent behavior of Betti numbers with
respect to the threshold $r$.

The family of Betti numbers $\beta_{k}$ parameterized by the threshold $r$ is
referred to as the $k$-th persistent Betti numbers.

######  Remark.

In terms of TDA, the $k$-th persistent Betti numbers are the dimensions of the
$k$-th persistent homology parameterized by the threshold $r$.

More than just counting topological structures with persistent Betti numbers,
we can detect at which threshold values a topological structure is born and
dead.

Figure 5: The dim-$1$ persistence diagram of the pretzel-shaped point cloud

If for each topological structure we record the threshold values at which it
is born and dead, and plot them in a two-dimensional diagram, we can then have
a visualization of the topological structures of the point cloud $P$. More
precisely, for instance, if when the threshold $r$ is equal to $b$, a
topological structure is born, and when $r$ is equal to $d$, this topological
structure is dead, we then record the point $(b,d)$ on the diagram. We refer
to this diagram as a persistence diagram. In the case of pretzel-shaped point
cloud (Figure 1), the persistent diagram with respect to the first persistent
Betti numbers is displayed in Figure 5. We can see that there are three
remarkable points that are far away from the diagonal line, which represent
the three circular structures in Figure 1. The rest of the points in Figure 5
are closer to the diagonal line, which means their surviving time is much
shorter, and they are treated as noise.

######  Remark.

In the context of MCA, the more common feature values will be concentrated
around the original point, there would be less chance for them to form some
significant circular structures (In an extreme case, if all the samples share
the same value for each feature, all the feature values in the MCA plot will
collapse into the original point, hence no any circular structure will be
formed.). Even though any circular structure appears, it would be gone in a
very short time, which should be treated as noise. On the other hand, since
those low frequency feature values are further away from the original point,
if any circular structures could be constructed from them, such circular
structures will be survived for longer time and they should be more meaningful
in the framework of TDA.

We hope to quantify the difference between two point clouds by comparing their
corresponding persistence diagrams, however, it is usually difficult to
compare two persistence diagrams by visualization, therefore we also need some
suitable methods to measure the difference between different persistence
diagrams, by “suitable” we mean a “slight” fluctuation of sub-point clouds
should result in only a “slight” difference between the corresponding
persistence diagrams (Technically, this is referred to as “Stability Theorem”
in TDA ). Fortunately there are many such suitable options, one option is the
$p$-Wasserstein distance with $p>0$ being a parameter. Especially when
$p=\infty$, we call the $\infty$-Wasserstein distance the bottleneck distance.
We skip the exact definition of $p$-Wasserstein distance here since it is too
technical, the readers can check Definition B.15 in Appendix B for details.
Once we compute out the (topological) distances between each pair of
languages, many different data analysis methods such as clustering analysis,
kernel methods, multidimensional scaling (MDS) can then be applied. In sub-
Section VI.1 and sub-Section VI.2 we will apply the MDS to display a
visualization of two South American language families: Nuclear-Macro-Jê and
Quechuan.

##  V Data analysis procedure

In this paper, we analyze222A repository of the images used for this paper is
available at $\href https://github.com/ruidongsmile/Lingusitics_vs_topology$
only the languages of South America in Grambank(v1.0.3)
grambank_dataset_zenodo_v1 . The Grambank dataset is the largest comparative
grammatical database available. It covers $2467$ language varieties and a wide
range of grammatical phenomena in $195$ core features, thus providing
morphosyntactic profiles of languages across several linguistic subdomains (it
lacks information on phonology) grambank . The data, coded by a large team of
linguists, are freely accessible at https://grambank.clld.org/. Almost all the
features in Grambank are binary, the only six ternary features (GB024, GB025,
GB065, GB130, GB193, GB203) all concern word order: “what is the order of
elements X and Y?” with the alternatives “XY”, “YX”, or “both”. There are many
different ways to analyze discrete data, the MCA method is one suitable option
for us.

###  V.1 Data preprocessing

We first clean the dataset using R package “glottospace”glottospace . The
Grambank dataset contains $224$ languages of South America covering $41$
different language families and $30$ isolates, among which $217$ languages are
classified as “language” and $7$ languages are classified as “dialect”. To
keep the analysis unbiased, we keep only one dialect per language. As a
result, the languages katu1276 and tere1281 are removed. Since we restrict our
analysis to South American languages, some features (GB165,GB291,GB319,GB320)
contain only one unique value, which do not play a significant role and are
ignored in our analysis.

Another issue of the Grambank database is that the values of many features are
not known. For instance, in the case of language Suruí (suru1262), $180$ out
of $195$ features are missing. To make imputation not excessive, we remove
features and languages with more than $20\%$ missing values. We then apply
random forest algorithm to impute the rest missing data implemented by the R
package “mice” mice ; VANBUUREN2018 . Finally we split each ternary feature
into two binary features, as the form of “Is the order XY?” and “Is the order
YX?”. After the data preprocessing procedure, we obtain $183$ languages and
$60$ binary features.

###  V.2 MCA processing

We implement the MCA algorithm using the R package “ca” ca_pkg , here we
choose the method “adjusted” mjca  so as to get a better approximation.

Figure 6: The Scree plot of variance.

The scree plot of variance is displayed in Figure 6. The first two components
account for $66.9\%$ of the total variance, we plot the corresponding two-
dimensional point cloud in Figure 7. Here each point represents a feature
value.

Figure 7: The first two components of MCA projection. Each point represents a
feature value.

For instance, the point “GB054:1” represents the value $1$ of feature GB054.
Since in our context we have $60$ binary features, hence the point cloud
obtained via MCA contains $120$ points, and each sub-point cloud contains $60$
points corresponding to the values of the $60$ features.

###  V.3 The shapes of languages

We can now visualize each language by filtering out the feature values. We
pick up three languages from each family of Quechuan and Nuclear-Macro-Jê
(NMJ) as examples, their corresponding sub-point clouds are displayed in
Figure 8. We can see that the languages in the same family have similar
shapes. Specifically, the sub-point clouds of languages in NMJ family have a
circular structure, while there is not such a circular structure in the
languages of the Quechuan family.

Figure 8: The sub-point clouds of some languages from families Quechuan and
Nuclear-Macro-Jê.

The persistence diagrams of corresponding sub-point clouds are displayed in
Figure 9. We observe that in the persistence diagrams of languages from
Quechuan family (the first row in Figure 9), there is not any point
significantly away from the diagonal line, while in the case of NMJ family
(the second row in Figure 9), there is a remarkable point away from the
diagonal line, which means that there is a significant circular structure in
the corresponding sub-point cloud. This phenomenon implies that the shape of a
sub-point cloud encodes some characteristic of the corresponding language,
although further research is needed to understand the linguistic meaning.

Figure 9: The persistence diagrams with respect to circular structures of some
languages from families Quechuan and Nuclear-Macro-Jê.

##  VI Applications of TDA

In this section we show how to apply the TDA method to analyze the Nuclear-
Macro-Jê (NMJ) family and the Quechuan family. According to the scree plot
(Figure 6) obtained from the MCA method, we pick the first four components,
which account for $76.8\%$ of the total variance, and we analyze the circular
structures (the first homology in the term of topology) of each sub-point
cloud. In this case, all the sub-point clouds will be four-dimensional, and it
is impossible to visualize them directly, however we can still detect there
circular structures via TDA.

###  VI.1 Nuclear-Macro-Jê

Despite some relatively recent advances, Nuclear-Macro-Jê (NMJ) remains one of
the lesser understood families of South America. The family is generally
represented as having a rake-like structure, consisting of a number of first-
order members, some of which have further branches. The first-order branch for
which most reconstructive work has been done is the Jê branch (or Jê-proper)
(e.g. nikulin_proto_mj, ), for which there is general consensus that they form
a genealogical sub-unit, in turn divided into a northern, central, and
southern branch. For the other first-order branches, the support ranges widely
michael_classifications , but the languages in our sample that are not from
Jê-proper are generally accepted as part of NMJ.

We will apply the TDA method to analyze the difference between Jê-proper
languages and the other NMJ languages. After the data-cleaning processing,
there are $11$ NMJ languages remained (Grambank contains $13$ languages of NMJ
in total), $5$ of which belong to the Jê-proper branch. The sample languages
are divided as follows over the first-order branches (classification follows
glottolog ):

  * •

Jê-proper: Apinayé (apin1244), Kanela-Krahô (cane1242), Panará (pana1307),
Xavánte (xava1240), Xokleng (xokl1240).

  * •

Jabutí: Djeoromitxi (djeo1235)

  * •

Karajá: Karajá

  * •

Maxakali-Borum: Borum/Krenák (kren1239), Maxakalí (maxa1247)

  * •

Ofayé: Ofayé (ofay1240)

  * •

Rikbaktsa: Rikbaktsa (rikb1245)

Figure 10: Geographical map of Nuclear-Macro-Jê (NMJ)

The persistence diagrams of Jê-proper and the other NMJ languages are
displayed in Figure 11 and Figure 12, respectively. We observe that there is a
special circular structure in the Jê-proper group when birth is around $0.10$
and death is around $0.15$, the only exception is the language pana1307 (
Panará), however there is no such a significant circular structure in the
group of non-Jê-proper languages. As mentioned, Panará (pana1307) forms an
exception within the Jê-proper group. The language has been described as
having a deviant morphosyntactic profile compared to its sister languages It
is at this point unclear what the sociohistorical explanation of this deviance
is, but it is possibly related to a turbulent past in the last century or so,
involving long-distance migrations and a severe population bottleneck as the
result of exposure to diseases and ecological disruption of their habitat
dourado_phd_panara , (nikulin_proto_mj, , page 11-12).

Figure 11: The persistence diagrams with respect to circular structures of Jê-
proper languages Figure 12: The persistence diagrams with respect to circular
structures of non-Jê-proper languages

We can detect which points contribute to this unique circular structure in
each sub-point cloud using the R package “TDA” r_TDA , the points are listed
as below:

  * •

apin1244: GB044:0, GB053:0, GB079:1, GB084:0, GB107:0, GB028:1, GB035:1,
GB093:0, GB039:0, GB318:1, GB089:0, ”GB136:1”.

  * •

cane1242: GB044:0, GB068:1, GB053:0, GB107:0, GB028:1, GB084:0, GB093:0,
GB039:0, GB318:1, ”GB136:1.

  * •

xava1240: GB044:0, GB042:0, GB053:0, GB039:0, GB079:1, GB084:0, GB107:0,
GB093:0, GB318:1, GB193_YX:1, GB089:0, GB136:1.

  * •

xokl1240: GB044:0, GB068:1, GB021:0, GB107:0, GB084:0, GB318:1, GB136:1.

If we take the intersection of all the points, we notice that the five points
GB044:0, GB084:0 GB107:0, GB318:1 and GB136:1 appear in all four languages,
while in all the other languages there is not such a combination. In Panará,
the values of these five features are: GB044:1, GB084:0, GB107:0, GB318:0 and
GB136:1. The mini profile of these five feature values can be connected to the
isolating tendencies of the Jê-proper languages, and their dependence on
syntactic rather than morphological strategies (rodrigues1999, , 180), as they
concern the overt marking of inflectional categories on nouns (plural) and
verbs (future tense, negation), and the presence of fixed clausal word order.
Panará, having morphologized its number marking on nouns, is slightly less
isolating than its sister languages, which use a prosodically free word to
mark plural.

We can also measure the bottleneck distance and $2$-Wasserstein distance
between different languages in NMJ. This can be implemented by using the R
package “TDApplied” r_TDApplied . We then can apply the MDS method to
visualize the corresponding distance matrix, which is displayed in Figure 13.
We notice that there is a clear gap between Jê-proper group and non-Jê-proper
group.

Figure 13: Left: MDS plot w.r.t. bottleneck distance of Nuclear-Macro-Jê.
Right: MDS plot w.r.t. $2$-Wasserstein distance of Nuclear-Macro-Jê.

Moreover, a hypothesis test for non-parametric statistical inference of
persistence diagrams can be conducted. We will use a permutation test method
introduced in MR3975554  to test if the Jê-proper group and the non-Jê-proper
groups are significantly distinguishable.

Our null hypothesis is that the sets of persistence diagrams associated with
each group are indistinguishable. The first step is to define a function that
can quantify the difference between two groups (In terms of mathematical
optimization or decision theory, this is called a loss function) based on the
chosen measurement between two persistence diagrams. Here we choose the
bottleneck distance measurement and $2$-Wasserstein distance measurement
separately. Then we calculate the value of the loss function, which is
referred to as the observed value. The next step is to randomly shuffle all
the persistence diagrams multiple times and calculate the corresponding values
of loss function. Finally we evaluate the ratio of values that are less than
or equal to the observed value, which accounts for the $p$-value of our
permutation test. There are many R packages like “TDApplied” and “TDAstats”
which can implement the permutation test, in this experiment we choose
“TDApplied” and shuffle $100$ times. in the case of $2$-Wasserstein distance,
we get a $p$-value of $0.0297(=3/101)$, and in the case of Bottleneck
distance, we get a $p$-value of $0.0099(=1/101)$. We observe that in both
cases our $p$-value is less than the threshold $0.05$, hence our null
hypothesis is rejected.

######  Remark.

In our computation, the $p$-value $0.0297$ associated with $2$-Wasserstein
distance is a round value of $3/101$, which means that among all the $100$
times of shuffles, only $2$ times we get a value of loss function less than or
equal to the observed value. And the $p$-value $0.0099$ associated with
bottleneck distance is a round value of $1/101$, which means that among all
the $100$ times of shuffles, there is not any value of loss function that is
less than or equal to the observed value.

In fact, since in the family of NMJ there are only $11$ samples which are
divided into two groups of size $5$ and $6$ separately, we have only $462$
different permutations, the calculation is relatively simple and we can get
the exact $p$-value.

By computing the loss function over all the $462$ permutations, we realize
that in the measurement of $2$-Wasserstein distance, $9$ out of $462$ values
of loss function are less than or equal to the observed value, hence the exact
$p$-value is equal to $0.0195$, and in the measurement of bottleneck distance,
$8$ out of $462$ values of loss function are less than or equal to the
observed value, the exact $p$-value is equal to $0.0173$.

###  VI.2 Quechuan

The Quechuan language family stretches along the Andean mountain range from
southern Colombia to northern Chile and Argentina, and contains about 40
languages. Most of its current spread was achieved only relatively recently,
with the expansion of the Inca empire (13th-16th century CE), of which Quechua
was the language of administration. The core of Quechuan, however, in central
Peru, is of greater time depth, estimated to about 2000 years (adelaar2004, ,
page 168). The family is divided into two main branches: Quechua I (in central
Peru) and Quechua II (corresponding to the later expansions north- and
southward). The northern Quechuan languages in north Peru, Ecuador and
Colombia have replaced previously existing languages, and have in the process
undergone substrate influence Muysken2021 , making them somewhat deviant
grammatically from the rest of the family.

Grambank contains $9$ Quechuan languages, of which $8$ are remained after the
data preprocessing procedure. Among these $8$ languages, the languages
Calderón Quechua (cald1236), Chimborazo Quechua (chim1302) and Imbabura
Quechua (imba1240) are in the north, and the languages San Martín Quechua
(sanm1289), Ayacucho Quechua (ayac1239), Cuzco Quechua (cusc1236), Huallaga
Quechua (hual1241), North Junín Quechua (nort2980) are in the south. North
Junín Quechua and Huallaga Quechua are part of the Quechua I branch, whereas
the others are part of the Quechua II branch. The northern (Ecuadorian)
Quechuan languages form a sub-branch within Quechua II.

Figure 14: Quechuan Languages

The persistence diagrams of Quechuan languages in north and south are shown in
Figure 15 and Figure 16 respectively.

Figure 15: The persistence diagrams with respect to circular structures of
Quechuan languages in the north Figure 16: The persistence diagrams with
respect to circular structures of Quechuan languages in the south

In the case of the Quechuan family it is relatively difficult to detect the
difference of the persistence diagrams just by visualization, however, we can
measure the $p$-Wasserstein distance of each pair of languages. We still
choose $2$-Wasserstein measurement, the Figure 17 displays the $3$-dimensional
MDS plot of Quechuan languages. As it is shown, the north Quechuan languages
(red) and south Quechuan languages (blue) are separated into two different
groups. We then do a permutation test to make a decision as well. Our null
hypothesis is that the north Quechuan languages and south Quechuan languages
are indistinguishable. In this example since the size of the two groups are
$3$ and $5$ respectively, there are only $56$ different permutations in total,
hence it is possible to calculate the loss function of each permutation.
Finally we obtain that the exact $p$-value is equal to $0.0179$. If we set our
threshold to be $0.05$, our null hypothesis is then denied. In fact, the
$p$-value $0.0179$ is equal to $\frac{1}{56}$, which means that every non-
trivial permutation gives a value of loss function larger than the observed
value.

Figure 17: The $3$-dim MDS visualization of $2$-Wasserstein distance of
Quechuan

######  Remark.

It is worth to emphasize it again that the persistence diagrams in Figure 9
are different from those in Figure 11, Figure 12, Figure 15 and Figure 16, The
reason is that in Figure 9 we take the first two components of MCA projection,
therefore all the point clouds are inside a $2$-dimensional Euclidean space,
while in the rest cases we take the first four components of MCA projection,
and the point clouds are inside a $4$-dimensional Euclidean space.

##  VII Discussion

We have presented the visualization of languages via MCA method. More than
that, we showed how to apply the TDA framework to analyze the shapes of
languages in the case of Nuclear-Macro-Jê (NMJ) family and Quechuan family
respectively. In addition, we measured the bottleneck distance and
$2$-Wasserstein distance of persistence diagrams. We see that there is a
significant difference between Jê-proper and non-Jê-proper via the $2$-dim MDS
plot. In the case of the Quechuan family, we can see the difference between
north Quechuan languages and south Quechuan languages via the $3$-dim MDS plot
with respect to $2$-Wasserstein measurement. Furthermore, if we apply the
permutation test, by computing the $p$-value we see that in the case of the
NMJ family, Jê-proper and non-Jê-proper are significantly different from each
other, and in the case of the Quechuan family, the north Quechuan languages
and south Quechuan languages are significantly different as well.

The drawback of persistence diagram is that it contains only topological
information. In an extreme case, suppose if the corresponding point clouds of
two languages are congruent, i.e., one point cloud can be transformed into the
other point cloud by a rigid transformation (i.e., a composition of rotations,
translations and reflections), the two persistence diagrams would be exactly
the same. Nevertheless, in the framework of topological data analysis there
are many more powerful tools like “persistent homology transform” MR3311455
that can deal with this issue. Another drawback is that the persistence
diagram can detect only the topological structures like clustering structures
and circular structures of a point cloud. In contrast to that, some recent
works book_jost ; pers_lap  indicate that the spectrum (or in other words, the
eigenvalues) of the persistent Laplacian operators contain much more deeper
geometrical information of a Vietoris-Rips filtration. It would be an
interesting topic to analyze what the spectrum of the persistent Laplacian
operators can tell about a language.

## Acknowledgements

This work is part of the ERC Consolidator project South American Population
History Revisited, funded by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 818854 - SAPPHIRE). This support is gratefully acknowledged. I
would like to express my sincere gratitude to Rik van Gijn for the valuable
discussion about South American languages and proofreading this work. The main
material of this work is based on a presentation I gave in the 2023 Grambank
Workshop in Leipzig, I feel thankful to Russell Gray and Hedvig Skirgård for
their insightful suggestions. In addition, I also would like to express my
thanks to Leonardo Arias Alvis and Sietze Norder for their inspiring feedback.

##  Appendix A Multiple correspondence analysis

In this work we use the multiple correspondence analysis (MCA) based on the
Burt Matrix with adjustment. We follow the notation of mjca . For a more
complete introduction to MCA, we refer the readers to mjca ; MCA .

Consider a table of $Q$ categorical features and $I$ samples, we denote by
$J_{q}$ the amount of categorical values the $q$-th feature contains, then in
total we have $J=\sum_{q=1}^{Q}J_{q}$ different categorical values. We then
create an $I\times J$ matrix $Z$, each row corresponds to a sample, each
column corresponds to a categorical value of a feature. The entries of matrix
$Z$ are either $0$ or $1$, depends on the feature values of samples, if a
sample has a feature value, the corresponding entry is $1$, otherwise, it is
$0$. The matrix $Z$ is called the indicator matrix. We then define the Burt
matrix $B$ as $B=\\{b_{ij}\\}=Z^{T}Z$.

The MCA algorithm is implemented in the following steps:

  1. 1.

Compute the correspondence matrix $P$:

$P=\\{p_{ij}\\}={b_{ij}\Big{/}\sum_{i,j}b_{ij}}.$

  2. 2.

Compute the row totals $r_{i}$:

$r_{i}=\sum_{j}p_{ij}.$

  3. 3.

Compute the standardized residuals matrix $S=\\{s_{ij}\\}$ with
$s_{ij}=(p_{ij}-r_{i}r_{j})\big{/}\sqrt{r_{i}r_{j}}.$

  4. 4.

Perform the singular value decomposition (SVD) on $S$: $S=V\Lambda V^{T},$
here $\Lambda=$diag$\left(\lambda_{1},\lambda_{2},\cdots,\lambda_{J}\right)$
with $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{J}$.

  5. 5.

Compute the adjusted variance $\widetilde{\lambda}_{j}^{2}$:

$\widetilde{\lambda}_{j}^{2}=\begin{cases}\left(\frac{Q}{Q-1}\right)^{2}\left(\lambda_{j}-\frac{1}{Q}\right)^{2}&\text{if
$\lambda_{j}>\frac{1}{Q}$},\\\ 0&\text{if
$\lambda_{j}\leq\frac{1}{Q}$}.\end{cases}$

  6. 6.

Compute the principal coordinates of all samples
$F=D_{r}^{-\frac{1}{2}}V\widetilde{\Lambda}$, here $D_{r}^{-\frac{1}{2}}$ and
$\widetilde{\Lambda}$ are diagonal matrices defined by
$D_{r}^{-\frac{1}{2}}=$diag$(r_{1}^{-1/2},r_{2}^{-1/2},\cdots,r_{J}^{-1/2})$
and
$\widetilde{\Lambda}=$diag$\left(\widetilde{\lambda}_{1},\widetilde{\lambda}_{2},\cdots,\widetilde{\lambda}_{J}\right)$.

  7. 7.

Compute the adjusted total variance:

$\tau=\frac{Q}{Q-1}\left(\sum_{j=1}^{J}\lambda_{j}^{2}-\frac{J-Q}{Q^{2}}\right).$

  8. 8.

Evaluate the percentage of variance for the $j$-th principal component:
$\widetilde{\lambda}^{2}_{j}\Big{/}\tau$.

##  Appendix B Background about topological data analysis

Suppose $P$ is a point cloud of finitely many points inside a Euclidean space
$\mathbb{R}^{d}$. In this appendix we will give the main definitions and
results necessary in this paper. We follow mainly the notation of the book
MR4381505 . We refer to the readers MR2572029 ; MR4381505 ; MR2549932  for
more details about the framework of TDA.

###  B.1 Simplices and complexes

######  Definition B.1.

We say a set of points $\\{p_{0},p_{1},...,p_{k}\\}\subset\mathbb{R}^{d}$ is
affinely independent if for any set of real numbers
$\\{t_{0},t_{1},...,t_{k}\\}\subset\mathbb{R}$,

$\sum_{i=0}^{k}t_{i}=0\textrm{ and }\sum_{i=0}^{k}t_{i}p_{i}=0$

implies $t_{0}=t_{1}=...=t_{k}=0$.

The set $\\{p_{0},p_{1},...,p_{k}\\}$ is affinely independent if and only if
the $k$ vectors $\\{p_{1}-p_{0},p_{2}-p_{1},...,p_{k}-p_{k-1}\\}$ are linearly
independent.

######  Definition B.2 (Simplex).

If the set of points $P=\\{p_{0},p_{1},...,p_{k}\\}$ is affinely independent,
the convex hull of $P$ is called a $k$-simplex.

######  Remark.

In more detail, the convex hull of $P$ is the smallest convex set that
contains $P$. Follow this definition, a $0$-simplex is a point, a $1$-simplex
is a line segment, a $2$-simplex is a triangle, and a $3$-simplex is a
tetrahedron.

A $k$-simplex $\sigma$ is said to have dimension $k$. For a subset
$P^{\prime}\subset P$ that contains $k^{\prime}$ points, the convex hull of
$P^{\prime}$ is a $k^{\prime}$-simplex, we call it a $k^{\prime}$-face of
$\sigma$.

######  Definition B.3 (Simplicial complex).

We call $K$ a simplicial complex if it is a finite collection of simplices and
satisfies the following two restrictions:

  * •

Each face of each simplex is contained in $K$.

  * •

Take any two simplices $\sigma_{1}$ and $\sigma_{2}$ in $K$, the intersection
$\sigma_{1}\cap\sigma_{2}$ is either a face of both $\sigma_{1}$ and
$\sigma_{2}$ or empty.

We define the dimension of a simplicial complex $K$ to be the maximum
dimension of all simplices in $K$, and denoted it as $\dim(K)$. We say $K$ is
a simplicial $k$-complex if $\dim(K)=k$.

###  B.2 Chains and boundary operators

Let $K$ be a simplicial $k$-complex with $n_{p}$ being the number of
$p$-simplices, according to Definition B.3 it is easy to see that $k\geq p\geq
0$.

######  Definition B.4 (Chains).

A $p$-chain $\gamma$ in a simplicial complex $K$ over the field
$\mathbf{F}_{2}$ is defined to be a formal sum of $p$-simplices, that is,
$\gamma=\sum_{i=1}^{n_{p}}c_{i}\sigma_{i}$ with $c_{i}\in\mathbf{F}_{2}$ and
each $\sigma_{i}$ being a $p$-simplex.

The set of all $p$-chains form a vector space over the field $\mathbf{F}_{2}$,
we can sum two $p$-chains $\gamma=\sum c_{i}\sigma_{i}$ and
$\gamma^{\prime}=\sum c_{i}^{\prime}\sigma_{i}$ up to obtain another $p$-chain
$\gamma+\gamma^{\prime}=\sum(c_{i}+c_{i}^{\prime})\sigma_{i}$, and we can also
multiply a $p$-chain $\gamma$ by a number $c^{\prime}\in\mathbf{F}_{2}$ to get
a new chain $c^{\prime}\gamma=\sum_{i=1}^{n_{p}}c^{\prime}c_{i}\sigma_{i}$. We
denote the set of all $p$-chains as $C_{p}(K)$. Since the simplicial complex
$K$ is clear from the context we always drop $K$ and use the notation $C_{p}$
to denote the set of all $p$-th chains.

There is a group homomorphism from $C_{p}$ to $C_{p-1}$ which is referred to
as a boundary operator. the definition is given as follows.

######  Definition B.5 (Boundary operator).

We take $\sigma=\langle v_{0},...,v_{p}\rangle$ to be a $p$-simplex and let

$\partial_{p}(\sigma):=\sum_{i=0}^{p}\langle
v_{0},...,\hat{v}_{i},...,v_{p}\rangle,$

where $\hat{v}_{i}$ indicates that $v_{i}$ is removed. Since $C_{p}$ is a
vector space, we can then extend $\partial_{p}$ to $C_{p}$ and obtain a group
homomorphism $\partial_{p}:C_{p}\to C_{p-1}$ as

$\partial_{p}\gamma=\sum_{i=1}^{m_{p}}c_{i}\partial_{p}(\sigma_{i})$

for a $p$-chain $\gamma=\sum c_{i}\sigma_{i}\in C_{p}$. The map $\partial_{p}$
is called a boundary operator.

The most important property of the boundary operator $\partial_{p}$ is that
the composition of two boundary operators is zero.

######  Theorem B.1.

For any $p$-chain $\gamma$, $\partial_{p-1}\cdot\partial_{p}(\gamma)=0$.

We can then obtain the following sequence of group homomorphisms:

$0=C_{k+1}\xrightarrow{\partial_{k+1}}C_{k}\xrightarrow{\partial_{k}}\cdots
C_{1}\xrightarrow{\partial_{1}}C_{0}\xrightarrow{\partial_{0}}C_{-1}=0.$

We refer to such a sequence as a chain complex.

######  Definition B.6 (Boundary group).

We call the image of $\partial_{p+1}$ the $p$-th boundary group, and denote it
as $\mathsf{B}_{p}$. That is, $\mathsf{B}_{p}=\big{\\{}\gamma\in
C_{p}\,\big{|}\,\partial_{p+1}\tilde{\gamma}=\gamma\textrm{ for some
}\tilde{\gamma}\in C_{p+1}\big{\\}}$

######  Definition B.7 (Cycle group).

We call the kernel of $\partial_{p}$ the $p$-th cycle group, and denote it as
$\mathsf{Z}_{p}$. In more detail, $\mathsf{Z}_{p}=\big{\\{}\gamma\in
C_{p}\,\big{|}\,\partial_{p}\gamma=0\big{\\}}.$

Since the composition of two boundary operators $\partial_{p}$ and
$\partial_{p-1}$ is zero by Theorem B.1, we can realize that
$\mathsf{B}_{p}\subset\mathsf{Z}_{p}$, and the quotient group plays an
important role in algebraic topology.

######  Definition B.8 (Homology group).

The $p$-th homology group is defined as the quotient group
$\mathsf{H}_{p}:=\mathsf{Z}_{p}/\mathsf{B}_{p}$.

In our context, the homology group $\mathsf{H}_{p}$ is a vector space over the
field $\mathbf{F}_{2}$, the dimension of $\mathsf{H}_{p}$ is called the $p$-th
Betti number, we denote it as $\beta_{p}$.

###  B.3 Basics about persistent homology

######  Definition B.9 (Vietoris-Rips complex).

We define the Vietoris-Rips complex of a point cloud $P$ at threshold $r$ to
be the simplicial complex

$\operatorname{\mathrm{Rips}}_{r}(P)=\big{\\{}\sigma\subset
P\big{|}\textrm{diam}(\sigma)\leq 2r\big{\\}},$

here $\textrm{diam}(\sigma)$ denotes the diameter of a subset $\sigma$, i.e.,
$\textrm{diam}(\sigma)=\max\limits_{p_{i},p_{j}\in\sigma}||p_{i}-p_{j}||$.

######  Definition B.10.

We define the Vietoris-Rips filtration of $P$ to be the nested collection of
Vietoris-Rips complexes of $P$:

$\operatorname{\mathbb{R}ips}(P)=\\{\operatorname{\mathrm{Rips}}_{r}(P)\\}_{r\in[0,\infty)}.$

The Definition B.10 indicates that if $r_{1}<r_{2}$, the Vietoris-Rips complex
$\operatorname{\mathrm{Rips}}_{r_{1}}(P)$ is a sub-complex of
$\operatorname{\mathrm{Rips}}_{r_{2}}(P)$, which induces a group homomorphism
$h_{p}^{r_{1},r_{2}}$ in the level of homology groups. In more detail,

$\begin{split}\operatorname{\mathrm{Rips}}_{r_{1}}(P)&\subset\operatorname{\mathrm{Rips}}_{r_{2}}(P)\\\
&\Big{\Downarrow}\\\
\mathsf{H}_{p}(\operatorname{\mathrm{Rips}}_{r_{1}}(P))&\xrightarrow{h_{p}^{r_{1},r_{2}}}\mathsf{H}_{p}(\operatorname{\mathrm{Rips}}_{r_{2}}(P)).\end{split}$
(4)

######  Definition B.11 (Persistent homology).

The $p$-th persistent homology groups are defined as the images of the group
homomorphisms $h_{p}^{r_{1},r_{2}}$, that is to say,
$\mathsf{H}_{p}^{r_{1},r_{2}}:=\operatorname{\mathrm{im}}\,h_{p}^{r_{1},r_{2}}$.

######  Definition B.12 (Persistent Betti number).

We define the dimensions $\dim\mathsf{H}_{p}^{r_{1},r_{2}}$ of the vector
spaces $\mathsf{H}_{p}^{r_{1},r_{2}}$ as the $p$-th persistent Betti numbers,
and denote it as $\beta_{p}^{r_{1},r_{2}}$.

######  Definition B.13 (Birth and death).

Let $\xi$ be a nontrivial $p-$th homology class in
$\mathsf{H}_{p}(\operatorname{\mathrm{Rips}}_{a}(P))$. For a number $r_{2}\leq
a$, we say $\xi$ is born at $\operatorname{\mathrm{Rips}}_{r_{2}}(P)$ if
$\xi\in\mathsf{H}_{p}^{r_{2},a}$ while $\xi\notin\mathsf{H}_{p}^{r_{1},a}$ for
any $r_{1}<r_{2}$. Similarly, for a number $s_{1}\geq a$, we say $\xi$ dies at
$\operatorname{\mathrm{Rips}}_{s_{1}}(P)$ if $h_{p}^{a,s_{1}}(\xi)$ does not
vanish while $h_{p}^{a,s_{2}}(\xi)=0$ for any $s_{2}>s_{1}$.

We can visualize the birth and death of homology classes via the so-called
persistence diagram. For a nontrivial $p-$th homology class
$\xi\in\mathsf{H}_{p}(\operatorname{\mathrm{Rips}}_{r}(P))$, if it is born at
$b$ and dies at $d$, we denote it as a point $(b,d)$ in the plane
$\mathbb{R}^{2}$, as is shown in Figure 5.

Suppose we have two point clouds $P_{1}$ and $P_{2}$, we hope to work out some
method to measure the difference between $P_{1}$ and $P_{2}$ via the
corresponding persistence diagrams. To start with, we review the definition of
matching between two persistence diagrams.

######  Definition B.14 (Matching).

Let $P_{1}$ and $P_{2}$ be two point clouds and $\operatorname{dgm}_{p,1}$,
$\operatorname{dgm}_{p,2}$ be the corresponding persistence diagrams with
respect to the $p$-th persistent homology. We denote by $A_{p,1},A_{p,2}$ the
non-diagonal points in the persistence diagrams $\operatorname{dgm}_{p,1}$ and
$\operatorname{dgm}_{p,2}$ respectively. For a point $a_{1}\in A_{p,1}$, let
$\overline{a}_{1}$ denote the nearest point (in the sense of Euclidean
distance) of $a_{1}$ on the diagonal. Define $\overline{a}_{2}$ for every
point $a_{2}\in A_{p,2}$ similarly. Let
$\overline{A}_{p,1}=\\{\overline{a}_{1}\\}_{a_{1}\in A_{p,1}}$ and
$\overline{A}_{p,2}=\\{\overline{a}_{2}\\}_{a_{2}\in A_{p,2}}$. Let
$\tilde{A}_{p,1}=A_{p,1}\bigcup\overline{A}_{p,2}$ and
$\tilde{A}_{p,2}=A_{p,2}\bigcup\overline{A}_{p,1}$. We call a bijective map
$\varphi$ from $\tilde{A}_{p,1}$ to $\tilde{A}_{p,2}$ a matching.

######  Definition B.15 (Wasserstein distance and bottleneck distance).

Let $\Phi$ be the set of all matchings as defined in Definition B.14. For any
$q\in[1,+\infty]$, the $q-$Wasserstein distance is defined as

$\displaystyle d_{q}(\operatorname{dgm}_{p,1},\operatorname{dgm}_{p,2})$
$\displaystyle=$
$\displaystyle\inf_{\varphi\in\Phi}\left(\sum_{x\in\operatorname{dgm}_{p,1}}(||x-\varphi(x)||_{q})^{q}\right)^{1/q}.$

Especially when $q=+\infty$ the distance $d_{\infty}(\cdot,\cdot)$ is referred
to as the bottleneck distance.

The significance of Wasserstein distance and bottleneck distance is that a
“slight modification” of the point cloud results in only a “slight
modification” of the corresponding persistence diagram. We denote by $|P|$ the
cardinality of the point cloud $P$. In the following two theorems we assume
$P_{1}$ and $P_{2}$ are two point clouds with $|P_{1}|=|P_{2}|=M<\infty$.

######  Theorem B.2 (Bottleneck stability theorem).

Let $P_{1}$ and $P_{2}$ be two point clouds, if there is a number $\epsilon>0$
and a bijection $\varphi:P_{1}\to P_{2}$ such that
$||x-\varphi(x)||\leq\epsilon$ for all $x\in P_{1}$. Then

$d_{\infty}(\operatorname{dgm}_{p,1},\operatorname{dgm}_{p,2})\leq\epsilon.$

######  Definition B.16.

Let $\Phi$ be the set of all bijections between point clouds $P_{1}$ and
$P_{2}$. We define the point cloud Wasserstein distance between $P_{1}$ and
$P_{2}$ as

$W_{p}(P_{1},P_{2})=\inf_{\varphi\in\Phi}\left(\sum_{x\in
P_{1}}||x-\varphi(x)||_{p}^{p}\right)^{1/p}.$

In the case of Wasserstein distance we also have the following Wasserstein
Stability theorem skraba2023wasserstein :

######  Theorem B.3 (Wasserstein stability theorem skraba2023wasserstein ).

$d_{q}(\operatorname{dgm}_{p,1},\operatorname{dgm}_{p,2})\leq
2^{M/(q+1)}W_{q}(P_{1},P_{2}).$

###  B.4 Permutation test

We briefly repeat the procedure of permutation test based on TDA. The reader
can check MR3975554  for a more detailed discussion.

Let $\mathbf{P}_{1}=\\{P_{1,1},...,P_{1,n_{1}}\\}$ and
$\mathbf{P}_{2}=\\{P_{2,1},...,P_{2,n_{2}}\\}$ be two groups of point clouds.
Our purpose is to detect if the shapes of the point clouds $\mathbf{P}_{1}$
and the point clouds $\mathbf{P}_{2}$ are significantly distinguishable. We
first compute the persistence diagrams with respect to the $p$-th persistent
homology of each point cloud, and denote the corresponding sets of persistence
diagrams by $\mathbf{D}_{1}=\\{D_{1,1},...,D_{1,n_{1}}\\}$ and
$\mathbf{D}_{2}=\\{D_{2,1},...,D_{2,n_{2}}\\}$.

We then need to pick up a Wasserstein distance $d_{q}(\cdot,\cdot)$ between
two persistence diagrams. Let $\operatorname{\mathrm{Avg}}(\mathbf{D}_{1})$ be
the average distance among all the persistence diagrams in the group
$\mathbf{D}_{1}$, i.e.,

$\operatorname{\mathrm{Avg}}(\mathbf{D}_{1})=\frac{1}{n_{1}(n_{1}-1)}\sum_{i,j=1}^{n_{1}}d_{q}(D_{1,i},D_{1,j}),$

and $\operatorname{\mathrm{Avg}}(\mathbf{D}_{2})$ is defined similarly for the
set of persistence diagrams $\mathbf{D}_{2}$.

We define a loss function $F$ of $\mathbf{D}_{1}$ and $\mathbf{D}_{2}$ as

$F(\mathbf{D}_{1},\mathbf{D}_{2})=\frac{1}{2}\left(\operatorname{\mathrm{Avg}}(\mathbf{D}_{1})+\operatorname{\mathrm{Avg}}(\mathbf{D}_{2})\right).$

We then randomly shuffle all the persistence diagrams in $\mathbf{D}_{1}$ and
$\mathbf{D}_{2}$ together, and form two new groups
$\widetilde{\mathbf{D}}_{1}$, $\widetilde{\mathbf{D}}_{2}$ of size $n_{1}$,
$n_{2}$ respectively, then compute the value
$F(\widetilde{\mathbf{D}}_{1},\widetilde{\mathbf{D}}_{2})$ and check whether
$F(\widetilde{\mathbf{D}}_{1},\widetilde{\mathbf{D}}_{2})\leq
F(\mathbf{D}_{1},\mathbf{D}_{2})$ or not.

We repeat this shuffling procedure for $N$ times, if among all the $N$ times
of experiments, $Z$ out of $N$ times it happens that
$F(\widetilde{\mathbf{D}}_{1},\widetilde{\mathbf{D}}_{2})\leq
F(\mathbf{D}_{1},\mathbf{D}_{2})$, the $p$-value is then defined to be
$\frac{Z+1}{N+1}$.

Intuitively speaking, imagine if all the point clouds in $\mathbf{D}_{1}$ are
identical and all the point clouds in $\mathbf{D}_{2}$ are identical as well
while $\mathbf{D}_{1}$ are different from $\mathbf{D}_{2}$, then
$F(\mathbf{D}_{1},\mathbf{D}_{2})=0$ while
$F(\widetilde{\mathbf{D}}_{1},\widetilde{\mathbf{D}}_{2})\neq 0$ for every
non-trivial permutation. Thus in this case the counting value $Z$ is equal to
$0$, and the corresponding $p$-value is then equal to $\frac{1}{N+1}$.

## References

  * [1] Hervé Abdi and Dominique Valentin.  Multiple correspondence analysis.  Encyclopedia of Measurement and Statistics, 01 2007.
  * [2] Willem Adelaar and Pieter Muysken.  The languages of the Andes.  Cambridge Language Surveys. Cambridge University Press, Cambridge, 2004\.
  * [3] Shael Brown and Dr. Reza Farivar.  TDApplied: Machine Learning and Inference for Topological Data Analysis, 2023.  R package version 2.0.4.
  * [4] Tamal Krishna Dey and Yusu Wang.  Computational topology for data analysis.  Cambridge University Press, Cambridge, 2022.
  * [5] Luciana Dourado.  Aspectos Morfossintáticos da Língua Panará (Jê).  PhD thesis, Universidade Estadual de Campinas, 2001.
  * [6] Herbert Edelsbrunner and John L. Harer.  Computational topology.  American Mathematical Society, Providence, RI, 2010.  An introduction.
  * [7] Brittany T. Fasy, Jisu Kim, Fabrizio Lecci, Clement Maria, David L. Millman, and Vincent Rouvreau.  TDA: Statistical Tools for Topological Data Analysis, 2023.  R package version 1.9.
  * [8] Michael Greenacre and Oleg Nenadic.  Computation of multiple correspondence analysis, with code in r.  Department of Economics and Business, Universitat Pompeu Fabra, Economics Working Papers, 09 2005.
  * [9] Harald Hammarström, Robert Forkel, Martin Haspelmath, and Sebastian Bank.  Glottolog 4.8, 2023.
  * [10] Parvaneh Joharinad and Jürgen Jost.  Mathematical Principles of Topological and Geometric Data Analysis.  01 2023.
  * [11] Facundo Mémoli, Zhengchao Wan, and Yusu Wang.  Persistent laplacians: Properties, algorithms and implications.  SIAM Journal on Mathematics of Data Science, 4(2):858–884, 2022\.
  * [12] Lev Michael.  The classification of south american languages.  Annual Review of Linguistics, 7:329–349, 2021.
  * [13] Pieter Muysken.  Substrate influence in northern quechua languages.  In Enoch O. Aboh and Cécile B. Vigouroux, editors, Variation Rolls the Dice: A Worldwide Collage in Honour of Salikoko S. Mufwene, pages 133–160. University Press of Colorado, Boulder, 2021.
  * [14] O. Nenadic and M. Greenacre.  Correspondence analysis in r, with two- and three-dimensional graphics: The ca package.  Journal of Statistical Software, 20(3):1–13, 2007.
  * [15] Andrey Nikulin.  Proto-Macro-Jê: um estudo reconstrutivo.  PhD thesis, Universidade de Brasilia, 2020.
  * [16] Sietze Norder, Laura Becker, Hedvig Skirgård, Leonardo Arias, Alena Witzlack-Makarevich, and Rik van Gijn.  glottospace: R package for the geospatial analysis of linguistic and cultural data.  Journal of Open Source Software, 7(77):4303, 2022.
  * [17] Alexander Port, Iulia Gheorghita, Daniel Guth, John M. Clark, Crystal Liang, Shival Dasu, and Matilde Marcolli.  Persistent topology of syntax.  Math. Comput. Sci., 12(1):33–50, 2018.
  * [18] Alexander Port, Taelin Karidi, and Matilde Marcolli.  Topological analysis of syntactic structures.  Math. Comput. Sci., 16(1):Paper No. 2, 68, 2022.
  * [19] Andrew Robinson and Katharine Turner.  Hypothesis testing for topological data analysis.  J. Appl. Comput. Topol., 1(2):241–261, 2017.
  * [20] Aryon D. Rodrigues.  Macro-jê.  In R.M.W. Dixon and Alexandra Y. Aikhenvald, editors, The Amazonian languages, pages 164–206. Cambridge University Press, Cambridge, 1999\.
  * [21] Hedvig Skirgård, Hannah J. Haynie, Damián E. Blasi, Harald Hammarström, Jeremy Collins, Jay J. Latarche, Jakob Lesage, Tobias Weber, Alena Witzlack-Makarevich, Sam Passmore, Angela Chira, Luke Maurits, Russell Dinnage, Michael Dunn, Ger Reesink, Ruth Singer, Claire Bowern, Patience Epps, Jane Hill, Outi Vesakoski, Martine Robbeets, Noor Karolin Abbas, Daniel Auer, Nancy A. Bakker, Giulia Barbos, Robert D. Borges, Swintha Danielsen, Luise Dorenbusch, Ella Dorn, John Elliott, Giada Falcone, Jana Fischer, Yustinus Ghanggo Ate, Hannah Gibson, Hans-Philipp Göbel, Jemima A. Goodall, Victoria Gruner, Andrew Harvey, Rebekah Hayes, Leonard Heer, Roberto E. Herrera Miranda, Nataliia Hübler, Biu Huntington-Rainey, Jessica K. Ivani, Marilen Johns, Erika Just, Eri Kashima, Carolina Kipf, Janina V. Klingenberg, Nikita König, Aikaterina Koti, Richard G. A. Kowalik, Olga Krasnoukhova, Nora L. M. Lindvall, Mandy Lorenzen, Hannah Lutzenberger, Tânia R. A. Martins, Celia Mata German, Suzanne van der Meer, Jaime Montoya Samamé, Michael Müller, Saliha Muradoglu, Kelsey Neely, Johanna Nickel, Miina Norvik, Cheryl Akinyi Oluoch, Jesse Peacock, India O. C. Pearey, Naomi Peck, Stephanie Petit, Sören Pieper, Mariana Poblete, Daniel Prestipino, Linda Raabe, Amna Raja, Janis Reimringer, Sydney C. Rey, Julia Rizaew, Eloisa Ruppert, Kim K. Salmon, Jill Sammet, Rhiannon Schembri, Lars Schlabbach, Frederick W. P. Schmidt, Amalia Skilton, Wikaliler Daniel Smith, Hilário de Sousa, Kristin Sverredal, Daniel Valle, Javier Vera, Judith Voß, Tim Witte, Henry Wu, Stephanie Yam, Jingting Ye, Maisie Yong, Tessa Yuditha, Roberto Zariquiey, Robert Forkel, Nicholas Evans, Stephen C. Levinson, Martin Haspelmath, Simon J. Greenhill, Quentin D. Atkinson, and Russell D. Gray.  Grambank reveals the importance of genealogical constraints on linguistic diversity and highlights the impact of language loss.  Science Advances, 9(16):eadg6175, 2023.
  * [22] Hedvig Skirgård, Hannah J. Haynie, Harald Hammarström, Damián E. Blasi, Jeremy Collins, Jay Latarche, Jakob Lesage, Tobias Weber, Alena Witzlack-Makarevich, Michael Dunn, Ger Reesink, Ruth Singer, Claire Bowern, Patience Epps, Jane Hill, Outi Vesakoski, Noor Karolin Abbas, Sunny Ananth, Daniel Auer, Nancy A. Bakker, Giulia Barbos, Anina Bolls, Robert D. Borges, Mitchell Browen, Lennart Chevallier, Swintha Danielsen, Sinoël Dohlen, Luise Dorenbusch, Ella Dorn, Marie Duhamel, Farah El Haj Ali, John Elliott, Giada Falcone, Anna-Maria Fehn, Jana Fischer, Yustinus Ghanggo Ate, Hannah Gibson, Hans-Philipp Göbel, Jemima A. Goodall, Victoria Gruner, Andrew Harvey, Rebekah Hayes, Leonard Heer, Roberto E. Herrera Miranda, Nataliia Hübler, Biu H. Huntington-Rainey, Guglielmo Inglese, Jessica K. Ivani, Marilen Johns, Erika Just, Ivan Kapitonov, Eri Kashima, Carolina Kipf, Janina V. Klingenberg, Nikita König, Aikaterina Koti, Richard G. A. Kowalik, Olga Krasnoukhova, Kate Lynn Lindsey, Nora L. M. Lindvall, Mandy Lorenzen, Hannah Lutzenberger, Alexandra Marley, Tânia R. A. Martins, Celia Mata German, Suzanne van der Meer, Jaime Montoya, Michael Müller, Saliha Murado$\breve{g}$lu, HunterGatherer, David Nash, Kelsey Neely, Johanna Nickel, Miina Norvik, Bruno Olsson, Cheryl Akinyi Oluoch, David Osgarby, Jesse Peacock, India O.C. Pearey, Naomi Peck, Jana Peter, Stephanie Petit, Sören Pieper, Mariana Poblete, Daniel Prestipino, Linda Raabe, Amna Raja, Janis Reimringer, Sydney C. Rey, Julia Rizaew, Eloisa Ruppert, Kim K. Salmon, Jill Sammet, Rhiannon Schembri, Lars Schlabbach, Frederick W. P. Schmidt, Dineke Schokkin, Jeff Siegel, Amalia Skilton, Hilário de Sousa, Kristin Sverredal, Daniel Valle, Javier Vera, Judith Voß, Daniel Wikalier Smith, Tim Witte, Henry Wu, Stephanie Yam, Jingting Ye, Maisie Yong, Tessa Yuditha, Roberto Zariquiey, Robert Forkel, Nicholas Evans, Stephen C. Levinson, Martin Haspelmath, Simon J. Greenhill, Quentin D. Atkinson, and Russell D. Gray.  Grambank v1.0, March 2023.  Dataset.
  * [23] Primoz Skraba and Katharine Turner.  Wasserstein stability for persistence diagrams, 2023.  arXiv:2006.16824.
  * [24] Katharine Turner, Sayan Mukherjee, and Doug M. Boyer.  Persistent homology transform for modeling shapes and surfaces.  Inf. Inference, 3(4):310–344, 2014.
  * [25] S. van Buuren.  Flexible Imputation of Missing Data. Second Edition.  Chapman & Hall/CRC Press, Boca Raton, FL, 2018.
  * [26] Stef van Buuren and Karin Groothuis-Oudshoorn.  mice: Multivariate imputation by chained equations in r.  Journal of Statistical Software, 45(3):1–67, 2011.
  * [27] Afra J. Zomorodian.  Topology for computing, volume 16 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, 2009.  Reprint of the 2005 original [MR2111929].