# From centre to centres:
polycentric structures in individual mobility

Rohit Sahasrabuddhe1, Renaud Lambiotte2,3 & Laura Alessandretti4,∗

1 Indian Institute of Science Education and Research, Pune, India
2 Mathematical Institute, University of Oxford, United Kingdom
3 Turing Institute, London, United Kingdom
4 Technical University of Denmark, Anker Engelundsvej 1, Kgs Lyngby, DK-2800,
Denmark

∗ To whom correspondence should be addressed<EMAIL_ADDRESS>

## Abstract

The availability of large-scale datasets collected via mobile phones has
opened up opportunities to study human mobility at an individual level. The
granular nature of these datasets calls for the design of summary statistics
that can be used to describe succinctly mobility patterns. In this work, we
show that the radius of gyration, a popular summary statistic to quantify the
extent of an individual’s whereabouts, suffers from a sensitivity to outliers,
and is incapable of capturing mobility organised around multiple centres. We
propose a natural generalisation of the radius of gyration to a polycentric
setting, as well as a novel metric to assess the quality of its description.
With these notions, we propose a method to identify the centres in an
individual’s mobility and apply it to two large mobility datasets with socio-
demographic features, showing that a polycentric description can capture
features that a monocentric model is incapable of.

The availability of individual-level mobility data collected from mobile
phones has taken the study of human travel to a new direction, leading to
fundamental discoveries on the mechanisms driving spatial behaviour [1, 2, 3,
4], and enhancing our ability to address societal challenge, from forecasting
the spread of epidemics [5, 6, 7], to designing sustainable transport
solutions [8, 9, 10]. While opening a wide range of new opportunities, the
need to store, process and anonymise large scale individual trajectory
datasets is also raising new challenges, and calling for the development of
novel analytical and computational solutions [11, 12, 13].

One important issue relates to designing summary statistics that, given an
individual sequence of raw spatial coordinates, capture relevant aspects of
spatial behaviour, such as the regularity and predictability of a person’s
movements [4, 14], the evolution of visitation patterns over long time-scales
[15] or the interplay between exploration and exploitation of locations [15].
It is critical to identify measures that provide concise and accurate
descriptions of these aspects, but still allow to capture the heterogeneity of
behaviours across individuals.

In this work, we address this challenge focusing on one of the most widely
studied aspects of spatial behaviour, the area characterising individuals’
whereabouts, e.g. the region within which people move during the course of
their daily activities. Following the geographical literature, we refer to
this concept as to the _activity space_ [16, 17, 18]. Note that similar ideas
are investigated in animal ecology, where researchers refer to the _home-
range_ as to the region that encompasses the resources an animal requires to
survive and reproduce [19, 20, 21].

A number of approaches were proposed to operationalise the concept of
_activity space_ [18, 21]. In recent years, the _radius of gyration_ approach
inspired by the physics of rigid bodies has gained widespread popularity in
the Human Mobility literature, because of its simple formulation and scalable
computation time [3, 22]. This model, which we call the monocentric model,
proposes that, for each individual, the activity space can be summarised as a
circle extending on the earth surface. The centre of the circle, the so-called
_centre of mass_ , and its radius, the so-called _radius of gyration_ , are
determined by the positions of the locations visited by the individual and
their importance in terms of visitation probability. Following its
introduction in the pioneer article by Gonzalez et al [3], the monocentric
model has been used for many applications, including characterising
regularities across individuals [23, 24, 14, 4, 25, 26], modelling the spread
of epidemics [5, 27], studying patterns of crime activity [28, 29], and
movements during natural disasters [30, 31].

While being simple and concise, the description of individual whereabouts as a
circle has recently been challenged by empirical findings revealing that
mobility is organised around multiple centres. Individual daily movements are
anchored around several focal points, including home, work, and relevant
points of interests such as schools, family homes etc [32, 33] It has further
been shown that the organisation of infrastructure, from cities to locations
within cities, is polycentric [34] and hierarchical [35], and this structure
affects human travel [2]. Due to its monocentric nature, the _radius of
gyration_ approach for describing activity spaces may, in some cases, fail to
accurately capture spatial behaviour. We identify at least two anomalies: (i)
the _position_ of the centre of mass may be far from a person’s actual
whereabouts, because, when locations are organised around multiple focal
points, the centre of mass may be located in an unfamiliar area (see Fig. 1a);
and (ii) comparing the _extent_ of the activity space, measured by the radius
of gyration, for individuals with polycentric mobility may be misleading. For
example, individuals with the same radius of gyration but different number of
centres have widely different spatial behaviour (see Fig. 1b for an
illustration).

Figure 1: Inadequacy of the monocentric description of activity spaces. Each
map shows an individual’s activity space displayed using a log-azimuthal
projection [36]. Individuals (all from dataset D1, see data description in the
_Methods_ section) are based in the Zealand region of Denmark. (a) Position of
the centre of mass. In all examples (left to right), the centre of mass (red
star) lies far from any of the visited locations (black dots). Our method
identifies multiple centres (white triangles) closer to individual activities.
(b) Extent of the activity space. The examples (left to right) show the
activity spaces for three individuals with similar radius of gyration (red
circles), but different number of centres, according to our method (white
triangles). The radius of gyration is unable to capture the wide differences
in the extent of whereabouts across the three individuals.

In this work, we show that a simple generalisation of the radius of gyration
description can capture the polycentricity of human mobility. We propose a new
metric, the _cover score_ , that allows us to assess the quality of a mobility
description, by quantifying the interplay between its extent and the time
spent within it. In doing so, we distinguish between two different spatial
scales in human mobility, one associated to the distance between centres, and
one to the distance around the centres. We propose a new algorithm to compute
polycentric activity spaces, and we validate it against synthetic data. We
define the description arising from our method as the _flexicentric model_ ,
because it accommodates both monocentric and polycentric activity spaces. We
further validate our approach using real-world data: first, we formulate a
simple time allocation model that, given an activity space, describes where
individuals spend time; then, we show that, under this simple time allocation
model, the flexicentric description outperforms the monocentric model at the
task of predicting unseen data without overfitting. Finally, we use the
proposed method to extensively characterise activity spaces across two large-
scale datasets describing the GPS trajectories of $\sim{50,000}$ individuals
(see the _Methods_ section), in relation to socio-demographic attributes.

## Results

From monocentric to polycentric description of individual whereabouts.
Consider the set $L=\\{l_{1},...,l_{i},...l_{N}\\}$ of $N$ locations visited
by an individual within a certain period of time, where each visited location
$l_{i}$ is characterised by its geographical position $\mathbf{x_{i}}$, and
its weight $w_{i}$, defined as the fraction of time spent in location $l_{i}$
within the period considered, such that $\sum_{i=1}^{N}w_{i}=1$. In the
monocentric model, one summarises the whereabouts of this individual as a
circle centred at the so-called _centre of mass_ :

$\mathbf{c.o.m.}=\sum_{i=1}^{N}w_{i}{\mathbf{x_{i}}},$ (1)

with radius equal to the so-called _radius of gyration_ :

$r_{g}=\sqrt{\sum_{i=1}^{N}w_{i}\;(\mathbf{x_{i}}-\mathbf{c.o.m.})^{2}}.$ (2)

The origin of this quantity is rooted in the physics of rigid bodies. In their
pioneering article[3], the authors drew an analogy between a 2-dimensional
rigid body and the set of individual locations $L$ described above. In
classical mechanics, a rigid body can be described as a set of point masses,
where each point is characterised by its position and mass [37]. Similarly, in
the study of mobility, each location $l_{i}$ is characterised by its position
$\mathbf{x_{i}}$ and its weight $w_{i}$. Further, the _moment of inertia_ of a
rigid body quantifies the spread of its mass around a given axis, and
physically relates to how difficult it is to rotate the body about the axis.
For a 2-dimensional body, the moment of inertia can be defined with respect to
a centre $\mathbf{c}$, where the axis of rotation is perpendicular to the
plane of the body, and intersects it at $\mathbf{c}$. Using our notation, the
moment of inertia is defined as:

$I(\mathbf{c})=\sum_{i=1}^{N}w_{i}\;(\mathbf{x_{i}}-\mathbf{c})^{2}.$ (3)

It can be proven [37] that, choosing
$\mathbf{c}=\sum_{i=1}^{N}w_{i}\mathbf{x_{i}}$ (the _centre of mass_)
minimises the moment of inertia. The corresponding ’characteristic distance’
is the so-called _radius of gyration_ ,
$r_{g}=\sqrt{\frac{I(\mathbf{c.o.m.})}{\sum_{i=1}^{N}w_{i}}}$ (see Eq. 2), and
describes a circle around the centre of mass, such that if the entire mass
were concentrated upon it, the mechanical properties of the body would be
preserved. Note that the radius of gyration also has a simple statistical
interpretation, as it is the standard deviation of the data points when
$w_{i}$ is interpreted as their frequency in the Euclidean space.

Building on this analogy, the quantities described in Eq. 1 and Eq. 2
implicitly model an individual’s activity space as a circle, and have been
used in a variety of contexts [23, 24, 14, 4, 25, 26, 5, 27, 28, 29, 30, 31].
Yet, recent observational studies have documented polycentricity in human
mobility, which an approach identifying only one centre and one spatial scale
fails to capture [34, 35, 2, 32]. To address this important limitation, the
monocentric model can be extended to account for the polycentric nature of
individual activity spaces. Specifically, we propose that an individual’s
activity space can be described as a set of $k$ ($k\geq 1$) circles, with
centres $\\{\mathbf{c_{1}},\mathbf{c_{2}},\dots\mathbf{c_{k}}\\}$. Under this
flexicentric model, the moment of inertia (see Eq. 3) can be generalised as:

$I(\\{\mathbf{c}_{1},\mathbf{c}_{2},\dots\mathbf{c}_{k}\\})=\sqrt{\sum_{i=1}^{N}w_{i}\,\left(\min_{\alpha\in\\{1,2,\dots
k\\}}|\mathbf{x_{i}}-\mathbf{c}_{\alpha}|\right)^{2}}$ (4)

Using this new definition, the inertia quantifies the weighted average squared
distance of each location from its closest centre. Clearly, setting $k=1$
reduces the formula above to the monocentric description of inertia (see Eq.
3). Analogous to the monocentric model, we can choose the best set of centres
as those which minimise the generalised inertia (Eq. 4).

Identifying centres. We note that the quantity defined in Eq. 4 corresponds
precisely to the loss function of the k-means clustering algorithm [38]. The
k-means algorithm aims to partition a set of weighted points into $k$ clusters
in such a way that each point belongs to the cluster with the nearest weighted
mean (see also the _Methods_ section). In our case, the set of points
corresponds to the set of locations $L$ defined in the previous subsection.
Given this correspondence, we could simply use the k-means algorithm to
identify the centres minimising the inertia defined in Eq. 4.

However, it is known that the solution of k-means clustering is sensitive to
outliers [39], which are quite common in our case due to the broad
distribution of the location weights [3, 4]. Hence, we suggest a modified
version of the k-means algorithm, the trimmed-k-means (t-k-means), which
accounts for the presence of outliers by trimming out the furthest points from
the identified centres. Note that similar improvements of the k-means
algorithm have been suggested by existing literature [40][41][42](see the
description of the algorithm in the _Methods_ section). While accounting for
outliers, the t-k-means algorithm does not minimise the inertia in Eq. 4.

A novel metric to assess the quality of a model. The t-k-means algorithm
provides us a way to find centres for a given value of $k$, but what is the
optimal value of $k$? This question is notoriously challenging and generally
arises in clustering problems [43, 44]. In this section, we introduce the
_cover score_ \- a new metric specifically designed for assessing the quality
of a description of an individual’s activity space.

Having identified the centres for a given choice of $k$, we describe
individual’s whereabouts as a set of $k$ circles around these points. For
simplicity, all the circles have the same radius $r_{k}$. Here, we note that
the choice of $r_{k}$ is arbitrary, but should allow for comparisons across
choices of $k$. Thus, we choose $r_{k}$ such that the sum of the areas of the
$k$ circles is equal to the area of the circle described by the radius of
gyration. By setting $\pi kr_{k}^{2}=\pi r_{g}^{2}$, we find
$r_{k}=\frac{r_{g}}{\sqrt{k}}$.

The cover score quantifies the interplay between the area of the circles and
the time spent within the space defined by the circles, as we increase their
size. It is computed for a given set of centres as follows:

  * •

Consider a set $C(r)$ of $k$ circles located at the centres, where each circle
has radius $r$. Set $r=0$.

  * •

Increase the radius $r$ of the circles at the same rate, from $r=0$ up to
$r=r_{k}$. As the circles grow, compute:

    * (i)

the fraction of time spent in locations located inside any of the circles
$T(r)=\sum_{i\in\\{j\,|\,l_{j}\in C(r)\\}}w_{i}$, where $w_{i}$ are the
weights of the locations located within any of the circles in $C(r)$ (in our
notation, $l_{i}\in C(r)$)

    * (ii)

the ratio $A(r)=\frac{r^{2}}{r_{g}^{2}}$ between the total area covered by the
set of circles $C(r)$, and the maximum area, obtained when the circles have
radius $r=r_{k}$.

  * •

Compute the cover score as the area under the curve $T(r)$ vs $A(r)$, i.e. the
definite integral of the curve describing the growth of time spent within the
activity space as its size increases (see Fig. 2)

Figure 2: Illustration of the cover score on synthetic data. In the top panels
we show the distribution of locations in space (black dots), the centres (star
markers), and the corresponding circles with radius increasing from $0$ (blue
circles, from dark to lighter) up to the radius of gyration (dashed red line).
In the bottom panels, we show the curve describing the fraction of time spent
inside the circles, $T(r)$ as a function of the fraction of the total area
within the circles, $A(r)$. The area under the curve (coloured in blue)
corresponds to the cover score. Results are shown in (a) for 100 locations
sampled from a 2-dimensional Gaussian distribution with $\sigma^{2}=4$ along
both axes; and in (b) and (c) for 100 locations sampled (50 each) from two
Gaussian distributions with $\sigma^{2}=4$ along both axes, centred at $(0,0)$
and $(20,0)$. A monocentric representation is presented in (b), as well as a
polycentric ($k=2$) one in (c). For simplicity, all locations have the same
weight.

Intuitively, given an activity space with $k$ centres, the cover score is high
when activity is concentrated close to the centres, and low if the activity is
concentrated far from the centres. In Fig. 2, we illustrate how the cover
score is computed using synthetic data. Considering a set of $100$ locations
sampled from a 2-dimensional Gaussian distribution, we find that the
monocentric description provides a good representation of the data, with a
cover score of $0.358$ (Fig. 2a). Instead, considering locations generated
from two Gaussian distributions, we find that the monocentric model is not
suitable (Fig. 2b, top panel), as we achieve a low cover score of $0.133$
(Fig. 2b, bottom panel), while using two centres to describe the same data
(Fig. 2c, top panel), results in a much higher score of $0.876$ (Fig. 2c,
bottom panel).

Importantly, using empirical data, we find that the centres identified by t-k-
means achieve a higher cover score compared to those identified by the simple
k-means algorithm (see Fig. 3 and Supplementary Information Fig. S2).

Choosing the number of centres. The cover score can be used to compare the
quality of the solutions obtained by the t-k-means algorithm for different
values of $k$. We note that increasing the number of centres increases the
complexity of the model. Thus, we expect the quality of description to
increase as well. In the extreme situation when $k=|L|$, the cover score is
equal to $1$ (see Supplementary Information for details).

Since the problem we face is similar to the problem of choosing the value of
$k$ in the k-means clustering, we adopt a widely used solution to the latter,
the _gap statistic_ method [45]. This method allows us to evaluate if the
increase in cover score for larger values of $k$ is significantly higher than
expected simply due to the increase in model complexity. The method works as
follows. Let $x_{k}$ be the cover score of the clustering with $k$ centres. To
provide a baseline, we generate $100$ ’reference’ datasets with the same set
of weights, each by sampling uniformly within the bounding box of the original
data. Let $y_{k}$ and $sd_{k}$ be the mean and standard deviation respectively
of these reference datasets’ cover scores. We define the gap statistic as
$gap(k)=x_{k}-y_{k}$. To find the optimal value of $k$, we first set an upper
bound (6 in our simulations), and choose the maximum value of $k^{*}$ such
that its gap statistic is a significant improvement over the highest gap
statistic of $k<k^{*}$. That is, $k^{*}$ is the maximum value of $k$ such that
$gap(k^{*})-sd_{k^{*}}>gap(k^{\prime})$, where
$k^{\prime}=\operatorname*{argmax}_{k<k^{*}}gap(k)$. Using synthetic data, we
verified that the gap-statistic method is reliable when clusters are
sufficiently disjoint, and is also sensitive to outliers (see Supplementary
Information Fig. S4).

Figure 3: Polycentric activity spaces in empirical data. (a) t-k-means
outperforms k-means. Distribution of the cover score for centres obtained
using the t-k-means vs. the k-means methods for $k=k^{*}$. The inset depicts
the same distribution restricted to polycentric individuals. (b) The
flexicentric description outperforms the monocentric description. Distribution
of cover score of centres obtained setting $k=k^{*}$ vs $k=1$ using the t-k-
means algorithm. Elements on the diagonal correspond to the case when
$k^{*}=1$. The inset depicts the same distribution restricted to polycentric
individuals. (c) Distribution of the optimal number of centers $k^{*}$.
Results are shown for dataset D2 (see Supplementary Information Fig. S2 for
dataset D1).

Empirical Results. We apply the method described above to the trajectories of
$\sim 50,000$ individuals from two datasets: $D1$ describes the movements of
$566$ students at the Technical University of Denmark, and $D2$ describes the
movements of $\sim 49,000$ individuals from Germany, Japan and the United
Kingdom (see further details in the _Methods_ section). Using the gap
statistic method, we find that $42.8\%$ and $22.3\%$ of individuals are
polycentric ($k^{*}\geq 2$) in D1 and D2 respectively (see Supplementary
Information Table S3). Further, the centres identified by t-k-means have
larger cover score than those identified by the standard k-means algorithm:
t-k-means strictly outperforms k-means for 43921 (90.5%) of the individuals,
including 8516 (78.6%) out of 10828 polycentric individuals (see Fig. 3 for
results on D2, and Supplementary Information Fig. S2 for results on D1).

Validation with empirical data. Given a set of locations with corresponding
weights, the algorithm described above identifies the number and position of
centres of activity. But does the flexicentric activity space description
capture individual whereabouts without overfitting? We validate our approach
by comparing the ability of the monocentric and the flexicentric descriptions,
as estimated using only a fraction of an individual’s mobility trace, to
describe unseen data points.

In line with previous literature [46, 47], we use a simple Gaussian mixture to
model how individuals allocate time in the activity space. In particular, we
propose that, given an activity space characterised by a set of $k$ centres
$\\{\mathbf{c}_{1},\mathbf{c}_{2},\dots\mathbf{c}_{k}\\}$ and radius of
gyration $r_{g}$, the occupation probability $P(x)$, characterising the
likelihood of finding the individual in a given location $x$ is described by a
superposition of $k$ 2-dimensional symmetric Gaussians. Each Gaussian is
centred in one of the centres and has equal variance of
$\sigma^{2}=(\frac{r_{g}}{\sqrt{k}})^{2}$. To get an accurate description, we
propose that each Gaussian distribution has weight equal to the total weight
of the locations associated with its centre. Thus, we have that the occupation
probability is:

$P(x)=\sum_{\alpha}W_{\alpha}\mathscr{N}(\mathbf{c}_{\alpha},\,\sigma^{2})(x)\
$ (5)

, where $\mathscr{N}(\mathbf{c}_{\alpha},\,\sigma^{2})$ is a Gaussian
distribution with parameters $\mathbf{c}_{\alpha}$ and $\sigma^{2}$, and
$W_{\alpha}$ is the weight associated to centre $\mathbf{c}_{\alpha}$. Note
that the model reduces to a simple monocentric Gaussian model for $k=1$.

Using dataset D1, for each individual, we select a random sample consisting of
$80\%$ of all records (the train set), to find: (i) the centre of mass and
radius of gyration characterising the monocentric description (see Eq. 1 and
Eq. 2); (ii) the centres (with associated weights) characterising our
flexicentric description using the t-k-means algorithm. In this case, for
consistency across replicates of the routine, we set the number of cluster
$k=k^{*}$ as found by the gap-statistic method using the entire dataset. By
doing a train-test split on each individual 20 times, we verified that in
$(91.7\pm 17.5)\%$ of the cases, $k^{*}$ of the train set matches that of the
entire dataset. We use the values found in (i) and (ii) as input parameters to
the occupation probability in Eq. 5 for the monocentric and flexicentric
model, respectively. Note that, for individuals identified as monocentric
($k^{*}=1$) by our method, the two competing models differ only in the
position of the centre of the Gaussian, because the t-k-means algorithm does
not necessarily identify the centre of mass as the centre.

We evaluate the goodness of the two models at the task of predicting
individuals positions in space by computing the log-likelihood [48] of the
remaining $20\%$ of the data (the test set) under the two models. We repeat
this procedure $20$ times for each individual, and use the Welch’s t-test [49]
(at significance value $\alpha=0.01$) to compare the log-likelihoods of the
two competing predictions. We find that, out of $242$ cases of individuals
with $k^{*}>1$, the flexicentric model performs better than the monocentric
model in $230$ cases. For $323$ out of $324$ individuals with $k^{*}=1$, the
test is not conclusive, while in one case the monocentric model performs
better (see Supplementary Table S4 and Fig. 4). To check that our approach is
not overfitting, we compare the log-likelihoods of the train and test data
under the flexicentric model, and find that they differ significantly only for
37 out of 566 individuals (see also Fig. 4b). As a further validation against
overfitting, we compare the Shannon entropy describing the allocation of
weight across clusters $\left(-\sum_{\alpha}W_{\alpha}\log W_{\alpha}\right)$
in the train and test sets. We find that the average ratio between the test
and train entropies across polycentric ($k^{*}\geq 2$) users is $0.99\pm
0.05$, suggesting that the allocation of time is relatively similar in the two
sets (see also Supplementary Information Fig. S5). Altogether, the results
above reveal that, for individuals with polycentric activity, the flexicentric
model captures well individual whereabouts without overfitting.

Figure 4: Polycentric activity spaces describe individuals whereabouts without
overfitting. (a) Log-likelihood of the flexicentric model vs the monocentric
model. Results are shown for data in the test set (larger subplot) and train
set (inset). (b) Log-likelihood of the data in the test vs the train set under
the flexicentric model. In both subplots, only points such that the difference
in the two log-likelihoods is significant under Welch’s t-test (significance
level $\alpha=0.01$) are depicted in colour (monocentric individuals in blue,
and polycentric in orange). The rest are plotted in grey.

Understanding heterogeneities of individuals’ activity spaces. Having
identified that there is wide heterogeneity across individuals with respect to
the number of activity centres (see Supplementary Information Table S3), it is
natural to wonder which aspects are related to polycentric mobility. In this
section, we use dataset D2 to study the relation between the number of centres
and individuals socio-demographic attributes, focusing on gender, age, country
of origin, and urbanisation level, where the latter is estimated based on
individuals home-locations (see the _Methods_ section).

First, we focus on gender alone. We find that there is a statistically
significant difference between females and males, with $\sim 24.2\%$ of males
being polycentric, compared to only $19.1\%$ of females ($p\sim 10^{-38}$
under Fisher exact test [50]). This result reveals that gender and
polycentricity are correlated, and it is robust across the three countries
under study (see Table 1).

Next, we turn to urbanisation, where, for simplicity, we group individuals
into _urban_ and _non-urban_ (see the _Methods_ section). We find no
significant difference between the two groups, with polycentric individuals
representing $\sim 22.7\%$ of non-urban, and $\sim 22.2\%$ of urban
individuals, respectively ($p=0.212$ under the Fisher exact test[50]). The
result reveals that level of urbanisation and polycentricity are not
correlated, and it is verified in each of the three countries under study,
excluding Germany, where the correlation between the two is significant (see
Table 1).

To further investigate the correlation between gender and polycentricity, we
control for age. We divide the population into three age-groups: 18 to 26, 27
to 50, and 51+ years old. We find that in UK and Germany, the correlation
between polycentricity and gender is significant only in the middle age group
(27 to 50 years old). Instead, in Japan there is a strong correlation for all
age groups (see Table 2).

country | % polycentric | p-value | % polycentric | p-value
---|---|---|---|---
| males | females |  | urban | non-urban |
Japan | 24.4 | 19.2 | 6.46e-30** | 22.5 | 22.9 | 0.555
UK | 22.6 | 19.2 | 5.31e-4** | 20.8 | 21.7 | 0.369
Germany | 24.4 | 18.8 | 3.31e-7** | 20.1 | 23.1 | 0.00628**
Overall | 24.2 | 19.1 | 1.78e-38** | 22.2 | 22.7 | 0.212
Table 1: Correlation between urbanisation level, gender and the number of centres. Percentage of polycentric individuals by country and gender (left table) or by country and urbanisation-level (right table). P-values smaller than 0.01 (marked with **) indicate that there is a significant correlation between gender (left) or urbanisation (right) and polycentricity, under the Fisher exact test [50]. country | age group | % polycentric | p-value
---|---|---|---
|  | males | females |
Japan | 18-26 | 23.7 | 18.5 | 9.01e-4**
| 27-50 | 23.7 | 18.9 | 1.68e-16**
| 51 + | 25.9 | 19.9 | 1.15e-11**
UK | 18-26 | 19.9 | 18.6 | 0.566
| 27-50 | 22.6 | 18.8 | 0.00257**
| 51 + | 24.8 | 20.4 | 0.0888
Germany | 18-26 | 25.3 | 19.2 | 0.0139
| 27-50 | 23.6 | 18.2 | 8.69e-5**
| 51 + | 26.4 | 21.3 | 0.0815
Table 2: Correlation between age and number of centres. Percentage of
polycentric individuals by country, age group and gender. P-values smaller
than 0.01 (marked with **) indicate that there is a significant correlation
between gender and polycentricity (given country and age group), under the
Fisher exact test [50].

## Discussion

The widely popular radius of gyration [3] for quantifying activity spaces has
recently been challenged by new empirical findings showing that locations
visited by individuals are organised in multiple clusters [32, 2]. In this
paper, we have proposed a generalisation of the radius of gyration, the
flexicentric model, that properly accounts for the polycentric nature of human
activities. We have also introduced the _cover score_ , a new metric to assess
the quality of a polycentric model describing individual whereabouts, and
presented a new algorithm to compute individual numbers of centres and their
locations. Analysing two real-world large-scale datasets of GPS trajectories,
we have found that around $42.8\%$ and $22.3\%$ of individuals in our two
datasets respectively have polycentric mobility behaviour. We showed that, for
these individuals, the flexicentric model trained on a subset of the entire
trajectory allows to better describe unseen data in comparison to the
monocentric model, ensuring that it does not overfit. Interestingly, we
discovered that having more than one centre is significantly more likely in
men compared to women, especially in the middle age group, across the
different countries under study. This result is in line with similar empirical
findings on the existence of gender gaps in human mobility [51, 2]. A
potential reason why the gender gap is more pronounced in the middle age-group
concerns the fact that employed individuals are likely to have at least two
centres, if their home and work locations are far enough [32]. The gender gap
in employment may thus partly explained the differences in travel behaviour
across genders [52]. In contrast, we found that the level of urbanisation
around an individual’s home location does not relate to the number of centres.
These observations would not have been possible from a monocentric point of
view. By opening the possibility to characterise individuals whereabouts in a
polycentric fashion, our work allows to enrich the description of human
mobility, going from a single quantity, i.e. the radius of gyration, to
multiple numbers, or spatial scales, associated to distances between centres
and around the centres. This new lens to study individual-level travel
behaviour may help to improve our quantification of heterogeneities across
individuals, which is critical to design accurate mobility models, and
ultimately improve our ability to forecast the spread of epidemics and plan
transportation that account for the diverse needs of all individuals. In
addition to these potential applications, future works include the possibility
to study the evolution of individuals whereabouts over time and investigate
human mobility in more general metric spaces, where the distance between
locations is not purely geographical, but may be built on commuting time
information or distance in transport networks.

## Methods

### Data Description

Our study is based on two datasets.

Dataset D1 was collected as part of an experiment that took place between
September 2013 and September 2015[53]. The experiment involved $851$ Technical
University of Denmark students (about 22% female and about 78% male),
typically aged between $19$ and $21$ years old. Participants’ GPS records over
time was sampled every 5 minutes. Data collection was approved by the Danish
Data Protection Agency. All participants provided informed consent. We
selected 5 months of data for $566$ individuals from the entire dataset (see
Supplementary Information Section S1).

Dataset D2 contains anonymised GPS location data collected by a global
smartphone and electronics company between 2017 and 2019. The data consist of
anonymised users with self-reported age, gender, height, weight and country of
residence. Data were extracted through a smartphone app. All data analysis was
carried out in accordance with the European Union’s General Data Protection
Regulation 2016/679 (GDPR) and the regulations set out by the Danish Data
Protection Agency. Individuals provided informed consent to the analysis of
their data. We estimated the level of urbanisation and country for each
individual from the position of the top-location. Data on the urbanisation
level in the area surrounding individuals’ home locations is based on the GHS
Settlement Model grid that delineates and classifies settlement typologies via
a logic of population size, population and built-up area densities. This
classification categorises areas in urban areas, towns and rural areas. In our
analysis, we merged towns and rural areas into a single category. Data can be
downloaded from: https://ghsl.jrc.ec.europa.eu/data.php. We selected
approximately 50,000 individuals from the top three countries in the dataset
(Japan, Germany and the UK) with at least six months of data, whose position
is known, every day, on average 75% of the time. Supplementary Tables S1 and
S2 summarise the distributions of urbanisation and gender. We categorise
individuals into three age groups: $18-26$ years old, $27-50$ years old, and
$\geq 51$ years old. The summary of the age distribution is included in
Supplementary Information Fig. S1.

### Data Availability

Derived data that support the findings of this study are available in DTU Data
with the identifier [AVAILABLE UPON PUBLICATION]. Additional data related to
this paper may be requested from the authors. Raw data for dataset D1 are not
publicly available due to privacy considerations, but are available to
researchers who meet the criteria for access to confidential data, sign a
confidentiality agreement and agree to work under supervision in Copenhagen.
Please direct your queries to the corresponding author. Raw data for dataset
D2 are not publicly available to preserve individuals’ privacy under the
European General Data Protection Regulation.

### Code Availability

Code is available at https://github.com/rohit-sahasrabuddhe/polycentric-
mobility.git.

### The trimmed k-means algorithm

The traditional k-means algorithm takes as input a set of (possibly weighted)
points embedded in a metric space and the number of clusters $k$ [38]. It
iteratively partitions the data points into $k$ clusters as follows. The
algorithm is initialised by choosing $k$ centres randomly amongst the data
points. In every iteration, points are assigned to their closest centre, and
each centre is updated to the (weighted) mean of the points assigned to it.
The algorithm terminates when the centres become stationary or when the number
of iterations exceeds a pre-defined limit. It is well-known that k-means
clustering is sensitive to the initial choice of centres, as well as to
outliers in the data. Two types of strategies to deal with these issues are
(a) choosing a better initialisation of centres (such as k-means++ [54]), and
(b) running the clustering several times with different initial centres and
choosing the best result [38].

The problem of outliers affects location data, because, due to the effects of
exploration and preferential return [4], the distribution of time spent per
location is highly skewed. In order to ensure that places where individuals
spend little time do not affect dramatically the solution of the clustering,
and inspired by previous modifications of the k-means algorithm [40], we
propose the trimmed k-means algorithm, _t-k-means_. At each iteration of t-k-
means, the centre is updated to the (weighted) mean of points corresponding to
$t$ fraction of weight closest to the centre. This reduces the effect of
points that are assigned to the cluster, but are far from the centre. We note
that setting the trimming parameter $t=1$ corresponds to running the
traditional k-means algorithm. The value of $t$ is set to $0.9$ for the
results presented in this paper, and are robust to variations of this
parameter (see Supplementary Information Fig. S3). Further, we deal with the
sensitivity to initialisation by clustering each dataset $50$ times, and then
picking the result that minimises the inertia.

Acknowledgements RS is supported by the DST-Inspire scholarship (Government of
India).
Author Contributions. RS, RL and LA designed the study. RS performed the
analyses and implemented the code. All authors analysed the results and wrote
the paper.
Competing Interests. The authors declare that they have no competing financial
interests.

Supplementary Information
From centre to centres: polycentric structures in individual mobility

##  S1 Data description and pre-processing

For both datasets, the raw data of mobility traces is first put through the
Infostop algorithm [55] with parameters $r_{1},r_{2}=30$ metres, $t_{min}=10$
mins. This gives us, for each individual, a set of records, each corresponding
to the individual staying within a 30 meter radius for at least 10 minutes.
For each record, we also have the latitude and longitude, and the amount of
time spent at that stop location.

Filtering of dataset D1. With this data, we perform the following steps to
clean the dataset D1 and choose a subset to analyse:

  * •

We first discard all records with time-spent $\geq$ 7 days.

  * •

We restrict the analysis to locations in the Sjælland and Hovedstaden
administrative regions (as per the shapefile obtained from https://diva-
gis.org/ accessed 20-Aug-2020). This region contains 88.8% of records
corresponding to 91.6% of the time spent.

  * •

We then consider records from the Spring semester (January - May) of 2014
which corresponds to 151 days. For each individual, we calculate the active
time as the total amount of time for which we know their stop location. We
note that the distribution of active times is bimodal, and thus restrict our
dataset to individuals of high active time by placing a cutoff at the 68 day
mark.

  * •

At the end, we arrive at a dataset of 566 individuals, each with more than 68
days worth of data in the region of interest.

Filtering of dataset D2. We perform the following steps to clean and filter
the dataset D2:

  * •

We discard all records with time-spent $\geq$ 7 days.

  * •

We select users whose average daily time-coverage (fraction of time an
individual’s location is known) larger than $75\%$

  * •

We select individuals whose country of origin is either Japan or Germany or
the UK. These are the top three countries by number of users in the entire
dataset.

  * •

For each individual, we select the first $180$ days of data.

For the users selected in dataset D2, we present the distribution of age in
Fig. S1, the distribution of urbanization level in Table S1 and the
distribution of gender in Table S2.

Figure S1: Age distribution, and classification into age groups of individuals in dataset D2. Country | Rural (1) | Semi-urban (2) | Urban (3) | Total
---|---|---|---|---
Japan | 1903 | 5061 | 28644 | 35608
UK | 918 | 1822 | 4289 | 7029
Germany | 1474 | 2075 | 2348 | 5897
Total | 4295 | 8958 | 35281 | 48534
Table S1: Number of individuals per country and urbanization-level in dataset D2. Country | Males | Females | Total
---|---|---|---
Japan | 23187 | 12421 | 35608
UK | 4091 | 2938 | 7029
Germany | 3280 | 2617 | 5897
Total | 30558 | 17976 | 48534
Table S2: Number of individuals per country and gender in dataset D2.

##  S2 The t-k-means algorithm

Empirical results on applying the t-k-means algorithm on dataset D1 are
presented in Fig. S2.
Considerations on the relation between the cover score and the number of
centres. Increasing the number of centres $k$ increases the complexity of the
model, and thus we expect an increase in the quality of description. However,
we note that the cover score can either increase or decrease when increasing
$k$. For example, creating a new centre in a sparse region could result in the
cover score decreasing, because (a) the area around the other centres is
reduced, and (b) the sparse area around the new centre is given more
importance. In this respect, there is a subtle but important difference with
respect to the k-means clustering algorithm, in which the inertia is used (see
Eq. 3) to assess the quality of a solution. The inertia, in fact, decreases
monotonically with increasing $k$, because adding a centre results necessarily
in a better description.

Figure S2: Results on D1. (a) t-k-means outperforms k-means Distribution of the cover score for centres obtained using the t-k-means vs. the k-means methods for $k=k^{*}$. The inset depicts the same distribution restricted to polycentric individuals. (b) The flexicentric description outperforms the monocentric description. Distribution of cover score of centres obtained setting $k=k^{*}$ vs $k=1$ using the t-k-means algorithm. The inset depicts the same distribution restricted to polycentric individuals. (c) Distribution of $k^{*}$. $k^{*}$ | number of individuals
---|---
| D1 | D2
1 | 324 | 37706
2 | 187 | 8634
3 | 43 | 1206
4 | 12 | 555
5 | 0 | 287
6 | 0 | 146
Table S3: Distribution of $k^{*}$ for empirical datasets D1 and D2.

Sensitivity to changes in trimming parameter $t$. We check the sensitivity of
the results to variation in $t$ by running t-k-means on all individuals in the
D1 dataset for $t\in[0.85,0.95]$, and call the value of $k^{*}$ as
$k^{*}_{t}$. Note that in this notation, the results in the main text
correspond to $k^{*}_{0.9}$. In Fig. S3, we plot the fraction of individuals
for whom the $k^{*}_{t}=k^{*}_{0.9}$.

Figure S3: Sensitivity of $k^{*}$ on the choice of trimming parameter $t$.

##  S3 Validation against synthetic data.

We first validate our approach for selecting $k^{*}$ using synthetic data. We
consider datasets of $100$ locations sampled from two bivariate Gaussian
distributions with variance $\sigma^{2}$ (along both axes) and distance
between their centres $d$, implying the number of centre is $k^{*}=2$. For
simplicity, we set all locations to have equal weights, but sample a different
number of points from each Gaussian. In particular, we sample points from the
50/50, 60/40, 70/30, and 80/20 splits. In each case, we vary $d$ and
$\sigma^{2}$ in the range $[1,10]$, and record the value of $k^{*}$ across
$25$ runs over randomly generated datasets for each pair. We find that the
value of $k^{*}$ identified by the algorithm has desirable properties (Fig
S4). First, for increasingly uneven splits, the algorithm is more likely to
predict $k^{*}=1$, implying the method is robust to outliers. Second, for each
split, the algorithm is more likely to predict $k^{*}=2$ when $d$ is high and
$\sigma^{2}$ is low, which is also desirable.

Figure S4: Extracting clusters from synthetic data: The plots depict the value
of $k^{*}$ identified by our algorithm for synthetic data comprising of points
sampled from 2 bivariate Gaussian distributions. The 4 panels show (from left
to right) the results for a 50/50, 60/40, 70/30, and 80/20 split. Each cell in
the heatmap corresponds to a particular value of $\sigma^{2}$ (variance of
each Gaussian along both axes) and $d$ (distance between their centres), and
depicts the average value of $k^{*}$ across 25 runs over randomly generated
datasets.

##  S4 Validation against empirical data

Results for the empirical validation of our method using a tran/test split are
presented here and summarized in Table S4.

Monocentric individuals: Of the 324 monocentric individuals, we find that the
traditional monocentric model outperforms our model for 84 individuals in the
train set, but just 1 individual in the test set. Our method never outperforms
the traditional method.

Polycentric individuals: Of the 242 polycentric individuals, our model
outperforms the traditional method in 237 and 230 individuals in the train and
test set respectively. In turn, the traditional method is significantly better
for 2 and 1 individual respectively.

Train vs test: Comparing the log-likelihoods of our model on the train and
test sets shows that their means differ statistically significantly only for
37 individuals out of 566. The log-likelihood of the train set is
significantly greater for 25 of these, and that of the test set is greater for
the remaining 12.

| Train set | Test set
---|---|---
| flexicentric | monocentric | flexicentric | monocentric
Monocentric (324) | 0 | 84 | 0 | 1
Polycentric (242) | 237 | 2 | 230 | 1
Table S4: Validation of the flexicentric method on D1. The values under the
’flexicentric’ (resp. ’monocentric’) correspond to the number of individuals
for whom the log-likelihood of data under the flexicentric (resp. monocentric)
model was significantly better. The significance of the difference is
established via Welch’s t-test, and the values in this table correspond to a
p-value $<$ 0.01. Figure S5: The distribution of the ratio of the Shannon
entropy of the weight distribution across clusters in the test and train sets
for polycentric users in D1. The mean is 0.99 $\pm$ 0.05.

## References

  * [1] Markus Schläpfer, Lei Dong, Kevin O’Keeffe, Paolo Santi, Michael Szell, Hadrien Salat, Samuel Anklesaria, Mohammad Vazifeh, Carlo Ratti, and Geoffrey B West.  The universal visitation law of human mobility.  Nature, 593(7860):522–527, 2021.
  * [2] Laura Alessandretti, Ulf Aslak, and Sune Lehmann.  The scales of human mobility.  Nature, 587(7834):402–407, 2020.
  * [3] Marta C Gonzalez, Cesar A Hidalgo, and Albert-Laszlo Barabasi.  Understanding individual human mobility patterns.  nature, 453(7196):779–782, 2008.
  * [4] Chaoming Song, Zehui Qu, Nicholas Blumm, and Albert-László Barabási.  Limits of predictability in human mobility.  Science, 327(5968):1018–1021, 2010.
  * [5] Emanuele Pepe, Paolo Bajardi, Laetitia Gauvin, Filippo Privitera, Brennan Lake, Ciro Cattuto, and Michele Tizzoni.  Covid-19 outbreak response, a dataset to assess mobility changes in italy following national lockdown.  Scientific data, 7(1):1–7, 2020.
  * [6] Ruth F Hunter, Leandro Garcia, Thiago Herick de Sa, Belen Zapata-Diomedi, Christopher Millett, James Woodcock, Esteban Moro, et al.  Effect of covid-19 response policies on walking behavior in us cities.  Nature Communications, 12(1):1–9, 2021.
  * [7] Serina Chang, Emma Pierson, Pang Wei Koh, Jaline Gerardin, Beth Redbird, David Grusky, and Jure Leskovec.  Mobility network models of covid-19 explain inequities and inform reopening.  Nature, 589(7840):82–87, 2021.
  * [8] Xuan Song, Hiroshi Kanasugi, and Ryosuke Shibasaki.  Deeptransport: Prediction and simulation of human mobility and transportation mode at a citywide level.  In Proceedings of the Twenty-Fifth International Joint Conference on Artificial Intelligence, pages 2618–2624, 2016.
  * [9] Sina Dabiri and Kevin Heaslip.  Inferring transportation modes from gps trajectories using a convolutional neural network.  Transportation research part C: emerging technologies, 86:360–371, 2018.
  * [10] Li Zhu, Fei Richard Yu, Yige Wang, Bin Ning, and Tao Tang.  Big data analytics in intelligent transportation systems: A survey.  IEEE Transactions on Intelligent Transportation Systems, 20(1):383–398, 2018.
  * [11] Hugo Barbosa, Marc Barthelemy, Gourab Ghoshal, Charlotte R James, Maxime Lenormand, Thomas Louail, Ronaldo Menezes, José J Ramasco, Filippo Simini, and Marcello Tomasini.  Human mobility: Models and applications.  Physics Reports, 734:1–74, 2018.
  * [12] Yves-Alexandre De Montjoye, Sébastien Gambs, Vincent Blondel, Geoffrey Canright, Nicolas De Cordes, Sébastien Deletaille, Kenth Engø-Monsen, Manuel Garcia-Herranz, Jake Kendall, Cameron Kerry, et al.  On the privacy-conscientious use of mobile phone data.  Scientific data, 5(1):1–6, 2018.
  * [13] David Lazer, Eszter Hargittai, Deen Freelon, Sandra Gonzalez-Bailon, Kevin Munger, Katherine Ognyanova, and Jason Radford.  Meaningful measures of human society in the twenty-first century.  Nature, pages 1–8, 2021.
  * [14] Xin Lu, Erik Wetter, Nita Bharti, Andrew J Tatem, and Linus Bengtsson.  Approaching the limit of predictability in human mobility.  Scientific reports, 3(1):1–9, 2013.
  * [15] Laura Alessandretti, Piotr Sapiezynski, Vedran Sekara, Sune Lehmann, and Andrea Baronchelli.  Evidence for a conserved quantity in human mobility.  Nature human behaviour, 2(7):485–491, 2018.
  * [16] Shannon N Zenk, Amy J Schulz, Stephen A Matthews, Angela Odoms-Young, JoEllen Wilbur, Lani Wegrzyn, Kevin Gibbs, Carol Braunschweig, and Carmen Stokes.  Activity space environment and dietary and physical activity behaviors: a pilot study.  Health & place, 17(5):1150–1161, 2011.
  * [17] David WS Wong and Shih-Lung Shaw.  Measuring segregation: An activity space approach.  Journal of geographical systems, 13(2):127–145, 2011.
  * [18] Jill E Sherman, John Spencer, John S Preisser, Wilbert M Gesler, and Thomas A Arcury.  A suite of methods for representing activity space in a healthcare accessibility study.  International journal of health geographics, 4(1):1–21, 2005.
  * [19] Don W Hayne.  Calculation of size of home range.  Journal of mammalogy, 30(1):1–18, 1949.
  * [20] Peter N Laver and Marcella J Kelly.  A critical review of home range studies.  The Journal of Wildlife Management, 72(1):290–298, 2008.
  * [21] Roger A Powell and Michael S Mitchell.  What is a home range?  Journal of mammalogy, 93(4):948–958, 2012.
  * [22] Chaoming Song, Tal Koren, Pu Wang, and Albert-László Barabási.  Modelling the scaling properties of human mobility.  Nature physics, 6(10):818–823, 2010.
  * [23] Bartosz Hawelka, Izabela Sitko, Euro Beinat, Stanislav Sobolevsky, Pavlos Kazakopoulos, and Carlo Ratti.  Geo-located twitter as proxy for global mobility patterns.  Cartography and Geographic Information Science, 41(3):260–271, 2014\.
  * [24] Luca Pappalardo, Filippo Simini, Salvatore Rinzivillo, Dino Pedreschi, Fosca Giannotti, and Albert-László Barabási.  Returners and explorers dichotomy in human mobility.  Nature communications, 6(1):1–8, 2015.
  * [25] Raja Jurdak, Kun Zhao, Jiajun Liu, Maurice AbouJaoude, Mark Cameron, and David Newth.  Understanding human mobility from twitter.  PloS one, 10(7):e0131469, 2015.
  * [26] Zhiyuan Cheng, James Caverlee, Kyumin Lee, and Daniel Sui.  Exploring millions of footprints in location sharing services.  In Proceedings of the International AAAI Conference on Web and Social Media, volume 5, 2011.
  * [27] Stefano Merler and Marco Ajelli.  The role of population heterogeneity and human mobility in the spread of pandemic influenza.  Proceedings of the Royal Society B: Biological Sciences, 277(1681):557–565, 2010.
  * [28] Wim Bernasco.  Adolescent offenders’ current whereabouts predict locations of their future crimes.  PloS one, 14(1):e0210733, 2019.
  * [29] Sophie Curtis-Ham, Wim Bernasco, Oleg N Medvedev, and Devon LL Polaschek.  A national examination of the spatial extent and similarity of offenders’ activity spaces using police data.  ISPRS International Journal of Geo-Information, 10(2):47, 2021.
  * [30] Qi Wang and John E Taylor.  Quantifying human mobility perturbation and resilience in hurricane sandy.  PLoS one, 9(11):e112608, 2014.
  * [31] Qi Wang and John E Taylor.  Patterns and limitations of urban human mobility resilience under the influence of multiple types of natural disaster.  PLoS one, 11(1):e0147299, 2016.
  * [32] James P Bagrow and Yu-Ru Lin.  Mesoscopic structure and social aspects of human mobility.  PloS one, 7(5):e37676, 2012.
  * [33] Christian M Schneider, Vitaly Belik, Thomas Couronné, Zbigniew Smoreda, and Marta C González.  Unravelling daily human mobility motifs.  Journal of The Royal Society Interface, 10(84):20130246, 2013.
  * [34] Rémi Louf and Marc Barthelemy.  Modeling the polycentric transition of cities.  Physical review letters, 111(19):198702, 2013.
  * [35] Aleix Bassolas, Hugo Barbosa-Filho, Brian Dickinson, Xerxes Dotiwalla, Paul Eastham, Riccardo Gallotti, Gourab Ghoshal, Bryant Gipson, Surendra A Hazarie, Henry Kautz, et al.  Hierarchical organization of urban mobility and its connection with city livability.  Nature communications, 10(1):1–10, 2019.
  * [36] John P Snyder.  “magnifying-glass” azimuthal map projections.  The American Cartographer, 14(1):61–68, 1987.
  * [37] Herbert Goldstein, Charles Poole, and John Safko.  Classical mechanics, 2002.
  * [38] Aristidis Likas, Nikos Vlassis, and Jakob J Verbeek.  The global k-means clustering algorithm.  Pattern recognition, 36(2):451–461, 2003.
  * [39] Ville Hautamäki, Svetlana Cherednichenko, Ismo Kärkkäinen, Tomi Kinnunen, and Pasi Fränti.  Improving k-means by outlier removal.  In Scandinavian Conference on Image Analysis, pages 978–987. Springer, 2005.
  * [40] Juan Antonio Cuesta-Albertos, Alfonso Gordaliza, Carlos Matrán, et al.  Trimmed $k$-means: An attempt to robustify quantizers.  Annals of Statistics, 25(2):553–576, 1997.
  * [41] Dau Pelleg and Andrew Moore.  X-means: Extending k-means with efficient estimation of the number of clusters.  In In Proceedings of the 17th International Conf. on Machine Learning, pages 727–734. Morgan Kaufmann, 2000.
  * [42] Partitioning Around Medoids (Program PAM), chapter 2, pages 68–125.  John Wiley & Sons, Ltd, 1990.
  * [43] Mark EJ Newman and Gesine Reinert.  Estimating the number of communities in a network.  Physical review letters, 117(7):078301, 2016.
  * [44] Greg Hamerly and Charles Elkan.  Learning the k in k-means.  Advances in neural information processing systems, 16:281–288, 2004\.
  * [45] Robert Tibshirani, Guenther Walther, and Trevor Hastie.  Estimating the number of clusters in a data set via the gap statistic.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(2):411–423, 2001.
  * [46] Christen H Fleming and Justin M Calabrese.  A new kernel density estimator for accurate home-range and species-range area estimation.  Methods in Ecology and Evolution, 8(5):571–579, 2017.
  * [47] Joni A Downs and Mark W Horner.  Effects of point pattern shape on home-range estimates.  The Journal of Wildlife Management, 72(8):1813–1818, 2008.
  * [48] Aaron Clauset, Cosma Rohilla Shalizi, and Mark EJ Newman.  Power-law distributions in empirical data.  SIAM review, 51(4):661–703, 2009.
  * [49] Bernard L Welch.  The generalization of ‘student’s’problem when several different population varlances are involved.  Biometrika, 34(1-2):28–35, 1947.
  * [50] Michel Raymond and François Rousset.  An exact test for population differentiation.  Evolution, pages 1280–1283, 1995.
  * [51] Laetitia Gauvin, Michele Tizzoni, Simone Piaggesi, Andrew Young, Natalia Adler, Stefaan Verhulst, Leo Ferres, and Ciro Cattuto.  Gender gaps in urban mobility.  Humanities and Social Sciences Communications, 7(1):1–13, 2020\.
  * [52] Wei-Shiuen Ng and Ashley Acker.  Understanding urban travel behaviour by gender for efficient and equitable transport policies.  International Transport Forum Discussion Paper, 2018.
  * [53] Arkadiusz Stopczynski, Vedran Sekara, Piotr Sapiezynski, Andrea Cuttone, Mette My Madsen, Jakob Eg Larsen, and Sune Lehmann.  Measuring large-scale social networks with high resolution.  PloS one, 9(4):e95978, 2014.
  * [54] David Arthur and Sergei Vassilvitskii.  k-means++: The advantages of careful seeding.  Technical report, Stanford, 2006.
  * [55] Ulf Aslak and Laura Alessandretti.  Infostop: scalable stop-location detection in multi-user mobility data.  arXiv preprint arXiv:2003.14370, 2020.