These authors contributed equally to this work.

[1]Paolo Girardi These authors contributed equally to this work.

These authors contributed equally to this work.

1]Dipartimento di Scienze Ambientali, Informatica e Statistica, Università Ca’
Foscari, Via Torino 155, Venezia-Mestre, 30172, Italy

2]Department of Pure and Applied Sciences, Kirinyaga University, C74, Kutus-
Kerugoya Rd, Kerugoya, 10300, Kenya

MSC Classification]91C20, 60G60, 62P12, 37M10, 62F15

# Spatial quantile clustering of climate data

Carlo Gaetan<EMAIL_ADDRESS><EMAIL_ADDRESS>Victor Muthama Musau
<EMAIL_ADDRESS>[  [

###### Abstract

In the era of climate change, the distribution of climate variables evolves
with changes not limited to the mean value. Consequently, clustering
algorithms based on central tendency could produce misleading results when
used to summarize spatial and/or temporal patterns. We present a novel
approach to spatial clustering of time series based on quantiles using a
Bayesian framework that incorporates a spatial dependence layer based on a
Markov random field. A series of simulations tested the proposal, then applied
to the sea surface temperature of the Mediterranean Sea, one of the first seas
to be affected by the effects of climate change.

###### keywords:

Asymmetric Laplace distribution, Markov random field, model-based clustering,
time series

###### pacs:

[

##  1 Introduction

The extraction of relevant information hidden in complex spatio-temporal
climate datasets is one of the main goals of statistical climatology.
Clustering algorithms are routinely used to summarize and visualize important
spatial and/or temporal patterns.

To identify spatial patterns, the most well-known statistical techniques are
based on the concept of intra- and inter-cluster variances, such as the
k-means algorithm or methods based on empirical orthogonal function analysis.

K-means is a widely used clustering procedure that identifies clusters by
minimizing the distance between each cluster’s representative center, called
centroid, and the observations it contains. The distance measure can be
defined using an appropriate norm. This approach is appropriate for
identifying patterns related to unknown mean behaviors. The principle is
ideally suited when the variables of interest follow a distribution that
spreads around a center, such as in the case of a mixture of Gaussian
distributions.

On the other hand, climate change refers to long-term shifts in environmental
conditions caused primarily by human activities [1]. This phenomenon
eventually leads to changes in the distribution of many climatic variables [2,
3, 4]. These changes emphasize the significance of utilizing suitable
statistical methods that take into account non-stationary patterns when
analysing time series observations with a flexible approach. Quantile
regression has gained popularity in detecting changes in various features of
climate variable distribution. This method enables handling of different
quantiles and the presence of outliers [5, 6, 7].

Several clustering methods for time series have also been proposed [8]. The
former methods are based on raw time series data and adapt distance/similarity
measures of multivariate data to time series data. Latter methods perform a
dimension reduction extracting some vectors of features or parameters. The
differences between these vectors can then be clustered using standard
clustering techniques. Dimension reduction is also used in clustering
functional data, where the time series is seen as a finite realization of a
time-varying random function at regular intervals [9].

Alternatively, we can use model-based clustering techniques [10, 11], but the
widespread application of the Gaussian distribution may affect the results in
the presence of asymmetry in the data distribution.

The presence of outliers in the time series can make it difficult to select
the appropriate clustering technique. The robustness of clustering algorithms
can be assessed through the application of trimming tools [12, 13, 14]. Their
main focus was to obtain robust clustering procedures (i.e. robust k-means)
developing mathematical methods that perform both classification and outliers
detection, as the introduction of restrictions on the dispersion matrices of
the groups to avoid detection of spurious clusters [15] or the inclusion of
constraints on the scatter matrices [16].

Identifying regions with common temporal patterns is one of the potential
objectives of clustering climate variables. However, this is complicated by
the spatial dimension. The literature proposes various workarounds to address
this issue.

In the context of model-based clustering, Nguyen et al. [17] modeled the
single time series as an autoregressive process, with a Markov random field
inducing a spatial structure for the cluster membership. Disegna et al. [18]
employed a copula-based clustering algorithm. Spatial time series can also be
treated as spatial functional data [19, 20, 21]. Still in the framework of
model-based clustering Jiang and Serban [22] utilized a Markov random field to
induce a spatial structure for the cluster memberships of time functions over
a site network. This idea has been extended in a Bayesian framework by Hu et
al. [23]. Instead, Secchi et al. [24] suggested a bagging-Voronoi strategy for
unsupervised classification of functional data. In a geostatistical context,
Giraldo et al. [25] represented time series as curves using basis functions
and weighted the dissimilarity matrix among the curves using variograms
calculated with the coefficients of the basis.

This paper introduces a mixture-model approach for time series clustering. The
focus is on identifying common long-term trends and seasonal patterns across a
collection of time series of sea surface temperature measurements on a grid
that spans the Mediterranean Sea. The dataset is particularly challenging due
to the skewness of the data distribution, time-varying heteroscedasticity, and
the presence of outliers. Our model considers time-varying quantiles, enabling
us to address the aforementioned issues and potentially attain more reliable
clustering outcomes.

To model a single time series at a grid point, we exploit the connection
between asymmetric Laplace regression and quantile regression [26]. Our
regression model is additive and includes a smooth component that represents
the multi-year trend in the time series, as well as modulated harmonic
components that estimate seasonal variation. We incorporated spatial
dependence over the grid points into the model through a Potts model [27]. The
Bayesian inference for the resulting hierarchical model is based on Markov
Chain Monte Carlo (MCMC) methods.

The article is structured as follows. In the next section (Section 2), we
provide additional details about the data. Section 3 introduces the proposed
model-based clustering procedure, and Section 4 presents a simulation study to
illustrate the performance and peculiarities of the proposed procedure.
Section 5 collects the results of the application to the Mediterranean Sea
surface temperature data. In Section 6, we analyze the relative strengths and
weaknesses of our proposal.

##  2 Real data

Surrounded by Europe, Africa, and Asia, the Mediterranean Sea is the largest
semi-enclosed European sea. The strait of Sicily separates the Western
Mediterranean and the Eastern Mediterranean, the major basins, while other
sub-basins such as the Alboran, Ligurian, Tyrrhenian, Adriatic, Ionian, and
Aegean Seas fragment the area in specific zones of interest. The Mediterranean
Sea is connected to other three bodies of water: the Atlantic Ocean, the Black
Sea, and the Red Sea through the Strait of Gibraltar, the Dardanelles, and the
Suez Canal, respectively. The presence of different conditions passing from
shallow waters to deep ones permits us to consider the Mediterranean Sea as a
scale model of the world’s oceans in terms of water circulation and deep water
formation [28, 29].

---
Figure 1: Average Sea Surface Temperature (SST) of the Mediterranean Sea over
the period 1982 - 2012.

The Mediterranean Sea was one of the first world seas where temperature
changes could be attributed to climate change [30]. Air temperature is
strongly influenced by atmospheric phenomena, not even connected to global
warming, but the increase of heat waves in the Mediterranean Sea [31] affected
the marine plant and animal world, in particular, the maintenance of
biodiversity [32]. However, the Sea Surface Temperature (SST) is a significant
indicator of climate change, since it can, in turn, impact atmospheric
circulation. Normally the SST is influenced by various factors, including
solar irradiation intensity, water bathymetry, presence of superficial and
deep currents, and the presence of rivers. The recent history of SST in the
Mediterranean Sea shows a slow decrease in the upper layer (150 m) until the
1980s and then a significant increase of about 0.4$\,{}^{o}$C until 2015 [33].

We have considered a dataset from a reanalysis of SST consisting of monthly
mean SST in Celsius degrees from January 1982 to December 2012 for a total of
$372$ time points (see
https://data.marine.copernicus.eu/product/SST_MED_SST_L4_REP_OBSERVATIONS_010_021/description
for a complete data description). Gridded data is available at 1-degree
resolution for each month, with a total of $251$ grid points extending to the
entire Mediterranean Sea. There are no missing values in the data set.

---
Figure 2: Three sites selected on the Mediterranean Sea (A) and seasonal
profile (B, C, and D).

In Figure 1 we show the overall mean SST for all sites. The mean SSTs range
between 16.5$\,{}^{o}$C and 22.3$\,{}^{o}$C with a northwest-southeast
gradient and, in particular, the coldest water masses are recorded near the
Gulf of Lion, while the warmest ones belong to the southern part of the
Levantine Sea, near the Egyptian coastal area. As reported by Figure 2, among
three selected sites (selected following an NW to SE direction) a strong
seasonal pattern is evident as well as the presence of different temperature
extremes. In addition, the highest values were recorded at all the three sites
during the years 2006 to 2012, especially at site 3 located in the Levantine
Sea between Egypt and Cyprus, indicating a general increase in SST.

The clustering of the different Mediterranean areas seems difficult not only
because of the presence of a shaded increasing trend and a strong seasonality
but also because of the asymmetric distribution of the SST, as reported for
site 3 in the Algerian Basin (Figure 3 \- D).

Figure 3: Average SST distribution by month of observation in the Algerian
Basin.

Classifying SST behavior using robust approaches and considering a specific
quantile of interest (e.g., 90th percentile) can help identify areas with
different trends and seasonal patterns.

##  3 Mixture-model clustering

We consider a finite set of sites $\mathcal{S}=\\{1,\ldots,N\\}$ . Each site
is associated with a given spatial localization. At each site
$i\in\mathcal{S}$, we observe a time series of $T$ observations,
$\boldsymbol{y_{i}}=(y_{i,1},\ldots,y_{i,T})^{\prime}$. The goal is to cluster
the $N$ time series into $K$ groups to detect homogenous regions.

We follow a mixture-model approach [34] to clustering, according to which the
cluster membership of the $i$-th time series is represented by a latent random
variable ${c}_{i}\in\\{1,\ldots,K\\}$, where $c_{i}=k$ indicates that the
$i$-th time series belongs to the group $k$.

###  3.1 Membership model

To model the memberships, we suppose that
$\boldsymbol{c}=(c_{1},\ldots,c_{N})^{\prime}$ are independent and identically
distributed multinomial variables with probability mass function for

$\displaystyle f(\boldsymbol{c}|\boldsymbol{\alpha})$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\frac{\operatorname{exp}\left(\sum_{k=1}^{K}\alpha_{k}I(c_{i}=k)\right)}{\exp\left(\sum_{k=1}^{K}\alpha_{k}\right)}=\frac{\operatorname{exp}\left(\sum_{k=1}^{K}\alpha_{k}n_{k}\right)}{\exp\left(N\sum_{k=1}^{K}\alpha_{k}\right)}$
(1) $\displaystyle=$
$\displaystyle\frac{\operatorname{exp}\left(\sum_{k=1}^{K}\alpha_{k}n_{k}\right)}{Z(\boldsymbol{\alpha})}$

where $n_{k}=\sum_{i=1}^{N}I(c_{i}=k)$,
$Z(\boldsymbol{\alpha})=\exp\left(N\sum_{k=1}^{K}\alpha_{k}\right)$, with
$\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\ldots,\alpha_{K})^{\prime}$, and
$\alpha_{1}=0$ for model identifiability, since $N=\sum_{k=1}^{K}n_{k}$.

However, when time series are collected over a spatial domain, cluster
memberships are characterized by some dependence. In this paper, the spatial
dependence among site memberships on a network is modeled using a Markov
Random Field (MRF).

We denote with $i\sim l$ if $i$ and $l$ are neighbors and the set $\partial
i=\\{l:i\sim l\\}$ indicates the neighborhood of $s_{i}$. A MRF is specified
in terms of a neighborhood system, i.e., the set of the neighborhood
$\mathcal{G}=\\{\partial i,\quad i=1,\ldots,n\\}$. The Potts model [27, 35] is
a MRF with probability mass function

$f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})=\frac{\exp\left(\sum_{k=1}^{K}\alpha_{k}n_{k}+\sum_{i\sim
j}\sum_{k\leq
l}\beta_{kl}I(c_{i}=k)I(c_{j}=l)\right)}{Z(\boldsymbol{\alpha},\boldsymbol{\beta})}$

with
$\boldsymbol{\beta}=(\beta_{11},\beta_{12},\ldots,\beta_{1K},\beta_{22},\beta_{23},\ldots,\beta_{2K},\ldots,\beta_{(K-1)K})^{\prime}$.
Here $\sum_{i\sim j}$ means that the summation is taken over all the pairs of
neighboring sites. The parameter $\boldsymbol{\beta}$ controls the spatial
dependence, allowing a spatial dependence inside each cluster and spatial
interaction between two clusters.

Note that, for a moderate or large number of sites, the calculation of the
normalizing constant

$Z(\boldsymbol{\alpha},\boldsymbol{\beta})=\sum_{c_{1}=1}^{K}\cdots\sum_{c_{n}=1}^{K}\exp\left(\sum_{k=1}^{K}\alpha_{k}n_{k}+\sum_{i\sim
j}\sum_{k\leq l}\beta_{kl}I(c_{i}=k)I(c_{j}=l)\right)$

is intractable because it involves $K^{N}$ terms.

In this paper, we examine a simplified variation of the Potts model [22, 36].
We achieve this by setting $\beta_{kl}=0$ for $k\neq l$ and
$\beta_{kk}=\beta_{k}$.

$\displaystyle f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})$
$\displaystyle=$
$\displaystyle\frac{\exp\left\\{\sum_{k=1}^{K}\left(\alpha_{k}n_{k}+\beta_{k}\sum_{i\sim
j}I(c_{i}=k)I(c_{j}=k)\right)\right\\}}{Z(\boldsymbol{\alpha},\boldsymbol{\beta})}$
(2) $\displaystyle=$
$\displaystyle\frac{\exp\left\\{\sum_{k=1}^{K}\left(\alpha_{k}n_{k}+\beta_{k}\sum_{i=1}^{N}n_{k,i}/2\right)\right\\}}{Z(\boldsymbol{\alpha},\boldsymbol{\beta})},$

where $n_{k,i}=\sum_{j\in\partial_{i}}I(c_{i}=k)I(c_{j}=k)$ is the number of
the neighbors of the site $i$ with membership equal to $k$. Consequently, we
explore distinct spatial dependence parameters ($\beta_{k}$) in each cluster.
The parameter $\beta_{k}$ can be interpreted as a regularization parameter
where higher values of $\beta_{k}$ favor clustering to cluster $k$. By setting
$\beta_{k}=0$ for all $k$, model (2) degenerates to the independent model (1).

###  3.2 Time series model

We are mainly interested in clustering the time series according to their
spatio-temporal pattern, so we model the $i$-th time series by
$y_{i,t}=\mu_{c_{i}}(t)+e_{i,t},$ where the deterministic function
$\mu_{c_{i}}(t)$ represents trend and/or seasonality for each group, and
$e_{i,t}$ represents an additive independent error. A flexible specification
of $\mu_{c_{i}}(t)$ is given by a linear combination of known functions, so
that

$\mu_{c_{i}}(t)=\sum_{j=0}^{J}\gamma_{c_{i},j}\psi_{j}(t)=\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}$
(3)

where $\boldsymbol{\psi}(t)=(\psi_{0}(t),\ldots,\psi_{J}(t)^{\prime}$ and
$\boldsymbol{\gamma}_{c_{i}}=(\gamma_{c_{i},0},\ldots,\gamma_{c_{i},J})^{\prime}$.
In this way, we result in a different parametric linear model for each group.
The parameters $\gamma_{c_{i},0}$ can introduce an offset in the data, by
setting the function $\psi_{0}(t)=1$. There are many possible choices for the
functions: polynomial, kernel, sigmoidal, spline, trigonometric, wavelet, etc.
A common assumption in clustering time series is that the additive error
$e_{i,t}$ is symmetric and Gaussian distributed with mean $0$ and standard
deviation $\sigma_{c_{i}}$. This implies that
$\mu_{c_{i}}(t)=\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}$ is
the mean regression function, which is fitted by evaluating a quadratic loss
function, i.e.
$(y_{i,t}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}$, in
the likelihood. However, we could be interested in considering a loss function
with more robust properties or modelling different parameters of the marginal
distribution of the time series.

Here we follow the regression quantile approach [37] and assume that $e_{i,t}$
has the density function $d(\cdot|p)$ and $p$-th, $0<p<1$, quantile equal to
zero, that is, $\int_{-\infty}^{0}d(e|p)de=p$. It turns out that the $p$-th
quantile of $y_{i,t}$ is simply
$Q(y_{i,t})=\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}$. In
particular, we assume that $e_{i,t}$ has Asymmetric Laplace Distribution,
$\operatorname{ALD}(p,\sigma)$, [38] with density

$f(e|\sigma^{2})=\frac{p(1-p)}{\sigma}\exp\left\\{-\rho_{p}\left(\frac{e}{\sigma}\right)\right\\}$
(4)

where $\rho_{p}(\cdot)$ is the loss function $\rho_{p}(u)=u(p-I(u<0))$ and
$\sigma>0$. Therefore, the density of $\boldsymbol{y}_{i}$ given $c_{i}$ is

$\displaystyle
f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma_{c_{i}}^{2})$
$\displaystyle=$
$\displaystyle\prod_{t=1}^{T}f({y}_{it}|\boldsymbol{\gamma}_{c_{i}},\sigma_{c_{i}}^{2})$
(5) $\displaystyle=$
$\displaystyle\prod_{t=1}^{T}\frac{p(1-p)}{\sigma_{c_{i}}}\exp\left\\{-\rho_{p}\left(\frac{y_{i,t}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}}{\sigma_{c_{i}}}\right)\right\\}$
$\displaystyle=$
$\displaystyle\left\\{\frac{p(1-p)}{\sigma_{c_{i}}}\right\\}^{T}\exp\left\\{-\sum_{t=1}^{T}\rho_{p}\left(\frac{y_{i,t}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}}{\sigma_{c_{i}}}\right)\right\\}$

with $\boldsymbol{y}_{i}=(y_{i1},\ldots,y_{iT})^{\prime}$. For example, when
$p=1/2$, the conditional median regression is more robust than the conditional
mean regression with outliers in the observations [39].

##  4 Model estimation and inference

We consider the full set of observations
$\boldsymbol{y}=\\{\boldsymbol{y}_{1}^{\prime},\ldots,\boldsymbol{y}_{N}^{\prime}\\}$
and we want to make inference on the model parameters
$\boldsymbol{\theta}=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}^{\prime},\boldsymbol{\sigma^{2}}^{\prime})^{\prime}$,
with
$\boldsymbol{\gamma}=(\boldsymbol{\gamma}_{1}^{\prime},\ldots,\boldsymbol{\gamma}_{K}^{\prime})^{\prime}$
and $\boldsymbol{\sigma^{2}}=(\sigma^{2}_{1},\ldots,\sigma^{2}_{K})$. The
joint distribution of $\boldsymbol{y},\boldsymbol{c}$ is given by

$f(\boldsymbol{y},\boldsymbol{c}|\boldsymbol{\theta})=\prod_{i=1}^{N}f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}})f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta}),$

and the likelihood is

$f(\boldsymbol{y}|\boldsymbol{\theta})=\sum_{\boldsymbol{c}}\prod_{i=1}^{N}f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}})f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})$

Its evaluation is computationally demanding due to the normalizing constant
$Z(\boldsymbol{\alpha},\boldsymbol{\beta})$. A workaround relies on a pseudo-
likelihood approximation [40] of
$f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})$, namely

$f^{*}(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})=\prod_{i=1}^{N}f(c_{i}|c_{j},j\sim
i,\boldsymbol{\alpha},\boldsymbol{\beta}),$

where the conditional distribution for a location $i$, due to the Markov
property of Potts model, is given by

$f(c_{i}|c_{j},j\sim
i,\boldsymbol{\alpha},\boldsymbol{\beta})=\cfrac{\exp(\alpha_{c_{i}}+\beta_{c_{i}}n_{c_{i},i})}{\sum_{k=1}^{K}\exp(\alpha_{k}+\beta_{k}n_{k,i})}.$
(6)

This approach has been criticized in the literature [41] despite the reduction
in complexity. In the following, we follow a Bayesian approach to estimation
and first specify a prior distribution $\pi(\boldsymbol{\theta})$ for the
parameter $\boldsymbol{\theta}$.

###  4.1 Prior distribution

We assume independent prior distributions for each component of
$\boldsymbol{\theta}$, namely

$\pi(\boldsymbol{\theta})=\pi(\boldsymbol{\alpha})\times\pi(\boldsymbol{\beta})\times\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k})\times\prod_{k=1}^{K}\pi(\sigma^{2}_{k}).$
(7)

and we remember that we have imposed $\alpha_{1}=0$ in the equation (2). The
prior distribution for the remaining parameters $\boldsymbol{\alpha}$ is the
$K-1$ multivariate Gaussian distribution
$\mathcal{N}(\boldsymbol{0},a\boldsymbol{I})$ with an improper prior for $a\gg
0$. The prior for each $\beta_{k}$ is a uniform distribution on $(0,B)$, where
$B$ is the maximum possible value of $\beta_{k}$. Note that our a priori
preference is for the spatial aggregation. The value of $B$ can greatly
influence the spatial aggregation. It should be noted that when the length of
the time series increases, the likelihood term in the a posteriori
distribution formula for the memberships becomes predominant defeating the
role of aggregation. As a consequence, if the model is used to favor spatial
aggregation and pattern simplification, the value of $B$ should be adjusted
according to the length of the time series. In fact, in numerical examples
where the time series have different lengths, we use two values for $B$.
Specifically, $B=10$ for the simulation experiments and $B=100$ for the real
data case.

We assume that the coefficients $\boldsymbol{\gamma}_{k}$, $k=1,\ldots,K$ in
(3) are independent and identically distributed according to the $K$
multivariate Gaussian distribution
$\mathcal{N}(\boldsymbol{0},g\boldsymbol{I})$, where $g\gg 0$. Finally,
similar to Benoit and Van den Poel [42], we propose the inverse gamma
distribution $\mathrm{IG}(s_{0},d_{0})$ as the prior for each $\sigma_{k}$,
where the shape and scale parameters $s_{0}$ and $d_{0}$, respectively, are
known.

###  4.2 Hybrid Gibbs sampler

The evaluation of the posterior distribution
$\pi(\boldsymbol{\theta}|\boldsymbol{y})\propto
f(\boldsymbol{y}|\boldsymbol{\theta})\pi(\boldsymbol{\theta})$ requires an
implementation of a Markov Chain Monte Carlo (MCMC) algorithm [43]. The
popular Gibbs sampler successively draws samples according to the full
conditional distributions associated with the distribution of
$\pi(\boldsymbol{\theta}|\boldsymbol{y})$. When a conditional distribution
cannot be easily sampled, one can resort to a Metropolis-Hastings (MH) move,
which generates samples according to an appropriate proposal and accept or
reject these generated samples with a given probability. The MCMC algorithm is
called hybrid Metropolis-within-Gibbs sampler and generates samples that are
asymptotically distributed according to the posterior distribution.

To develop our MCMC algorithm we exploit a mixture representation [44, 45] of
$e_{i,t}$ in the equation (4) for sampling from the conditional distributions.
This trick is known as completion [43, Definition 10.3]. More specifically,
let $\nu_{it}$ and $w_{it}$ be mutually independent random variables, with
$\nu_{it}$ being a standard Gaussian random variable and $w_{it}$ being a
standard exponential random variable. Given $c_{i}$, $y_{it}$ can be rewritten
as

$\displaystyle
y_{it}=\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}+\tau\sigma_{c_{i}}w_{it}+\omega\sigma_{c_{i}}\sqrt{w_{it}}\,\nu_{it}$
(8)

where $\tau=(1-2p)/\\{p(1-p)\\}$, $\omega^{2}={2}/\\{p(1-p)\\}$.

The joint density of $\boldsymbol{y}_{i}$ given
$\boldsymbol{w}_{i}=(w_{i1},...,w_{iT})^{\prime}$ and $c_{i}$ is

$\displaystyle
f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},\boldsymbol{w}_{i})\propto{(\sigma^{2}_{c_{i}})}^{-{T}/{2}}\prod_{t=1}^{T}w_{it}^{-{1}/{2}}\exp\left\\{-\frac{1}{2}\sum_{t=1}^{T}\frac{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}-\tau\sigma_{c_{i}}w_{it})^{2}}{\sigma^{2}_{c_{i}}\omega^{2}w_{it}}\right\\}$
(9)

Therefore instead of devising a MCMC algorithm for sampling from
$\pi(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}}|\boldsymbol{y})$,
we devise an algorithm for sampling from

$\displaystyle\pi(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w}|\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle\prod_{i=1}^{N}f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},\boldsymbol{w}_{i})f(\boldsymbol{w}_{i})\times
f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})\times$
$\displaystyle\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k})\times\pi(\sigma^{2}_{k})\times\pi(\boldsymbol{\alpha})\times\pi(\boldsymbol{\beta})$

with
$\boldsymbol{w}=(\boldsymbol{w}_{1}^{\prime},\ldots,\boldsymbol{w}_{N}^{\prime})^{\prime}$.
Our MCMC algorithm is described below (Algorithm 1).

Algorithm 1 Hybrid MCMC algorithm

1:initial
$\boldsymbol{\alpha}^{(0)},\boldsymbol{\beta}^{(0)},\boldsymbol{c}^{(0)},\boldsymbol{\gamma}^{(0)},\boldsymbol{\sigma^{2}}^{(0)},\boldsymbol{w}^{(0)}$,
number of iterations $M$

2:$m\leftarrow 1$

3:while $m\leq M$ do

4: Generate $\boldsymbol{\gamma}^{(m)}$ from
$\pi(\boldsymbol{\gamma}|\boldsymbol{\alpha}^{(m-1)},\boldsymbol{\beta}^{(m-1)},\boldsymbol{c}^{(m-1)},\boldsymbol{\sigma^{2}}^{(m-1)},\boldsymbol{w}^{(m-1)},\boldsymbol{y})$

5: Generate $\boldsymbol{\sigma^{2}}^{(m)}$ from
$\pi(\boldsymbol{\sigma^{2}}|\boldsymbol{\alpha}^{(m-1)},\boldsymbol{\beta}^{(m-1)},\boldsymbol{c}^{(m-1)},\boldsymbol{\gamma}^{(m)},\boldsymbol{w}^{(m-1)},\boldsymbol{y})$

6: Generate $\boldsymbol{w}^{(m)}$ from
$\pi(\boldsymbol{w}|\boldsymbol{\alpha}^{(m-1)},\boldsymbol{\beta}^{(m-1)},\boldsymbol{c}^{(m-1)},\boldsymbol{\gamma}^{(m)},\boldsymbol{\sigma^{2}}^{(m)},\boldsymbol{y})$

7: Generate $\boldsymbol{c}^{(m)}={c}_{1}^{(m)},\ldots,{c}_{N}^{(m)}$
according to Algorithm 2

8: Generate $\boldsymbol{\alpha}^{(m)}$ and $\boldsymbol{\beta}^{(m)}$
according to Algorithm 3

9:end while

10:$\boldsymbol{\alpha}^{(M)},\boldsymbol{\beta}^{(M)},\boldsymbol{c}^{(M)},\boldsymbol{\gamma}^{(M)},\boldsymbol{\sigma^{2}}^{(M)},\boldsymbol{w}^{(M)}$

First, we describe how to sample from the conditional distribution
$\pi(\boldsymbol{\gamma}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})$,
$\pi(\boldsymbol{\sigma^{2}}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{w},\boldsymbol{y})$
and
$\pi(\boldsymbol{w}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{y})$
since the conjugate priors of the parameters lead to the availability of
closed forms of these conditional distributions (further details are reported
in Section A of the Supplementary Materials).

#### Conditional distribution of $\boldsymbol{\gamma}$

We have
$\pi(\boldsymbol{\gamma}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})=\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k}|{\sigma}^{2}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$,
where
$\pi(\boldsymbol{\gamma}_{k}|{\sigma}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$
is the density of a multivariate Gaussian distribution with vector mean

$\boldsymbol{m}_{k}=\left(\frac{1}{g}\boldsymbol{I}+\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{\Psi}I(c_{i}=k)\right)^{-1}\left(\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\mathbf{u}_{i}I(c_{i}=k)\right)$

and covariance matrix

$\boldsymbol{S}_{k}=\left(\frac{1}{g}\boldsymbol{I}+\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{\Psi}I(c_{i}=k)\right)^{-1}.$

Here $\boldsymbol{W}_{i}$ is a diagonal matrix with entries
$[\sigma^{2}_{c_{i}}\omega^{2}w_{it}]^{-1}$,
$\boldsymbol{\Psi}=[\boldsymbol{\psi}(1),\ldots,\boldsymbol{\psi}(T)]^{\prime}$
is the matrix of covariates and $\mathbf{u}_{i}=(u_{i1},...,u_{iT})^{\prime}$,
with $u_{it}=y_{it}-\tau\sigma_{c_{i}}w_{it}$.

#### Conditional distribution of $\boldsymbol{\sigma}^{2}$

We note that
$\pi(\boldsymbol{\sigma^{2}}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{w},\boldsymbol{y})=\prod_{k=1}^{K}\pi(\sigma_{k}^{2}|\boldsymbol{\gamma}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$.
It turns out (see the Supplementary material) that
$\pi(\sigma_{k}^{2}|\boldsymbol{\gamma}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$
is the distribution of an inverse Gamma random variable,
$\operatorname{IG}(s_{1},d_{1})$, with shape parameter and scale parameter

$s_{1}=s_{0}+n_{k}T/2,\qquad
d_{1}=d_{0}+\frac{1}{2\omega^{2}}\sum_{i=1}^{N}\left\\{\sum_{t=1}^{T}\frac{\left(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}\right)^{2}}{w_{it}}\right\\}I(c_{i}=k).$

Random samples from the inverse Gamma distribution can be extracted from a
Gamma distribution according to the relationship that if
$v\sim\mathrm{Gamma}(a,1/b)$ then $1/v\sim\mathrm{InvGamma}(a,b)$.

#### Conditional distribution of $\boldsymbol{w}$

We note that
$\pi(\boldsymbol{w}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{y})=\prod_{i=1}^{N}\prod_{t=1}^{T}\pi({w}_{it}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},y_{it})$.
It turns out (see again the Supplementary Material) that the conditional
distribution of $w_{it}$ is the Generalized Inverse Gaussian distribution,
$\operatorname{GIG}(d,a,b),$ with density

$f(w;a,b,d)=\frac{(a/b)^{d/2}}{2K_{d}(\sqrt{ab})}w^{(d-1)}e^{-(aw+b/w)/2},\qquad
x>0,$

where $K_{d}$ is a modified Bessel function of the second kind, $d=1/2$,
$\displaystyle a={(\tau^{2}+2\omega^{2})}/{\omega^{2}}$ and $\displaystyle
b={(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}}/{(\sigma^{2}_{c_{i}}\omega^{2})}$.

We can use this trick to sample from this posterior distribution [46]. The
density of the Inverse Gaussian distribution,
$\operatorname{InvGIG}(\lambda,\mu)$,

$\displaystyle f(w;\lambda,\mu)$ $\displaystyle\propto$ $\displaystyle
w^{-{3}/{2}}\exp\left\\{\frac{-\lambda(w-\mu)^{2}}{2\mu^{2}w}\right\\}$
$\displaystyle\propto$ $\displaystyle
w^{-{3}/{2}}\exp\left\\{-\frac{1}{2}\left(\frac{\lambda}{\mu^{2}}w+\frac{\lambda}{w}\right)\right\\}$

is equal to the $\operatorname{GIG}(-1/2,{\lambda}/{\mu^{2}},\lambda)$. We
also know that if $v\sim\operatorname{GIG}(-1/2,b,a)$, then
$1/v\sim\operatorname{GIG}(1/2,a,b)$, i.e. if
$v\sim\operatorname{InvG}({\lambda},\mu)$ then
$1/v\sim\operatorname{GIG}(1/2,\lambda,{\lambda}/{\mu^{2}})$.

Therefore, we can sample from the conditional posterior of $w_{it}$ by drawing
$v_{it}$ from $\operatorname{InvG}({\lambda}/{\mu^{2}},\lambda)$ with
$\displaystyle\lambda=\frac{\tau^{2}+2\omega^{2}}{\omega^{2}}$ and
$\displaystyle\frac{\lambda}{\mu^{2}}=\frac{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}}{\sigma^{2}_{c_{i}}\omega^{2}}$
that implies
$\displaystyle\mu=\frac{\sigma_{c_{i}}\sqrt{\tau^{2}+2\omega^{2}}}{|y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}|}$
and considering $w_{it}=1/v_{it}$.

#### Conditional distribution of $\boldsymbol{c}$

Generating $\boldsymbol{c}^{(m)}$ is difficult because of the normalizing
constant $Z(\boldsymbol{\alpha},\boldsymbol{\beta})$. However, it is possible
to embed another Gibbs step (see Algorithm 2) into Algorithm 1 and construct
an ergodic Markov chain whose stationary distribution is
$\pi(\boldsymbol{\theta}|\boldsymbol{y})$.

Algorithm 2 Gibbs step

initial $\boldsymbol{c}^{m-1}$

1:
$c_{1}^{m}\sim\pi(c_{1}|c_{2}^{m-1},\ldots,c_{N}^{m-1},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})$

2:
$c_{2}^{m}\sim\pi(c_{2}|c_{1}^{m},c_{3}^{m-1}\ldots,c_{N}^{m-1},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})$

$\vdots$

N:
$c_{N}^{m}\sim\pi(c_{N}|c_{1}^{m},\ldots,c_{N-1}^{m},\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})$

$\boldsymbol{c}^{m}$

The conditional distribution $\pi(c_{i}|c_{j},j\neq
i,\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})$
is given by

$\displaystyle\pi(c_{i}|c_{j},j\neq
i,\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})=\cfrac{\exp\left\\{-\frac{1}{2}\sum_{t=1}^{T}\frac{(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}}{w_{it}\sigma^{2}_{c_{i}}\omega^{2}}+\alpha_{c_{i}}+\beta_{c_{i}}n_{c_{i},i}\right\\}}{\sum_{k=1}^{K}\exp\left\\{-\frac{1}{2}\sum_{t=1}^{T}\frac{(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{k})^{2}}{w_{it}\sigma^{2}_{k}\omega^{2}}+\alpha_{k}+\beta_{k}n_{k,i}\right\\}},$

exploiting equation (6).

#### Conditional distribution of $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$

The last step, i.e. the simulation from the conditional distributions of
$\boldsymbol{\alpha}^{(m)}$ and $\boldsymbol{\beta}^{(m)}$ is performed using
a Metropolis-Hastings algorithm. However, the direct implementation of the
algorithm is still difficult due to the presence of the normalizing constant
$Z(\boldsymbol{\alpha},\boldsymbol{\beta})$.

Marjoram et al. [47] have shown that it is possible to define a valid
likelihood-free MH algorithm for posterior distributions with intractable
normalizing constants by introducing a carefully selected auxiliary variable
and a tractable sufficient statistic on the target density.

Since exact likelihood-free MH algorithms have several shortcomings described
in Marin et al. [48], we modify the algorithm described in Pereyra et al. [49]
that leads to an Approximate Bayesian Computation (ABC) with a probability-
free Metropolis-Hastings step. The considered algorithm is described in the
next box (Algorithm 3).

Algorithm 3 ABC likelihood-free MH step for $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$.

1:initial $\boldsymbol{\alpha}^{m-1},\boldsymbol{\beta}^{m-1}$

2:Fix $\alpha_{1}^{*}=0$

3:Generate proposals
$\alpha^{*}_{k}\sim\mathcal{N}(\alpha_{k}^{m-1},\sigma^{2}_{\alpha})$ for
$k=2,\ldots,K$

4:Generate
$\beta_{k}^{*}\sim\mathcal{NT}_{0,B}(\beta_{k}^{m-1},\sigma^{2}_{\beta})$ for
$k=1,\ldots,K$

5:Generate an auxiliary $\boldsymbol{c}^{*}$ from the distribution
$f(\boldsymbol{c}|\boldsymbol{\alpha}^{*},\boldsymbol{\beta}^{*})$ in (2)

6:if
$\rho(\boldsymbol{\eta}(\boldsymbol{c}^{*}),\boldsymbol{\eta}(\boldsymbol{c}^{m-1}))<\mbox{tolerance}$
then

7: Set $\displaystyle
r=\frac{\pi(\boldsymbol{\alpha}^{*})\pi(\boldsymbol{\beta^{*}})\,\prod_{k=2}^{K}\phi(\alpha_{k}^{m-1}|\alpha_{k}^{*},\sigma^{2}_{\alpha})\prod_{k=1}^{K}\phi_{[0,5]}(\beta^{m-1}_{k}|\beta^{*}_{k},\sigma^{2}_{\beta})}{\pi(\boldsymbol{\alpha}^{m-1})\pi(\boldsymbol{\beta}^{m-1})\,\prod_{k=2}^{K}\phi(\alpha_{k}^{*}|\alpha_{k}^{m-1},\sigma^{2}_{\alpha})\prod_{k=1}^{K}\phi_{[0,5]}(\beta^{*}_{k}|\beta^{m-1}_{k},\sigma^{2}_{\beta})}$

8: Draw $u\sim\mathcal{U}(0,1)$

9: if $r<u$ then

10:
$\boldsymbol{\alpha}^{m}=\boldsymbol{\alpha}^{*},\boldsymbol{\beta}^{m}=\boldsymbol{\beta^{*}}$

11: else

12:
$\boldsymbol{\alpha}^{m}=\boldsymbol{\alpha}^{m-1},\boldsymbol{\beta}^{m}=\boldsymbol{\beta}^{m-1}$

13: end if

14:else

15:
$\boldsymbol{\alpha}^{m}=\boldsymbol{\alpha}^{m-1},\boldsymbol{\beta}^{m}=\boldsymbol{\beta}^{m-1}$

16:end if

17:$\boldsymbol{\alpha}^{m},\boldsymbol{\beta}^{m}$

Here $\mathcal{NT}_{a,b}(\mu,\sigma^{2})$ is the truncated Gaussian
distribution with mean $\mu$ and variance $\sigma^{2}$, within the interval
$(a,b)$, with $-\infty\leq a<b\leq+\infty$. Here, the vector
$\boldsymbol{\eta}(\cdot)$ serves as a vector of statistics, while the
function $\rho(\cdot)$ functions as a distance metric.

The entire MCMC algorithm necessitates the specification of initial values in
Algorithm 1 and variances for the proposals, the vector
$\boldsymbol{\eta}(\cdot)$, the distance $\rho(\cdot)$ and the tolerance in
Algorithm 3.

In all numerical experiments, the initial values $\boldsymbol{w}^{(0)}$ were
randomly assigned. The algorithm outlined in Sottile and Adelfio [50] yields
the memberships $\boldsymbol{c}^{(0)}$, while the values of
$\boldsymbol{\alpha}^{(0)},\boldsymbol{\gamma}^{(0)},\boldsymbol{\sigma^{2}}^{(0)}$
are estimated by a multinomial regression model ($\boldsymbol{\alpha}^{(0)}$)
and K different quantile regressions
($\boldsymbol{\gamma}^{(0)},\boldsymbol{\sigma^{2}}^{(0)}$), respectively. The
initial value for the spatial dependence parameter differs in the simulations
and in the real data application. Specifically, we set ${\beta}_{k}^{(0)}$
equal to 2 for $k=1,\ldots,K$ in the simulations and 10 in the real data
application. This choice was arbitrary and motivated solely by the sensitivity
of the sampler’s mixing properties to the starting values, resulting in
prolonged burn-in periods.

We set the variance of proposals $\sigma^{2}_{\alpha}=0.1$ for
$\boldsymbol{\alpha}$. The variance of proposals $\sigma^{2}_{\beta}$ depends
on the choice of the uniform $\mathcal{U}(0,B)$ as prior for
$\boldsymbol{\beta}$. Since we set $B=10$ for the simulation experiments and
$B=100$ for the real data example, we specify $\sigma^{2}_{\beta}$ as 0.04 and
1, respectively. We select
$\boldsymbol{\eta}(\boldsymbol{c})=(\eta_{1}(\boldsymbol{c}),\ldots,\eta_{2K-1}(\boldsymbol{c}))^{\prime}=\left(n_{2},\ldots,n_{K},\sum_{i=1}^{N}n_{1,i},\ldots,\sum_{i=1}^{N}n_{K,i}\right)^{\prime},$
since the vector represents sufficient statistics in the Potts model (2). The
function
$\rho(\boldsymbol{\eta}(\boldsymbol{c}^{*}),\boldsymbol{\eta}(\boldsymbol{c}^{m-1}))=\frac{1}{2K-1}\sum_{j=1}^{2K-1}\frac{\left|{\eta}_{j}(\boldsymbol{c}^{*})-{\eta}_{j}(\boldsymbol{c}^{m-1})\right|}{\eta_{j}(\boldsymbol{c}^{m-1})},$
measures the distance between $\boldsymbol{\eta}(\boldsymbol{c}^{*})$ and
$\boldsymbol{\eta}(\boldsymbol{c}^{m-1})$.

The tolerance is adjusted to achieve a 5% acceptance rate. This value is
calculated using the 5th percentile of the observed tolerances
$\rho(\boldsymbol{\eta}(\boldsymbol{c}^{*}),\boldsymbol{\eta}(\boldsymbol{c}^{m-1}))$
during the burning phase.

##  5 Simulations

We conducted simulations to illustrate how our spatial clustering method can
distinguish three different groups of time series, taking into account both
the temporal pattern and the marginal error distribution. Two membership
allocations (concentric circles and flaglike rectangles) of $N=900$ sites on a
$30\times 30$ regular grid will be tested. We consider three groups for each
allocation, as in Figure 4. Concentric circles have different group sizes,
namely $440$, $244$, and $216$. The flaglike rectangle is divided into groups
of the same size.

Figure 4: Spatial pattern of concentric circles (A) and flaglike rectangles
(B).

###  5.1 First simulation example

The first simulation example attempts to measure the accuracy of clustering
and parameter estimation when in the joint distribution of a time series is
(5) a spatial pattern for membership is included. For each site in the group
$k$ (with $k=1,2,3$) a time series of $T=100$ observations is simulated based
on the model

$y_{t,k}=\mu_{k}(t)+e_{t,k},\qquad t=1,\ldots,T.$ (10)

The function

$\mu_{k}(t)=\frac{\theta_{0k}}{100}t+\theta_{1k}\cos\left(\frac{3\pi}{100}t\right)+\theta_{2k}\sin\left(\frac{3\pi}{100}t\right)$
(11)

represents an overall trend and a cyclic pattern. We set the parameters for
three clusters as follows:
$(\theta_{01},\theta_{11},\theta_{21})=(1,0.5,0.25)$,
$(\theta_{02},\theta_{12},\theta_{22})=(1.25,0.25,0.5)$ and
$(\theta_{03},\theta_{13},\theta_{23})=(1.5,-0.25,0.5)$, and it is assumed
that $e_{t,k}\sim\operatorname{ALD}(p,\sigma)$ with $p=0.5$ and $\sigma=0.5$.
It is worth mentioning that the presence of random ALD errors complicates the
time series clustering situation, as shown in Figure 5. The simulation
experiment is repeated $100$ times.

Figure 5: Simulated time series with ALD random error. The mean trend for each
cluster is reported with a solid line.

We evaluated the effectiveness of two clustering techniques: Asymmetric
Laplace Model-based Clustering with Spatial Dependence (ALC-S) and the
spatially-independent procedure (ALC), where $\beta_{k}=0$ for $k=1,\ldots,K$
in equation (2).

In the estimation stage the vector $\mathbf{\psi}(t)$ in (3) is
$(t/100,\cos(3\pi\,t/100),\sin(3\pi\,t/100))^{\prime}$ to facilitate a
comparison between the estimated and actual parameters. We use the MCMC
algorithm for 1000 iterations to obtain posterior estimates, discarding the
first $300$ as burn-in. Based on trace plots, this number of iterations
appears to be adequate for Markov chain convergence. The group membership is
determined by the posterior mode of $c_{i}$, while the regression parameter
estimates are calculated using the posterior mean. For the ALC-S method, three
nearest neighbor networks have been considered: specifically, the 4-NN, 8-NN,
and 12-NN.

The performance of both clustering algorithms in reconstructing the three
clusters is evaluated by gauging the degree of agreement between the estimated
partition and the actual partition, as determined by the Adjusted Rand Index
(ARI) [51]. The results are outlined in Table 1. The results suggest that the
inclusion of the spatial layer improves the consistency between the clustering
and the true spatial arrangement. On average, this enhancement leads to a
13-14% improvement in the ARI values for both the 4-NN spatial configurations
and results in values greater than 90%. The concordance in the concentric
circle scenario decreases as the number of neighbors increases. As reported in
Table 2, both methods show exceptional performance in estimating the actual
curve parameters. In particular, ALC-S exhibits higher accuracy, although it
may be slightly less accurate in certain situations such as the concentric
circle scenario.

Table 1:  First simulation example. Average ARI values of $100$ replications for each experiment and clustering method (standard error between parentheses). | Concentric circles | Flaglike rectangles
---|---|---
ALC-S with 4-NN | $0.914\,(0.062)$ | $0.923\,(0.075)$
ALC-S with 8-NN | $0.861\,(0.157)$ | $0.934\,(0.122)$
ALC-S with 12-NN | $0.801\,(0.194)$ | $0.941\,(0.145)$
ALC | $0.794\,(0.043)$ | $0.762\,(0.028)$
Table 2:  First simulation example. Average estimates for $\theta$ were obtained from 100 replications for each clustering method. Standard errors are reported in parentheses. Concentric circles |
---|---
| Cluster 1
| $\theta_{10}=1$ | $\theta_{11}=0.5$ | $\theta_{12}=0.25$
ALC-S with 4-NN | 1.001 (0.012) | 0.504 (0.011) | 0.244 (0.015)
ALC-S with 8-NN | 1.002 (0.016) | 0.500 (0.012) | 0.248 (0.018)
ALC-S with 12-NN | 1.026 (0.035) | 0.475 (0.046) | 0.271 (0.035)
ALC | 1.001 (0.010) | 0.508 (0.010) | 0.240 (0.013)
| Cluster 2
| $\theta_{20}=1.25$ | $\theta_{21}=0.25$ | $\theta_{22}=0.5$
ALC-S with 4-NN | 1.252 (0.012) | 0.239 (0.023) | 0.510 (0.017)
ALC-S with 8-NN | 1.253 (0.024) | 0.240 (0.029) | 0.510 (0.025)
ALC-S with 12-NN | 1.256 (0.066) | 0.224 (0.097) | 0.509 (0.070)
ALC | 1.244 (0.021) | 0.237 (0.021) | 0.504 (0.015)
| Cluster 3
| $\theta_{30}=1.5$ | $\theta_{31}=-0.25$ | $\theta_{32}=0.5$
ALC-S with 4-NN | 1.502 (0.014) | -0.256 (0.017) | 0.497 (0.010)
ALC-S with 8-NN | 1.500 (0.011) | -0.253 (0.016) | 0.500 (0.013)
ALC-S with 12-NN | 1.481 (0.047) | -0.212 (0.078) | 0.497 (0.023)
ALC | 1.500 (0.012) | -0.256 (0.015) | 0.500 (0.011)
|  |  |
Flaglike rectangles |
| Cluster 1
| $\theta_{10}=1$ | $\theta_{11}=0.5$ | $\theta_{12}=0.25$
ALC-S with 4-NN | 0.995 (0.017) | 0.511 (0.020) | 0.237 (0.021)
ALC-S with 8-NN | 0.999 (0.013) | 0.503 (0.009) | 0.246 (0.009)
ALC-S with 12-NN | 1.015 (0.043) | 0.487 (0.043) | 0.264 (0.042)
ALC | 1.005 (0.018) | 0.507 (0.018) | 0.242 (0.017)
| Cluster 2
| $\theta_{20}=1.25$ | $\theta_{21}=0.25$ | $\theta_{22}=0.5$
ALC-S with 4-NN | 1.247 (0.021) | 0.248 (0.025) | 0.501 (0.021)
ALC-S with 8-NN | 1.250 (0.011) | 0.249 (0.009) | 0.503 (0.009)
ALC-S with 12-NN | 1.251 (0.038) | 0.244 (0.031) | 0.510 (0.057)
ALC | 1.251 (0.022) | 0.237 (0.028) | 0.504 (0.018)
| Cluster 3
| $\theta_{30}=1.5$ | $\theta_{31}=-0.25$ | $\theta_{32}=-0.5$
ALC-S with 4-NN | 1.501 (0.012) | -0.257 (0.013) | 0.501 (0.011)
ALC-S with 8-NN | 1.498 (0.010) | -0.253 (0.009) | 0.501 (0.008)
ALC-S with 12-NN | 1.502 (0.011) | -0.252 (0.009) | 0.501 (0.011)
ALC | 1.496 (0.012) | -0.255 (0.019) | 0.500 (0.009)

###  5.2 Second simulation example

The second, more extensive example compares the proposed procedure with nice
competing procedures that have been implemented in R. The goal is still to
provide a comprehensive analysis and evaluation of the effectiveness of the
proposed procedure relative to the alternatives. We keep the two spatial
patterns previously analyzed, the same number of time series observations, and
the same functions (11) for the clusters.

For the distribution of the error term in (10), we defined two different
scenarios (A and B). In scenario A, it is assumed that $e_{t,k}$ follows a
Gaussian distribution with mean 0 and standard deviation 1. This scenario is
ideal for model-based Gaussian clustering to achieve high efficacy.

In scenario B, for each cluster we define the following error distribution for

  1. 1.

$e_{t,1}=N_{t}$;

  2. 2.

$e_{t,2}=U_{t}N_{t}+(1-U_{t})(G_{t}-1)$;

  3. 3.

$e_{t,3}=G_{t}-1$,

with $U_{t}\sim Ber(0.5)$, $N_{t}\sim\mathcal{N}(0,1)$ and
$G_{t}\sim\mbox{Gamma}(1,1)$ independent random variables. The three error
distributions share the same mean of $0$ and standard deviation of $1$.
However, the skewness of the error distribution increases with the cluster
number.  An example of a simulated series for both scenarios is reported in
Figure 6.

---
Figure 6: Simulated time series under scenario A (top) and B (bottom).

We estimate the function $\mu_{k}(t)$ for each cluster using a $J$ cubic
B-spline basis in (3), with equally spaced knots over $[0,T]$. The preliminary
selection of $J$ is done using a BIC-like overall criterion as follows.

  1. 1.

We estimate the parameter ${\boldsymbol{\gamma}}_{i}$ of model (5) for each
time series ${\boldsymbol{y}}_{i}$ by minimizing
$\sum_{t=1}^{T}\rho_{p}(y_{i,t}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{i})$.

  2. 2.

We evaluate the BIC-like criterion
$BIC_{i}(J)=\sum_{t=1}^{T}\rho_{p}(y_{it}-\boldsymbol{\psi}(t)^{\prime}\widehat{\boldsymbol{\gamma}}_{i})+J\log(T)$.

  3. 3.

We find the value $J$ that minimizes the overall value
$\overline{BIC}(J)=\sum_{i=1}^{N}BIC_{i}(J)$.

The minimum value was achieved at $J=6$ in the experiment. We continue to
utilize the MCMC algorithm for estimation, completing $1000$ iterations (with
$300$ as burn-in) for the ALC-S procedure utilizing an 8-NN network. Nine
clustering methods were assessed as competitors.

  1. 1)

Hierarchical Clustering based on the raw data (HC-R). We use the DIvisive
ANAlysis Clustering (DIANA) algorithm [52], with the dissimilarity measure
between two time series $\boldsymbol{y}_{i}$ and $\boldsymbol{y}_{i^{\prime}}$
$d_{ij}^{R}=\|\boldsymbol{y}_{i}-\boldsymbol{y}_{j}\|$ and we extract three
clusters.

  2. 2)

Hierarchical Clustering based on the Least Squares regression estimates (HC-
LS). For each time series, we get the least squares estimates of
${\boldsymbol{\gamma}}_{i}$ in (3), $\widehat{\boldsymbol{\gamma}}^{LS}_{i}$,
and we use the DIANA algorithm with the dissimilarity measure
$d_{ij}^{LS}=\|\boldsymbol{\gamma}^{LS}_{i}-\boldsymbol{\gamma}^{LS}_{j}\|$

  3. 3)

Hierarchical Clustering based on the Least Absolute Deviation regression
estimates (HC-AD): we get the absolute deviation regression estimates of
${\boldsymbol{\gamma}}_{i}$ in (3), $\widehat{\boldsymbol{\gamma}}^{AD}_{i}$
and we still use the DIANA algorithm with
$d_{ij}^{AD}=\|\boldsymbol{\gamma}^{AD}_{i}-\boldsymbol{\gamma}^{AD}_{j}\|$.

  4. 4)

Geostatistical Spectral Clustering (GSC) which is based on the spectral
decomposition of the graph Laplacian [53];

  5. 5)

Functional Principal Components Analysis Clustering (FPCAC) based on the
principal component analysis for functional data obtained by quantile
regression [50] as implemented in the R package clustEff;

  6. 6)

Functional Data Clustering Using Adaptive Density Peak Detection (FADP) which
the density is estimated by functional k-nearest neighbor density estimation
based on a proposed semi-metric between functions [54] as implemented in the R
package FADPclust;

  7. 7)

The FunEM algorithm is based on a discriminative functional mixture model that
enables data clustering in a unique and discriminative functional subspace
[55]. The procedure is implemented using the R package funFEM;

  8. 8)

Functional data clustering method for High-Dimensional Data Clustering
(FunHDDC), based on a functional latent mixture model [56]. The algorithm is
implemented in the R package funHDDC;

  9. 9)

Asymmetric Laplace model-based Clustering without spatial dependence (ALC).

To address possible spatial dependence, we modify the measure of dissimilarity
between competitors 1-3. Specifically, we define $d^{w}_{ij}$ as $d_{ij}$
multiplied by $\gamma(h_{ij})$, where $\gamma(h_{i,j})$ is the variogram
computed for the distance $h_{i,j}$ between site $i$ and site $j$. We consider
an exponential model with a nugget effect for the variogram. We fit the
variogram using regression scores derived from the first principal component
extracted from the raw data [57].

Table 3: Average ARI values of 100 replications for each experiment and clustering method (standard error between parentheses). | Concentric circles | Flaglike rectangles
---|---|---
| Scenario A | Scenario B | Scenario A | Scenario B
HC-R | $0.467\,(0.111)$ | $0.517\,(0.103)$ | $0.736\,(0.176)$ | $0.488\,(0.148)$
HC-LS | $0.392\,(0.074)$ | $0.392\,(0.076)$ | $0.421\,(0.040)$ | $0.418\,(0.045)$
HC-AD | $0.318\,(0.074)$ | $0.261\,(0.056)$ | $0.322\,(0.032)$ | $0.335\,(0.037)$
GSC | $0.348\,(0.067)$ | $0.537\,(0.351)$ | $0.459\,(0.181)$ | $0.559\,(0.417)$
FPCAC | $0.828\,(0.020)$ | $0.826\,(0.026)$ | $0.838\,(0.021)$ | $0.844\,(0.024)$
FADP | $0.730\,(0.062)$ | $0.667\,(0.131)$ | $0.728\,(0.065)$ | $0.728\,(0.072)$
FunFEM | $0.169\,(0.020)$ | $0.172\,(0.017)$ | $0.179\,(0.024)$ | $0.168\,(0.019)$
FunHDDC | $0.370\,(0.082)$ | $0.418\,(0.062)$ | $0.421\,(0.025)$ | $0.412\,(0.052)$
ALC | $0.752\,(0.038)$ | $0.703\,(0.023)$ | $0.717\,(0.038)$ | $0.685\,(0.032)$
ALC-S | $0.815\,(0.066)$ | $0.866\,(0.135)$ | $0.808\,(0.077)$ | $0.968\,(0.043)$

We evaluated the performance of clustering utilizing the adjusted Rand index
(ARI) [51]. Table 3 displays the adjusted Rand Index (ARI) values for the ten
clustering approaches in each scenario. From these findings, it is evident
that our clustering algorithm demonstrated an overall satisfactory performance
(ARI $>0.80$). In contrast, distance matrix-based methods (competitors 1-3)
showed low levels of agreement, especially in the presence of a complex
spatial pattern (concentric circles). A similar behavior was observed in the
GSC procedure. Non-spatial methods based on robust quantile regression (FPCAC
and FADP) worked moderately well, while the two considered procedures for
functional data clustering (FunFEM and FunHDDC) reported poor classification
performances. Taking into account scenario B, ALC-S reported the highest ARI
values (0.866 and 0.968, for concentric circles and flag-like rectangles,
respectively). Our proposal exhibited an increase in the average
classification accuracy in each scenario upon the addition of spatial
information, as previously stated in the previous simulation subsection.

##  6 Clustering of the Mediterranean Sea SST

In this section, we perform a spatial clustering of the Mediterranean SST
monthly time series, considering a temporal window between 1982 and 2012; the
reference data are previously presented in Section 2. As shown in Figure 2, we
can observe a strong and time-varying seasonal pattern, i.e. non-stationarity
in time, with the presence of a maximum in late summer and a minimum in the
spring months. Another distinctive element of this dataset is the dependence
of the mean sea surface temperature on latitude. Higher average values are
observed at sites close to the African continent, while near the Gulf of Lion
we can find the coldest water masses.  To eliminate the effect of
latitude/longitude, we convert each time series to its normalized anomalies
$(y_{it}-\bar{y}_{i})/s_{i}$, $t=1,\ldots,T$. Here, $\bar{y}_{i}$ and $s_{i}$
represent the sample mean and standard deviation, respectively, of the time
series $i$. After the transformation, the Potts model specified in equation
(2) seems appropriate for modeling spatial membership clustering. The network
specification based on sites and their neighbors allows the consideration of
non-convex domains. To account for edge effects, the number of neighbors of a
site at the edge has been modified. Apart from edge effects, the spatial
dependence between clusters is assumed to be different for each group and
isotropic. This fact may be a limitation, but a parameterization that would
allow for anisotropy would require a large number of parameters.

We specify a seasonal modulation model [58], namely

$\mu(t)=m(t)+s(t)$

where the component $m(t)$ defines an overall trend while the seasonality
$s(t)$ can be expressed as the sum of $D$ harmonics

$s(t)=\sum_{d=1}^{D}\left\\{g_{2d}(t)\cos(\pi dt/6)+g_{3d}(t)\sin(\pi
dt/6)\right\\}.$

In our example, we consider $D=2$, to capture intra-annual variations in the
seasonal component. The model can be structured by a linear combination of
$J_{1}$, $2\times J_{2}$, and $2\times J_{3}$ cubic B-splines, with equally
spaced knots over the time interval $[1,372]$. It is easy to see that the
resulting $\mu(t)$, namely

$\displaystyle\mu(t)$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{J_{1}}\gamma_{1l}\psi_{1j}(t)+\sum_{d=1}^{2}\left\\{\left[\sum_{j=1}^{J_{2}}\gamma_{2dj}\psi_{2dj}(t)\right]\cos(\pi
dt/6)+\left[\sum_{j=1}^{J_{3}}\gamma_{3dj}\psi_{3dj}(t)\right]\sin(\pi
dt/6)\right\\}$ (12) ,

can be written as a linear combination of known functions.

A preliminary data analysis was performed to determine the number of functions
in (12) and the number of clusters. Following the strategy outlined in Section
5, we minimize the total $\operatorname{BIC}(J_{1},J_{2},J_{3})$ using a
median (i.e. $p=0.5$) regression, varying the value of $J_{1}$, $J_{2}$, and
$J_{3}$ from 4 to 8.

The total value $\operatorname{BIC}(J_{1},J_{2},J_{3})$ is reported in the
Supplementary Materials (Table 1), where the minimum value is reported for
$J_{1}=4$ for the trend (including an internal intercept), $J_{2}=8$, and
$J_{3}=4$. As a neighborhood network for Potts model 2, we have considered the
4-NN and 8-NN network.

Then, we identify the number of clusters that minimize the Widely Applicable
Information Criterion ($\operatorname{WAIC}$) [59], computed as follows:

$\operatorname{WAIC}=-2\left(\sum_{i=1}^{N}\log\mathbb{E}_{\boldsymbol{\theta},c_{i}}[f(\boldsymbol{y}_{i}|\boldsymbol{\theta},c_{i})|\boldsymbol{y}]-\sum_{i=1}^{N}\mathbb{V}_{\boldsymbol{\theta},c_{i}}[\log
f(\boldsymbol{y}_{i}|\boldsymbol{\theta},c_{i})|\boldsymbol{y}]\right).$

The expectation
$\mathbb{E}_{\boldsymbol{\theta},c_{i}}[f(\boldsymbol{y}_{i}|\boldsymbol{\theta},c_{i})|\boldsymbol{y}]$
and the variance $\mathbb{V}_{\boldsymbol{\theta},c_{i}}[\log
f(\boldsymbol{y}_{i}|\boldsymbol{\theta},c_{i})|\boldsymbol{y}]$ are estimated
by the sample mean of
$f(\boldsymbol{y}_{i}|\boldsymbol{\theta}^{(m)},c_{i}^{(m)})$, and the sample
variance of $\log
f(\boldsymbol{y}_{i}|\boldsymbol{\theta}^{(m)},c_{i}^{(m)})$, $m=1\ldots,M$ of
the posterior simulations derived from running the MCMC algorithm after
convergence.

The WAIC index is calculated for a range of clusters, specifically from 2 to
8. As shown in the Supplementary Materials, we found this minimum to be equal
to 8 and 6 clusters for the 4-NN and 8-NN networks, respectively. The final
membership ${c}_{i}$ and coefficients $\boldsymbol{\gamma_{k}}$ were estimated
using the mode and mean, respectively, of the posterior distribution after
3000 MCMC iterations, with the first 500 iterations discarded as burn-in.

We present the clustering results obtained from grouping data into $K=6$
clusters with the 8-NN network and compare them to the results obtained using
the ALC method, disregarding spatial dependence.

---
Figure 7: Spatial clustering by the ALC and ALC-S for K=6 clusters.

Cluster 1 characterizes a small area that overlaps with the Gulf of Genoa,
whereas Cluster 2 in dark yellow encompasses the shallow waters of the
Adriatic Sea, the western Aegean Sea, and the Nile’s mouths. Cluster 3 extends
to the central-western part of the Mediterranean Sea, while Cluster 4 in light
blue defines another small area that coincides with the Strait of Gibraltar.
Cluster 5, depicted in dark blue, encompasses the portion of the Mediterranean
Sea that lies between the coasts of Turkey and those of Greece, as well as
Cyprus and the island of Crete. Cluster 6 denotes a distinct area within the
central-southern Mediterranean Sea, extending from Crete to the Egyptian
coasts. Compared to the ALC procedure, the ALC-S provides a more consistent
categorization in terms of the spatial extent and characterization of the
water masses (Figure 7). In particular, ALC-S can delineate two small regions
(Gulf of Genoa and Strait of Gibraltar) with distinct temperature patterns.

We further analyzed the classification results of several clustering methods
that exhibited excellent performance in the previous simulations section
(Section 5). We selected the hierarchical clustering method based on raw data,
the GSC algorithm, and the FPCAC procedure using $93$ B-spline bases (3 bases
for each year). Although it is not a spatial method, it obtained high ARI
values in the simulations. Since these non-model-based methods were used, we
selected the number of clusters that maximized the average silhouette width by
utilizing the fviz_nbclust function, which is available in the R package
factoextra. Please refer to Figure 1 in the Supplementary Materials for
additional information.

Figure 8 shows the spatial pattern of the three clustering methods considered.
The partitions obtained by the HC-R procedures suggest 3 clusters and largely
separate the Mediterranean basin into three parts, resulting in a coarse
classification (Figure 1A). The results obtained by GSC and FPCAC suggest
dividing the Mediterranean Sea into two groups based on their approaches:
eastern vs. western for GSC and southern vs. northern for FPCAC. However, the
GSC results appear to be more spatially consistent.

---
Figure 8: Spatial clustering by the three alternative clustering approaches.  Table 4: Summary statistics of the SST (in $\,{}^{o}$C) for each of the clusters identified by ALC-S procedure based on 8-NN network. Cluster | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
Size | 2 | 21 | 107 | 2 | 24 | 95
Average | 18.2 | 18.3 | 19.0 | 17.9 | 20.2 | 20.8
Standard Deviation | 4.2 | 4.7 | 4.2 | 2.4 | 4.0 | 4.0
Median | 17.7 | 17.7 | 18.3 | 17.4 | 19.7 | 20.4
Median Abs. Deviation | 5.6 | 6.0 | 5.4 | 2.7 | 4.9 | 5.4
Min | 12.4 | 8.2 | 11.4 | 14.2 | 12.5 | 12.9
Max | 27.5 | 30.8 | 28.8 | 25.6 | 29.9 | 29.8
Average 10-year increase °C | 0.32 | 0.53 | 0.56 | 0.43 | 0.48 | 0.61
Median intra-annual $\operatorname{IQR}$ | 3.76 | 8.69 | 7.64 | 7.79 | 7.41 | 7.00
Figure 9: Estimated median values for the standardized temperature anomalies:
trend (A) and seasonal (B) components for each cluster resulted by ALC-S
procedure with 8-NN.

To better understand the differences between the defined groups of time
series, the position and variability indices, including the median intra-
annual IQR, are reported in Table 4. For instance, Cluster $1$ and $4$ have a
similar average SST value, but a marked difference in terms of range and
intra-annual variation. Clusters $5$ and $6$ are spatially close in terms of
average temperature, but show different behavior in terms of average
temperature increase (0.48 vs. 0.61 °C decade-1). The clusters $2$ and $3$ are
very similar in terms of median/mean, but report a different seasonal pattern.
It is important to note that our classification is similar to that obtained by
other studies considering a wide range of physical indicators [60] or looking
at the trophic status based on chlorophyll type-a concentration [61], showing
how changes in temperature are related to many other environmental indicators.

Figure 10: Amplitude and phase function of each harmonic estimated by ALC-S
procedure with 8-NN.

Looking at the median trend and seasonality curves shown in Figure 9, the
overall temporal behavior is similar for all the identified clusters, with a
steadily increasing trend with an overall increase close to 1 standard
deviation. The most obvious difference concerns the seasonality, but it seems
to be too difficult to distinguish the main differences from Figure 9. By
analyzing the seasonality decomposition in the first and second harmonics, we
can extract the amplitude and phase variation through the following
transformation:

$A_{d}(t)=\sqrt{(\hat{g}_{2d}(t)^{2}+\hat{g}_{3d}(t)^{2})},\quad
B_{d}(t)=\arctan{(\hat{g}_{3d}(t)/\hat{g}_{2d}(t))}$

where $A_{d}(t)$ and $B_{d}(t)$ represent the functions related to the
amplitude and phase, respectively, of the $d$-th harmonic. From the data in
Figure 10, it is clear that the primary difference between the clusters is
based on the phase functions that integrate the shift information of the
seasonal component and the amplitude of the second harmonic.

##  7 Discussion

In this paper, we propose a model for clustering time series based on
quantiles. The model exploits a link between quantile regression and the
likelihood of an asymmetric Laplace distribution. The model can capture the
characteristics of data dispersion and detect the spatial extent of clusters,
taking into account different spatial membership structures. In this paper,
spatial dependence was incorporated using a Potts model of neighborhood
membership. A Bayesian framework for deriving a posteriori estimates of the
parameters is facilitated by the hierarchical structure. In this regard, it is
not difficult to extend the formulation to allow time series, within the same
cluster, to be spatially and temporally dependent by assuming space-time
correlated Gaussian errors in Equation (8). Furthermore, our model formulation
is not limited to temporal data but can be effortlessly expanded to spatial
longitudinal data with covariates.

The analysis results of the Mediterranean Sea Surface Temperature are
consistent with the literature, which classified sub-basins [33] or used
abiotic characteristics and spatial hierarchical clustering of raw data [60].
Nevertheless, our clustering method provides a more coherent spatio-temporal
classification of the water masses when compared to prior classification
attempts. Our method allows the characterization of both well-defined basins,
such as the Adriatic Sea, and clusters with small spatial extent and specific
temporal behavior, as indicated by the presence of fresh Atlantic waters in
the Strait of Gibraltar [62] or the unique conformation of the Gulf of Genoa
[63]. Moreover, the classification presented here may help to identify of
water masses most affected by the effects of climate change [64]. Although our
clustering results are based on the median to estimate the central tendency of
the time series, an important advantage of our approach is the possibility to
obtain different classifications by varying the quantile of interest according
to the objectives of the study[65].

Although the model has been adapted to cluster time series, the semi-
parametric formulation makes it suitable for application to functional data as
well. In this respect, a work closely related to ours on functional data
clustering is Kim and Oh [66]. However, the strength of our proposal is that
it accounts for all uncertainties in the clustering procedure. An additional
feature included is a flexible spatial dependence layer that allows coherent
classifications and the grading of spatial dependence strength within each
cluster.

The presented method has some limitations, mainly represented by the time
required for the MCMC iterations. In fact, in our motivating example, which
consists of 251 time series with 372 time points, it takes about 1,500 seconds
to run 3,000 iterations with $K=2$ clusters and a 31-column regressor matrix.
Using an ad hoc expectation maximization algorithm can be time-saving but with
a reduced number of features that can be included in the model specification
compared to our Bayesian framework. Another interesting aspect pertains to the
application of real data where the spatial domain is non-convex and
characterized by a water distance that may not be Euclidean due to the
presence of coasts and land between two locations. However, the specification
of anisotropic models can be challenging and context dependent. By
implementing an inter-cluster, variable spatial dependence layer, one can
achieve general adaptability in numerous real-data applications, despite the
use of a standard isotropic Potts model. Another limitation is that only one
quantile can be considered. Extending this work to a context similar to
Bondell et al. [67]’s is a research direction to be explored.

## Acknowledgement

The authors would like to thank Noémie Le Carrer for retrieving the dataset.

## Declarations

###  7.1 Funding

This research received no specific grant from any funding agency in the
public, commercial, or not-for-profit sectors.

###  7.2 Conflict of interest

The authors have no conflicts of interest to declare.

###  7.3 Code and data availability

Code and data are available upon request from the authors

## References

  *   * Schneider [2001] Schneider, S.H.: What is ’dangerous’ climate change? Nature 411, 17–19 (2001)
  * Portmann et al. [2009] Portmann, R.W., Solomon, S., Hegerl, G.C.: Spatial and seasonal patterns in climate change, temperatures, and precipitation across the United States. Proceedings of the National Academy of Sciences 106, 7324–7329 (2009)
  * Katz [2010] Katz, R.W.: Statistics of extremes in climate change. Climatic Change 100, 71–76 (2010)
  * Sun et al. [2018] Sun, F., Roderick, M.L., Farquhar, G.D.: Rainfall statistics, stationarity, and climate change. Proceedings of the National Academy of Sciences 115, 2305–2310 (2018)
  * Reich [2012] Reich, B.J.: Spatiotemporal quantile regression for detecting distributional changes in environmental processes. Journal of the Royal Statistical Society Series C: Applied Statistics 61, 535–553 (2012)
  * Cannon [2018] Cannon, A.J.: Non-crossing nonlinear regression quantiles by monotone composite quantile regression neural network, with application to rainfall extremes. Stochastic Environmental Research and Risk Assessment 32, 3207–3225 (2018)
  * Vandeskog et al. [2022] Vandeskog, S.M., Thorarinsdottir, T.L., Steinsland, I., Lindgren, F.: Quantile based modeling of diurnal temperature range with the five-parameter lambda distribution. Environmetrics 33, 2719 (2022)
  * Liao [2005] Liao, T.W.: Clustering of time series data—a survey. Pattern Recognition 38, 1857–1874 (2005)
  * Zhang and Parnell [2023] Zhang, M., Parnell, A.: Review of clustering methods for functional data. ACM Transactions on Knowledge Discovery from Data 17 (2023)
  * Grün and Leisch [2007] Grün, B., Leisch, F.: Fitting finite mixtures of generalized linear regressions in r. Computational Statistics & Data Analysis 51, 5247–5252 (2007)
  * Jiang and Serban [2012] Jiang, H., Serban, N.: Clustering random curves under spatial interdependence with application to service accessibility. Technometrics 54, 108–119 (2012)
  * Cuesta-Albertos et al. [1997] Cuesta-Albertos, J.A., Gordaliza, A., Matrán, C.: Trimmed $k$-means: an attempt to robustify quantizers. The Annals of Statistics 25, 553–576 (1997)
  * Gallegos [2002] Gallegos, M.T.: Maximum likelihood clustering with outliers. In: Jajuga, K., Sokołowski, A., Bock, H.-H. (eds.) Classification, Clustering, and Data Analysis, pp. 247–255. Springer, Berlin, Heidelberg (2002)
  * Amovin-Assagba et al. [2022] Amovin-Assagba, M., Gannaz, I., Jacques, J.: Outlier detection in multivariate functional data through a contaminated mixture model. Computational Statistics & Data Analysis 174, 107496 (2022)
  * Garcia-Escudero et al. [2015] Garcia-Escudero, L.A., Gordaliza, A., Matrán, C., Mayo-Iscar, A.: Avoiding spurious local maximizers in mixture modeling. Statistics and Computing 25, 619–633 (2015)
  * Fritz et al. [2013] Fritz, H., Garcia-Escudero, L.A., Mayo-Iscar, A.: A fast algorithm for robust constrained clustering. Computational Statistics & Data Analysis 61, 124–136 (2013)
  * Nguyen et al. [2016] Nguyen, H.D., McLachlan, G.J., Ullmann, J.F., Janke, A.L.: Spatial clustering of time series via mixture of autoregressions models and Markov random fields. Statistica Neerlandica 70, 414–439 (2016)
  * Disegna et al. [2017] Disegna, M., D’Urso, P., Durante, F.: Copula-based fuzzy clustering of spatial time series. Spatial Statistics 21, 209–225 (2017)
  * Delicado et al. [2010] Delicado, P., Giraldo, R., Comas, C., Mateu, J.: Statistics for spatial functional data: some recent contributions. Environmetrics 21, 224–239 (2010)
  * Vandewalle et al. [2022] Vandewalle, V., Preda, C., Dabo-Niang, S.: Clustering spatial functional data. In: Mateu, J., Giraldo, R. (eds.) Geostatistical Functional Data Analysis, pp. 155–174. Wiley, New York (2022)
  * Koner and Staicu [2023] Koner, S., Staicu, A.-M.: Second-generation functional data. Annual Review of Statistics and Its Application 10, 547–572 (2023)
  * Jiang and Serban [2012] Jiang, H., Serban, N.: Clustering random curves under spatial interdependence with application to service accessibility. Technometrics 54, 108–119 (2012)
  * Hu et al. [2022] Hu, G., Geng, J., Xue, Y., Sang, H.: Bayesian spatial homogeneity pursuit of functional data: an application to the us income distribution. Bayesian Analysis 1, 1–27 (2022)
  * Secchi et al. [2013] Secchi, P., Vantini, S., Vitelli, V.: Bagging Voronoi classifiers for clustering spatial functional data. International Journal of Applied Earth Observation and Geoinformation 22, 53–64 (2013)
  * Giraldo et al. [2012] Giraldo, R., Delicado, P., Mateu, J.: Hierarchical clustering of spatially correlated functional data. Statistica Neerlandica 66, 403–421 (2012)
  * Bera et al. [2016] Bera, A.K., Galvao Jr, A.F., Montes-Rojas, G.V., Park, S.Y.: Asymmetric Laplace regression: maximum likelihood, maximum entropy and quantile regression. Journal of Econometric Methods 5, 79–101 (2016)
  * Potts [1952] Potts, R.B.: Some generalized order-disorder transformations. Mathematical Proceedings of the Cambridge Philosophical Society 48, 106–109 (1952)
  * Nykjaer [2009] Nykjaer, L.: Mediterranean Sea surface warming 1985–2006. Climate Research 39, 11–17 (2009)
  * Lejeusne et al. [2010] Lejeusne, C., Chevaldonné, P., Pergent-Martini, C., Boudouresque, C.F., Pérez, T.: Climate change effects on a miniature ocean: the highly diverse, highly impacted Mediterranean Sea. Trends in Ecology & Evolution 25, 250–260 (2010)
  * Bethoux et al. [1990] Bethoux, J.-P., Gentili, B., Raunet, J., Tailliez, D.: Warming trend in the western Mediterranean deep water. Nature 347, 660–662 (1990)
  * Ibrahim et al. [2021] Ibrahim, O., Mohamed, B., Nagy, H.: Spatial variability and trends of marine heat waves in the eastern Mediterranean Sea over 39 years. Journal of Marine Science and Engineering 9, 643 (2021)
  * Nunes et al. [2019] Nunes, S., Perez, G.L., Latasa, M., Zamanillo, M., Delgado, M., Ortega-Retuerta, E., Marrasé, C., Simó, R., Estrada, M.: Size fractionation, chemotaxonomic groups and bio-optical properties of phytoplankton along a transect from the Mediterranean Sea to the SW Atlantic Ocean. Scientia Marina 83, 87–109 (2019)
  * Pastor et al. [2019] Pastor, F., Valiente, J.A., Palau, J.L.: Sea surface temperature in the mediterranean: Trends and spatial patterns (1982–2016). Meteorology and Climatology of the Mediterranean and Black Seas 175, 297–309 (2019)
  * McLachlan and Peel [2000] McLachlan, G.J., Peel, D.: Finite Mixture Models. Wiley, New York (2000)
  * Strauss [1977] Strauss, D.J.: Clustering on coloured lattices. Journal of Applied Probability 14, 135–143 (1977)
  * Gaetan et al. [2017] Gaetan, C., Girardi, P., Pastres, R.: Spatial clustering of curves with an application of satellite data. Spatial Statistics 20, 110–124 (2017)
  * Koenker and Bassett [1978] Koenker, R., Bassett, G.: Regression quantiles. Econometrica 46, 33–50 (1978)
  * Koenker and Machado [1999] Koenker, R., Machado, J.A.: Goodness of fit and related inference processes for quantile regression. Journal of the American Statistical Association 94, 1296–1310 (1999)
  * Huber [1981] Huber, P.J.: Robust Statistics. Wiley, New York (1981)
  * Besag [1986] Besag, J.: On the statistical analysis of dirty pictures. Journal of the Royal Statistical Society: Series B (Methodological) 48, 259–279 (1986)
  * Cucala et al. [2009] Cucala, L., Marin, J.-M., Robert, C.P., Titterington, D.M.: A Bayesian reassessment of nearest-neighbor classification. Journal of the American Statistical Association 104, 263–273 (2009)
  * Benoit and Van den Poel [2017] Benoit, D.F., Van den Poel, D.: bayesQR: A Bayesian approach to quantile regression. Journal of Statistical Software 76, 1–32 (2017)
  * Robert and Casella [2004] Robert, C.P., Casella, G.: Monte Carlo Statistical Methods. Springer, New York (2004)
  * Kotz et al. [2001] Kotz, S., Kozubowski, T., Podgorski, K.: The Laplace Distribution and Generalizations: A Revisit with Applications to Communications, Economics, Engineering, and Finance. Springer, New York (2001)
  * Kozumi and Kobayashi [2011] Kozumi, H., Kobayashi, G.: Gibbs sampling methods for Bayesian quantile regression. Journal of Statistical Computation and Simulation 81, 1565–1578 (2011)
  * Jorgensen [1982] Jorgensen, B.: Statistical Properties of the Generalized Inverse Gaussian Distribution. Springer, New York (1982)
  * Marjoram et al. [2003] Marjoram, P., Molitor, J., Plagnol, V., Tavaré, S.: Markov chain Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences 100, 15324–15328 (2003)
  * Marin et al. [2012] Marin, J.-M., Pudlo, P., Robert, C.P., Ryder, R.J.: Approximate Bayesian computational methods. Statistics and Computing 22, 1167–1180 (2012)
  * Pereyra et al. [2013] Pereyra, M., Dobigeon, N., Batatia, H., Tourneret, J.-Y.: Estimating the granularity coefficient of a Potts-Markov random field within a Markov chain Monte Carlo algorithm. IEEE Transactions on Image Processing 22, 2385–2397 (2013)
  * Sottile and Adelfio [2019] Sottile, G., Adelfio, G.: Clusters of effects curves in quantile regression models. Computational Statistics 34, 551–569 (2019)
  * Hubert and Arabie [1985] Hubert, L., Arabie, P.: Comparing partitions. Journal of Classification 2, 193–218 (1985)
  * Kaufman and Rousseeuw [2009] Kaufman, L., Rousseeuw, P.J.: Finding Groups in Data: an Introduction to Cluster Analysis. John Wiley & Sons, New York (2009)
  * Romary et al. [2015] Romary, T., Ors, F., Rivoirard, J., Deraisme, J.: Unsupervised classification of multivariate geostatistical data: Two algorithms. Computers & Geosciences 85, 96–103 (2015)
  * Wang and Xu [2017] Wang, X.-F., Xu, Y.: Fast clustering using adaptive density peak detection. Statistical Methods in Medical Research 26, 2800–2811 (2017)
  * Bouveyron et al. [2015] Bouveyron, C., Côme, E., Jacques, J.: The discriminative functional mixture model for a comparative analysis of bike sharing systems. The Annals of Applied Statistics 9, 1726–1760 (2015)
  * Bouveyron and Jacques [2011] Bouveyron, C., Jacques, J.: Model-based clustering of time series in group-specific functional subspaces. Advances in Data Analysis and Classification 5, 281–300 (2011)
  * Oliver and Webster [1989] Oliver, M., Webster, R.: A geostatistical basis for spatial weighting in multivariate classification. Mathematical Geology 21, 15–35 (1989)
  * Eilers et al. [2008] Eilers, P.H., Gampe, J., Marx, B.D., Rau, R.: Modulation models for seasonal time series and incidence tables. Statistics in Medicine 27, 3430–3441 (2008)
  * Watanabe and Opper [2010] Watanabe, S., Opper, M.: Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research 11, 3571–3594 (2010)
  * de la Hoz et al. [2018] de la Hoz, C.F., Ramos, E., Puente, A., Méndez, F., Menéndez, M., Juanes, J.A., Losada, Í.J.: Ecological typologies of large areas. an application in the Mediterranean Sea. Journal of Environmental Management 205, 59–72 (2018)
  * D’Ortenzio and Ribera d’Alcalà [2008] D’Ortenzio, F., Ribera d’Alcalà, M.: On the trophic regimes of the Mediterranean Sea: a satellite analysis. Biogeosciences Discussions 5, 139–148 (2008)
  * Bouzinac et al. [2003] Bouzinac, C., Font, J., Johannessen, J.: Annual cycles of sea level and sea surface temperature in the western Mediterranean Sea. Journal of Geophysical Research: Oceans 108(C3), 3059 (2003)
  * Cutroneo and Capello [2023] Cutroneo, L., Capello, M.: The cold waters in the port of Genoa (NW Mediterranean Sea) during the marine heatwave in summer 2022. Journal of Marine Science and Engineering 11, 1568 (2023)
  * Shaltout and Omstedt [2014] Shaltout, M., Omstedt, A.: Recent sea surface temperature trends and future scenarios for the Mediterranean Sea. Oceanologia 56, 411–443 (2014)
  * Cade and Noon [2003] Cade, B.S., Noon, B.R.: A gentle introduction to quantile regression for ecologists. Frontiers in Ecology and the Environment 1, 412–420 (2003)
  * Kim and Oh [2020] Kim, J., Oh, H.-S.: Pseudo-quantile functional data clustering. Journal of Multivariate Analysis 178, 104626 (2020)
  * Bondell et al. [2010] Bondell, H.D., Reich, B.J., Wang, H.: Noncrossing quantile regression curve estimation. Biometrika 97, 825–838 (2010)

Supplementary material for
“Spatial quantile clustering of climate data”

##  Appendix A Conditional distributions in the MCMC algorithm

We consider the full set of observations
$\boldsymbol{y}=\\{\boldsymbol{y}_{1}^{\prime},\ldots,\boldsymbol{y}_{N}^{\prime}\\}$
of $N$ vectors $\boldsymbol{y}_{i}=(y_{i1},...,y_{i\,T})^{\prime}$ of
independent observations. Each observation comes from the model

$\displaystyle
y_{it}=\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}+\tau\sigma_{c_{i}}w_{it}+\omega\sigma_{c_{i}}\sqrt{w_{it}}\,\nu_{it}$

where $\tau=(1-2p)/\\{p(1-p)\\}$, $\omega^{2}={2}/\\{p(1-p)\\}$ where
$\nu_{it}\sim\mathcal{N}(0,1)$ and $w_{it}\sim\mathcal{E}(1)$ are mutually
independent random variables.

The pmf of the random vector
$\boldsymbol{c}=(\boldsymbol{c}_{1},\ldots,\boldsymbol{c}_{N})^{\prime}$ (a
Potts model) is given by

$\displaystyle f(\boldsymbol{c}|\boldsymbol{\alpha},\beta)$ $\displaystyle=$
$\displaystyle\frac{\exp\left\\{\sum_{k=1}^{K}\left(\alpha_{k}n_{k}+\beta_{k}\sum_{i=1}^{N}n_{k,i}/2\right)\right\\}}{Z(\boldsymbol{\alpha},\boldsymbol{\beta})},$

where $n_{k,i}=\sum_{j\in\partial_{i}}I(c_{i}=k)I(c_{j}=k)$ is the number of
the neighbors of the site $i$ with membership equal to $k$.

The joint density of $\boldsymbol{y}_{i}$ given
$\boldsymbol{w}_{i}=(w_{i1},...,w_{iT})^{\prime}$ and $c_{i}$ is

$\displaystyle
f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},\boldsymbol{w}_{i})\propto{(\sigma_{c_{i}}^{2})}^{-{T}/{2}}\prod_{t=1}^{T}w_{it}^{-{1}/{2}}\exp\left\\{-\frac{1}{2}\sum_{t=1}^{T}\frac{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}-\tau\sigma_{c_{i}}w_{it})^{2}}{\sigma^{2}_{c_{i}}\omega^{2}w_{it}}\right\\}$

with $\boldsymbol{\psi}(t)=(\psi_{0}(t),\ldots,\psi_{J}(t)^{\prime}$ and
$\boldsymbol{\gamma}_{c_{i}}=(\gamma_{c_{i},0},\ldots,\gamma_{c_{i},J})^{\prime}$.

We want to make inference on the model parameters
$\boldsymbol{\theta}=(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma}^{\prime},\boldsymbol{\sigma^{2}}^{\prime})^{\prime}$,
with
$\boldsymbol{\gamma}=(\boldsymbol{\gamma}_{1}^{\prime},\ldots,\boldsymbol{\gamma}_{K}^{\prime})^{\prime}$
and $\boldsymbol{\sigma^{2}}=(\sigma^{2}_{1},\ldots,\sigma^{2}_{K})$.

The prior distribution for $\boldsymbol{\theta}$ is given by
$\pi(\boldsymbol{\theta})=\pi(\boldsymbol{\alpha})\times\pi({\boldsymbol{\beta}})\times\pi(\boldsymbol{\gamma})\times\pi(\boldsymbol{\sigma}^{2})$
where
$\pi(\boldsymbol{\alpha})=\prod_{k=2}^{K}\pi({\alpha}_{k})\propto\prod_{k=2}^{K}\exp\left\\{-\frac{\alpha_{k}^{2}}{2a}\right\\}$,
$a>0$;

$\pi(\boldsymbol{\beta})={\color[rgb]{0,0,0}\prod_{k=1}^{K}I_{(0,b)}(\beta_{k})}$,
$b>0$

$\pi(\boldsymbol{\gamma})=\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k})\propto\prod_{k=1}^{K}\exp\left\\{-\frac{1}{2g}\boldsymbol{\gamma}_{k}^{\prime}\boldsymbol{\gamma}_{k}\right\\}$,
$g>0$;

$\pi(\boldsymbol{\sigma}^{2})=\prod_{k=1}^{K}\pi(\sigma_{k}^{2})\propto\prod_{k=1}^{K}{(\sigma^{2}_{k})}^{-s_{0}-1}\exp\left(-{d_{0}}/{\sigma^{2}_{k}}\right)$,
$s_{0}>0$ and $d_{0}>0$

Applying the Bayes theorem we obtain the posterior distribution as

$\displaystyle\pi(\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\gamma},\boldsymbol{\sigma^{2}},\boldsymbol{w}|\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle\prod_{i=1}^{N}f(\boldsymbol{y}_{i}|c_{i},\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},\boldsymbol{w}_{i})f(\boldsymbol{w}_{i})\times
f(\boldsymbol{c}|\boldsymbol{\alpha},\boldsymbol{\beta})\times$
$\displaystyle\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k})\times\pi(\sigma^{2}_{k})\times\pi(\boldsymbol{\alpha})\times\pi(\boldsymbol{\beta})$

#### Conditional distribution for $\boldsymbol{\gamma}$

We have
$\pi(\boldsymbol{\gamma}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{c},\boldsymbol{\sigma^{2}},\boldsymbol{w},\boldsymbol{y})=\prod_{k=1}^{K}\pi(\boldsymbol{\gamma}_{k}|{\sigma}^{2}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$,
where the full conditional posterior distribution of the quantile regression
parameters $\boldsymbol{\gamma}_{c_{i}}$ is obtained as follows

$\displaystyle\pi(\boldsymbol{\gamma}_{k}|{\sigma}^{2}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2g}\boldsymbol{\gamma}_{k}^{\prime}\boldsymbol{\gamma}_{k}\right\\}\times$
$\displaystyle\prod_{i=1}^{N}\left[\exp\left\\{-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}}\sum_{t=1}^{T}\frac{\left(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}\right)^{2}}{w_{it}}\right\\}\right]^{I(c_{i}=k)}$

with $u_{it}=y_{it}-\tau\sigma_{c_{i}}w_{it}$.

The previous formula can be written in compact form using the matrix notation,
namely

$\displaystyle\pi(\boldsymbol{\gamma}_{k}|\sigma^{2}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2g}\boldsymbol{\gamma}_{k}^{\prime}\boldsymbol{\gamma}_{k}\right\\}\times\exp\left\\{-\frac{1}{2}\sum_{i=1}^{N}(\boldsymbol{u}_{i}-\boldsymbol{\Psi}\boldsymbol{\gamma}_{c_{i}})^{\prime}\boldsymbol{W}_{i}(\boldsymbol{u}_{i}-\boldsymbol{\Psi}\boldsymbol{\gamma}_{c_{i}})I(c_{i}=k)\right\\}$
$\displaystyle\propto$
$\displaystyle\exp\left\\{-\frac{1}{2}\left[\boldsymbol{\gamma}_{k}^{\prime}\left(\frac{1}{g}\boldsymbol{I}+\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{\Psi}I(c_{i}=k)\right)\boldsymbol{\gamma}_{k}-\right.\right.$
$\displaystyle\left.\left.2\boldsymbol{\gamma}_{k}^{\prime}\left(\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{u}_{i}I(c_{i}=k)\right)\right]\right\\}$

Here $\boldsymbol{u}_{i}=(u_{i1},\ldots,u_{iT})^{\prime}$,
$\boldsymbol{W}_{i}$ is a diagonal matrix with entries
$[\sigma^{2}_{c_{i}}\omega^{2}w_{it}]^{-1}$,
$\boldsymbol{\Psi}=(\boldsymbol{\psi}(1),\ldots,\boldsymbol{\psi}(T))^{\prime}$
is the matrix of covariates. Thus, the conditional posterior density of
$\boldsymbol{\gamma}_{k}$ is proportional to the density of a multivariate
Gaussian distribution with vector mean

$\boldsymbol{m}_{k}=\left(\frac{1}{g}\boldsymbol{I}+\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{\Psi}I(c_{i}=k)\right)^{-1}\left(\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\mathbf{u}_{i}I(c_{i}=k)\right)$

and covariance matrix

$\boldsymbol{S}_{k}=\left(\frac{1}{g}\boldsymbol{I}+\sum_{i=1}^{N}\boldsymbol{\Psi}^{\prime}\boldsymbol{W}_{i}\boldsymbol{\Psi}I(c_{i}=k)\right)^{-1}$

#### Full conditional for $\boldsymbol{\sigma}^{2}$

We note that

$\displaystyle\pi(\boldsymbol{\sigma}^{2}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})\propto\prod_{k=1}^{K}\pi(\sigma^{2}_{k}|\boldsymbol{\gamma}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$

To obtain the full conditional posterior of $\sigma^{2}_{k}$ we proceed as
follows

$\displaystyle\pi(\sigma_{k}^{2}|\boldsymbol{\gamma}_{k},\boldsymbol{c},\boldsymbol{w},\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle(\sigma^{2}_{k})^{-s_{0}-1}\exp\left(-\frac{d_{0}}{\sigma^{2}_{k}}\right)\times$
(1)
$\displaystyle\prod_{i=1}^{N}\left[{(\sigma^{2}_{c_{i}})}^{-{T}/{2}}\exp\left\\{-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}}\sum_{t=1}^{T}\frac{\left(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}\right)^{2}}{w_{it}}\right\\}\right]^{I(c_{i}=k)}$
$\displaystyle\propto$
$\displaystyle{(\sigma_{k}^{2})}^{-(s_{0}+n_{k}T/2)-1}\times$
$\displaystyle\exp\left[-\frac{1}{\sigma^{2}_{k}}\left[d_{0}+\frac{1}{2\omega^{2}}\sum_{i=1}^{N}\left\\{\sum_{t=1}^{T}\frac{\left(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}\right)^{2}}{w_{it}}\right\\}I(c_{i}=k)\right]\right]$

where $n_{k}=\sum_{i=1}^{N}I(c_{i}=k)$, which represents the number of vectors
of $\boldsymbol{y}$ with membership $k$.

Therefore, (1) is the distribution of an inverse Gamma random variable with
shape parameter and scale parameter,

$a=s_{0}+n_{k}T/2,\qquad
b=d_{0}+\frac{1}{2\omega^{2}}\sum_{i=1}^{N}\left\\{\sum_{t=1}^{T}\frac{\left(u_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}\right)^{2}}{w_{it}}\right\\}I(c_{i}=k).$

We note that random samples from inverse Gamma distribution can be drawn from
a Gamma distribution observing the relationship that if
$G\sim\mathrm{Gamma}(a,1/b)$ then $1/G\sim\mathrm{InvGamma}(a,b)$.

#### Conditional distribution for $\boldsymbol{w}$

We note that

$\displaystyle\pi(\boldsymbol{w}|\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\sigma}^{2},\boldsymbol{c},\boldsymbol{y})$
$\displaystyle\propto$
$\displaystyle\prod_{i=1}^{N}f(\boldsymbol{w}_{i}|\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},c_{i},\boldsymbol{y}_{i})=\prod_{i=1}^{N}\prod_{t=1}^{T}f({w}_{it}|\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},c_{i},y_{it})$

and

$\displaystyle
f({w}_{it}|\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},c_{i},y_{it})$
$\displaystyle\propto$ $\displaystyle
f(y_{it}|\boldsymbol{\gamma}_{c_{i}},\sigma^{2}_{c_{i}},c_{i},w_{it})\times
f(w_{it})$ $\displaystyle\propto$ $\displaystyle
w_{it}^{-{1}/{2}}\exp\left\\{-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}w_{it}}\left(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}}-\tau\sigma_{c_{i}}w_{it}\right)^{2}-w_{it}\right\\}$
$\displaystyle\propto$ $\displaystyle
w_{it}^{-{1}/{2}}\exp\left[-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}w_{it}}\left\\{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}+\tau^{2}\sigma^{2}_{c_{i}}w_{it}^{2}\right\\}-w_{it}\right]$
$\displaystyle\propto$ $\displaystyle
w_{it}^{-{1}/{2}}\exp\left\\{-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}w_{it}}(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}-\frac{\tau^{2}w_{it}}{2\omega^{2}}-w_{it}\right\\}$
$\displaystyle\propto$ $\displaystyle
w_{it}^{-{1}/{2}}\exp\left\\{-\frac{1}{2\sigma^{2}_{c_{i}}\omega^{2}w_{it}}(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}-\frac{w_{it}}{2\omega^{2}}(\tau^{2}+2\omega^{2})\right\\}$
$\displaystyle\propto$ $\displaystyle
w_{it}^{-{1}/{2}}\exp\left[-\frac{1}{2}\left\\{\frac{\tau^{2}+2\omega^{2}}{\omega^{2}}w_{it}+\frac{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}}{\sigma^{2}_{c_{i}}\omega^{2}w_{it}}\right\\}\right]$

This expression resembles a Generalized Inverse Gaussian distribution
$\operatorname{GIG}(d,a,b)$ with the density

$\displaystyle f(w;a,b,d)$ $\displaystyle\propto$ $\displaystyle
w^{(d-1)}\exp\left\\{-\frac{1}{2}\left(aw+\frac{b}{w}\right)\right\\}$

where $d=0.5$, $\displaystyle a=\frac{\tau^{2}+2\omega^{2}}{\omega^{2}}$,
$\displaystyle
b=\frac{(y_{it}-\boldsymbol{\psi}(t)^{\prime}\boldsymbol{\gamma}_{c_{i}})^{2}}{\sigma^{2}_{c_{i}}\omega^{2}}$.

##  Appendix B Selection of the number of elements for the B-spline basis

Table 1: Overall $\operatorname{BIC}$ for different values of $J_{1}$, $J_{2}$, and $J_{3}$. | $J_{3}$
---|---
$J_{2}=4$ | 4 | 5 | 6 | 7 | 8
$J_{1}=4$ | 5129 | 7626 | 9791 | 12092 | 14527
$J_{1}=5$ | 5796 | 8223 | 10419 | 12664 | 15095
$J_{1}=6$ | 4802 | 7179 | 9348 | 11450 | 13638
$J_{1}=7$ | 5074 | 7494 | 9632 | 11721 | 13852
$J_{1}=8$ | 5614 | 8042 | 10190 | 12315 | 14459
| $J_{3}$
$J_{2}=5$ | 4 | 5 | 6 | 7 | 8
$J_{1}=4$ | -9044 | -6861 | -5066 | -2993 | -737
$J_{1}=5$ | -8385 | -6204 | -4435 | -2416 | -120
$J_{1}=6$ | -9998 | -7865 | -6025 | -4044 | -1836
$J_{1}=7$ | -9996 | -7847 | -5990 | -4094 | -2105
$J_{1}=8$ | -9555 | -7426 | -5575 | -3723 | -1745
| $J_{3}$
$J_{2}=6$ | 4 | 5 | 6 | 7 | 8
$J_{1}=4$ | -15650 | -13529 | -11702 | -9573 | -7432
$J_{1}=5$ | -14936 | -12847 | -11019 | -8943 | -6807
$J_{1}=6$ | -16313 | -14273 | -12358 | -10281 | -8277
$J_{1}=7$ | -16492 | -14460 | -12573 | -10668 | -8736
$J_{1}=8$ | -16077 | -14065 | -12162 | -10251 | -8358
| $J_{3}$
$J_{2}=7$ | 4 | 5 | 6 | 7 | 8
$J_{1}=4$ | -21116 | -19031 | -17445 | -15082 | -12877
$J_{1}=5$ | -20299 | -18240 | -16678 | -14304 | -12083
$J_{1}=6$ | -21213 | -19147 | -17541 | -15076 | -12865
$J_{1}=7$ | -21716 | -19716 | -18138 | -15778 | -13607
$J_{1}=8$ | -21459 | -19451 | -17934 | -15579 | -13423
| $J_{3}$
$J_{2}=8$ | 4 | 5 | 6 | 7 | 8
$J_{1}=4$ | -23877 | -21808 | -20437 | -18240 | -16081
$J_{1}=5$ | -23084 | -21036 | -19738 | -17507 | -15358
$J_{1}=6$ | -24149 | -22158 | -20856 | -18480 | -16268
$J_{1}=7$ | -24547 | -22643 | -21478 | -19185 | -17028
$J_{1}=8$ | -24417 | -22532 | -21357 | -19095 | -16945

##  Appendix C Selection of number of clusters

|
---|---
|
Figure 1: Silhouette profile for HC procedures and $\operatorname{WAIC}$
values for ALC-S method with 4 and 8 nearest neighbors by number of clusters.