funOCLUST: Clustering Functional Data with Outliers

Clustering is a subtype of classification, which aims to group similar sets of observations together without any prior knowledge of the cluster memberships. While non-parametric methods exist, e.g., hierarchical and k-means clustering, model-based clustering employs a parametric approach. In general, the density of a finite mixture model is

$$f(\mathbf{x}\mid\bm{\vartheta})=\sum_{g=1}^{G}\pi_{g}f_{g}(\mathbf{x}\mid\bm{\theta}_{g}),$$

(1)

where $\bm{\vartheta}=\{\pi_{1},\dots,\pi_{G},\bm{\theta}_{1},\dots\bm{\theta}_{G}\}$, $\pi_{g}>0$ is the $g$th mixing proportion with $\sum_{g=1}^{G}\pi_{g}=1$, and $f_{g}(\mathbf{x}\mid\bm{\theta}_{g})$ is the $g$th component density with parameters $\bm{\theta}_{g}$. Typically, each component corresponds to a cluster (see McNicholas, 2016, for a discussion). While each cluster can be modelled with various component distributions, the Gaussian distribution remains popular due to its simplicity. The parameters and cluster memberships are estimated by maximizing the log-likelihood, often with the expectation-maximization (EM) algorithm (Dempster et al., 1977).

Clark and McNicholas (2024) develop the OCLUST algorithm, which simultaneously clusters data and trims outliers. Data are modelled with mixtures of Gaussian distributions, and outliers are trimmed iteratively one-by-one. A subset log-likelihood is defined to be the log-likelihood of the data with a single data point removed, and there are $n$ such subsets. After each outlier removal, the distribution of the subset log-likelihoods is compared to the mixture of shifted and scaled beta distributions derived in the paper. Specifically, let $l_{\mathcal{X}}$ be the complete-data log-likelihood of the entire dataset, let $l_{\mathcal{X}\setminus\mathbf{x}_{j}}$ be the complete-data log-likelihood of the $j$th subset, i.e., with observation $\mathbf{x}_{j}$ removed, and let $D_{j}=l_{\mathcal{X}\setminus\mathbf{x}_{j}}-l_{\mathcal{X}}$. Then, if $\mathbf{x}_{j}$ belongs to the $h$th cluster, i.e. with indicator variable $z_{jh}=1$,

$$D_{j}\mid(z_{jh}=1)\sim f_{\text{beta}}\left(\frac{2n_{h}}{(n_{h}-1)^{2}}(d_{j}-c)~\bigg{|}~\frac{p}{2},\frac{n_{h}-p-1}{2}\right)$$

(2)

for $c<d_{j}<\frac{(n_{h}-1)^{2}}{2n_{h}}+c,n_{h}>p+1$, where $c=-\log\hat{\pi}_{h}+\frac{p}{2}\log(2\pi)+\frac{1}{2}\log\lvert\mathbf{S}_{h}\rvert$, $n_{h}$ is the number of points in cluster $h$, $\hat{\pi}_{h}=n_{h}/n$,

$$\mathbf{S}_{h}=\frac{1}{n_{h}-1}\sum_{i=1}^{n}z_{ih}(\mathbf{x}_{i}-\bar{\mathbf{x}}_{h})(\mathbf{x}_{i}-\bar{\mathbf{x}}_{h})^{\prime}$$

is the sample covariance matrix of cluster $h$, and $\bar{\mathbf{x}}_{h}=\frac{1}{n_{h}}\sum_{i=1}^{n}z_{ih}\mathbf{x}_{i}.$ The trimmed model chosen minimizes the Kullback-Leibler (KL) divergence between the observed differences and the derived distribution.

In multivariate data analysis, a random sample ($\mathbf{x}_{1},\ldots\mathbf{x}_{n})$ is typically a collection of vectors in $\mathbb{R}^{p}$. In functional data, a random sample of functions ($Y_{1}(t),\ldots Y_{n}(t)),t\in\mathcal{T}\subset\mathbb{R}$ exists in infinite-dimensional space. In practice, however, these functions are evaluated at a series of discrete values $\{t_{i1},\ldots,t_{ij_{i}}\}$, resulting in observations of the form $(\mathbf{t}_{i},\mathbf{y}_{i})$, where $\mathbf{t}_{i}=(t_{i1},\ldots,t_{ij_{i}})$ is a vector of points in the domain and $\mathbf{y}_{i}=(y_{i1},\ldots,y_{ij_{i}})$ is the vector of observed values at the points $\mathbf{t}_{i}$. To return to infinite-dimensional space, we can reconstruct the functional form. A popular choice is a basis expansion, which is well explained by Ramsay and Silverman (2005).

In a basis expansion, the function $Y(t)$ is approximated by a linear combination of basis functions so that $\hat{Y}(t)=\sum_{k}\beta_{k}B_{k}(t)$, where $B_{k}(t)$ is the value at $t$ for the $k$th basis function. Basis functions must be independent from each other and must be able to approximate any function with reasonable accuracy. Some choices for basis functions are Fourier series, wavelet functions, power series, exponential, polynomial and power bases, as well as B-spline basis functions. Functional principal components analysis (FPCA) is a special type of basis expansion, where $Y(t)$ is approximated by a linear combination of mutually orthogonal and normalized weight functions that explain decreasing amounts of variation in the sample ($Y_{1}(t),\ldots Y_{n}(t))$.

Jacques and Preda (2014) organize functional data clustering into three paradigms — two-stage methods, nonparametric methods and model-based methods. In model-based approaches, the filtering and clustering steps are completed simultaneously. James and Sugar (2003) introduce a Gaussian mixture model for functional clustering by assuming that the coefficients in the basis expansion are normally-distributed. These coefficients are not fixed– they are estimated within the EM algorithm. Similarly, Jacques and Preda (2013) generate a Gaussian mixture model with FPCA scores. Bouveyron and Jacques (2011) combine the model of James and Sugar (2003) with the assumption that each cluster lies in a lower-dimensional subspace, thus reducing the dimension when there are large numbers of basis functions. The resultant algorithm is called funHDDC.

A similar, more general, approach is to treat the functions as finite mixtures of regressions:

$$Y_{i}(t)=\mathbf{\mbox{\boldmath$\beta$}}_{g}\mathfrak{X}+\sigma_{g}E_{i}(t),$$

(3)

where $\mathfrak{X}$ is a regression matrix, which may be constructed from polynomials, splines, or B-splines, etc. Extensions to this method are regularized regression mixtures and the regression mixtures with mixed-effects (Chamroukhi and Nguyen, 2019).

Abraham et al. (2003) propose a two-stage approach, first filtering the data with a B-spline basis and then clustering the resulting coefficients with the $k$-means algorithm. In this case, because the filtering is completed first, the coefficients are fixed. Similarly, Peng and Müller (2008) first use FPCA to decompose the functions. The FPCA scores are then scaled and clustered using $k$-means clustering. Nguyen et al. (2018) propose a two-part algorithm, first fitting the B-spline basis coefficients to the functions, and then applying a Gaussian mixture model to those coefficients. They suppose that $\beta$ is a random variable described by a Gaussian mixture model. Then the distribution of the fitted coefficients, $\hat{\mathbf{\mbox{\boldmath$\beta$}}}_{i}$ can be described as a mixture model with density

$$f_{\hat{\mathbf{\mbox{\boldmath$\beta$}}}}(\mathbf{b}_{i}\mid\bm{\vartheta})=\sum_{g=1}^{G}\pi_{g}\phi(\mathbf{b}\mid\bm{\mu}_{g},\bm{\Sigma}_{g}+\sigma^{2}(\mathbf{B}_{i}^{\top}\mathbf{B}_{i})^{-1}).$$

(4)

Hubert et al. (2015) separate functional outliers into two types: isolated outliers and persistent outliers. While isolated outliers are atypical in a narrow region, e.g., a spike, persistent outliers affect a larger portion of the domain, e.g., shift, amplitude, or shape. In functional data clustering, algorithms that handle outliers also belong to the three paradigms in Section 1.3. Using a two-stage approach, Garcia-Escudero and Gordaliza (2005) propose a method to cluster functional data with outliers by first decomposing the functions using a B-spline basis, and then clustering the resulting coefficients with trimmed $k$-means clustering.

In the model-based realm, funHDDC has been extended to be more robust to outliers. Amovin-Assagba et al. (2022) develop C-funHDDC, which combines a contaminated mixture model with funHDDC by assuming that the functions reside in a lower-dimensional cluster-specific subspace. Similarly, Anton and Smith (2023), use a t-distribution with the funHDDC foundation to generate a more robust model in the presence of outliers, called T-funHDDC.

In a non-parametric approach, a functional isolation forest (FIF; Staerman et al., 2019) is a collection of isolation trees which randomly split the functional space into partitions, iteratively, until each function is isolated. Functions isolated in fewer iterations tend to be more outlying. An outlier score is calculated using a forest of these isolation trees. Finally, Hubert et al. (2015) determine functional outlyingness by first calculating outlyingness at each time point and then taking the weighted average to create an outlier score for each curve.

The OCLUST algorithm hinges on the the fact that log-likelihoods of the data subsets have an approximate shifted and scaled beta distribution. Extending the OCLUST algorithm to functional data will require something similar: a model with an associated log-likelihood function that can be evaluated on data subsets, and for which the distribution of these subset log-likelihoods can be derived.

While functions exist in infinite-dimensional space, they are observed at a series of discrete points $(t_{1},...,t_{j})$. If each curve is sampled at the same points, then the data for each curve can be considered as a vector. If these vectors arise from a finite Gaussian mixture, then there will be an associated log-likelihood for the data and associated parameters $\pi_{g},\bm{\mu}_{g},\bm{\Sigma}_{g}$, for each component. An obvious choice would be to apply the OCLUST algorithm to the discretized values when modelling a finite Gaussian mixture model of functional data with outliers. However, these vectors will be the same dimension as the number of points at which each function is sampled. To retain enough information about the curve, the vectors of observations are likely to be very high-dimensional.

Instead, consider applying a filtering approach to the raw data, specifically their decomposition with cubic B-splines,

$$Y_{i}(t)=\sum_{k=1}^{K+4}\beta_{i,k}B(t)+\varepsilon_{i}=\mathbf{B}(t)\mathbf{\mbox{\boldmath$\beta$}}_{i}+\varepsilon_{i},$$

(5)

where $\mathbf{B}(t)$ is the cubic B-spline basis with $K$ interior knots evaluated at point $t$, $\mathbf{\mbox{\boldmath$\beta$}}_{i}$ is the vector of coefficients for the basis expansion of function $Y_{i}(t)$, and $\varepsilon_{i}\sim N(0,\sigma^{2})$ is the error term. This retains the infinite-dimensional properties of the functions, and each function has a unique $(K+4)$-dimensional representation given by the vector $\mathbf{\mbox{\boldmath$\beta$}}_{i}$. In effect, each $Y_{i}(t)$ is transformed to a vector in $\mathbb{R}^{K+4}$, the dimension of which is governed only by the choice in number of interior knots. Transforming $Y_{i}(t)$ maintains the principle of working with functions while creating a vector representation in a much more manageable dimension. If we assume that the vectors $\mathbf{\mbox{\boldmath$\beta$}}_{i}$ come from a multivariate normal distribution, then the extension of the OCLUST algorithm to functional data becomes straightforward.

Consider a sample of i.i.d. functions ($Y_{1}(t),\ldots,Y_{n}(t)$). Suppose that each $Y_{i}(t)$ can be decomposed into a linear combination of cubic B-splines, i.e., $Y_{i}(t)=\mathbf{B}(t)\mathbf{\mbox{\boldmath$\beta$}}_{i}+\varepsilon_{i}$, where $\varepsilon_{i}\sim N(0,\sigma^{2})$. To capture heterogeneity across functional observations and enable clustering of similar functional patterns, Nguyen et al. (2018) propose that each $\mathbf{\mbox{\boldmath$\beta$}}_{i}$ arise from a finite Gaussian mixture model with density

$$f_{\mathbf{\mbox{\boldmath$\beta$}}}(\mathbf{b}\mid\bm{\vartheta})=\sum_{g=1}^{G}\pi_{g}\phi(\mathbf{b}\mid\bm{\mu}_{g},\bm{\Sigma}_{g}).$$

(6)

In addition, they show that the fitted coefficients ($\hat{\mathbf{\mbox{\boldmath$\beta$}}}_{1},\ldots,\hat{\mathbf{\mbox{\boldmath$\beta$}}}_{n}$), corresponding to functions ($Y_{1}(t),\ldots,Y_{n}(t)$) follow their own Gaussian mixture model with density

$$f_{\hat{\mathbf{\mbox{\boldmath$\beta$}}}}(\mathbf{b}\mid\bm{\vartheta})=\sum_{g=1}^{G}\pi_{g}\phi(\mathbf{b}\mid\bm{\mu}_{g},\bm{\Sigma}_{g}+\sigma^{2}(\mathbf{B}^{\top}\mathbf{B})^{-1}),$$

(7)

where $\mathbf{B}=\mathbf{B}(\mathbf{t})$ is the matrix of basis functions evaluated at points $\mathbf{t}=(t_{1},\ldots,t_{j})$, the fitted functions are $\hat{Y}_{i}(t)=\mathbf{B}\hat{\mathbf{\mbox{\boldmath$\beta$}}}_{i}$, and the OLS-estimated coefficients are given by $\hat{\mathbf{\mbox{\boldmath$\beta$}}}_{i}=(\mathbf{B}^{\top}\mathbf{B})^{-1}\mathbf{B}^{\top}\mathbf{Y}_{i}$, with $\mathbf{Y}_{i}=(Y_{i}(t_{1}),\ldots,Y_{i}(t_{j}))$. Note that all functions $Y_{i}(t)$ are observed at the same time points $\mathbf{t}=(t_{1},\ldots,t_{j})$.

These principles allow us to formulate Lemma 1. Consider the complete set of fitted coefficients $\mathcal{X}=\{\mathbf{b}_{1},\ldots,\mathbf{b}_{n}\}$ generated by decomposing the functions $(Y_{1}(t),\ldots,Y_{n}(t))$ with a cubic B-spline basis with $K$ interior knots. Define the $j$th subset

$$\mathcal{X}\setminus\mathbf{b}_{j}=\{\mathbf{b}_{1},\ldots,\mathbf{b}_{j-1},\mathbf{b}_{j+1},\ldots,\mathbf{b}_{n}\}$$

as the set of fitted coefficients for the functions $(Y_{1}(t),\ldots,Y_{j-1}(t),Y_{j+1}(t),\ldots,Y_{n}(t))$, i.e., the entire sample with the $j$th function removed. Let $l_{\mathcal{X}}$ denote the complete-data log-likelihood

$$l_{\mathcal{X}}=\sum_{i=1}^{n}\sum_{g=1}^{G}z_{ig}\left[\log\pi_{g}+\log\phi(\mathbf{b}_{i}\mid\bm{\mu}_{g},\bm{\Sigma}_{g})\right],$$

(8)

where the complete-data comprise the coefficients $\{\mathbf{b}_{1},\ldots,\mathbf{b}_{n}\}$ and the cluster memberships $\{\mathbf{z}_{1},\ldots,\mathbf{z}_{n}\}$. Note that $z_{ig}=1$ when coefficient $\mathbf{b}_{i}$ belongs to cluster $g$, and $z_{ig}=0$ otherwise.

Using the density from (9) we can generate the unconditional density of $D$ using a mixture model, i.e.,

$$f(d\mid\bm{\vartheta})=\sum_{g=1}^{G}{\pi}_{g}f_{g}(d\mid\bm{\theta}_{g}),$$

(10)

where $f_{g}(d\mid\bm{\theta}_{g})$ is a shifted and scaled beta density described in (9), and $\bm{\theta}_{g}=\{n_{g},K,\hat{\pi}_{g},\bar{\mathbf{b}}_{g},\mathbf{S}_{g}\}$.

While the distribution in Section 2.2 used the complete-data log-likelihood, we now change to using the traditional log-likelihood, given by

$$\ell_{\mathcal{X}}(\bm{\vartheta})=\sum_{i=1}^{n}\log\left[\sum_{g=1}^{G}\pi_{g}\phi(\mathbf{b}\mid\bm{\mu}_{g},\bm{\Sigma}_{g}+\sigma^{2}(\mathbf{B}^{\top}\mathbf{B})^{-1})\right].$$

(11)

This choice is made because it is outputted by most clustering algorithms. In addition, Clark and McNicholas (2024) show that $\ell_{\mathcal{X}}\rightarrow l_{\mathcal{X}}$ as the clusters separate. The funOCLUST algorithm removes candidate outliers iteratively until the subset log-likelihoods of the B-spline coefficients conform to the derived distribution. It is thus necessary to define a candidate outlier.

The most outlying coefficient is the one whose corresponding log-likelihood is greatest, i.e., the model improved the most in its absence. When treating the coefficients as a multivariate normal dataset, this is consistent with the OCLUST algorithm. However, before proceeding, it is necessary to establish that a function $Y_{i}(t)$ is outlying when its set of coefficients $\mathbf{b}_{i}$ is outlying in terms of multivariate normality.

Relating this back to the taxonomy of Hubert et al. (2015), isolated outliers have $\mathbf{c}=c\mathbf{e}_{j}$, where $\mathbf{e}_{j}$ is a unit vector aligned with the $j$-axis. Shift outliers have $\mathbf{c}=c\mathbf{1}$, where $\mathbf{1}$ is the vector with each entry equal to 1. Amplitude outliers correspond to $a\neq 1,a>0$, with $a<1$ having a dampening effect and $a>1$ having an amplifying effect. Finally, shape outliers correspond to any combination thereof, including where $a<0$.

The funOCLUST algorithm proceeds in two steps: filtering and then clustering. First, the functions $Y_{1}(t),\ldots,Y_{n}(t)$ are filtered by decomposing them into linear combinations of cubic B-splines. In the second step, the funOCLUST algorithm clusters the fitted coefficients, $\mathbf{b}_{1},\ldots,\mathbf{b}_{n}$ and detects outliers. In each iteration of the second step, subsets of the coefficients are clustered and their corresponding log-likelihood recorded. The distribution of these subset log-likelihoods is measured against the theoretical distribution with the KL divergence. The candidate outlier is then removed before the next iteration begins. The model is free of outliers when KL is minimized. Formally, the funOCLUST algorithm is given in Algorithm 1.

This work was supported by an NSERC Canada Graduate Scholarship, an Ontario Graduate Scholarship, an NSERC Discovery Grant, the Canada Research Chairs program, and a Dorothy Killam Fellowship.