Multivariate Functional Data Modeling with Time-varying Clustering

TL;DR

This paper develops a model-based clustering approach for multivariate functional data collected over time at multiple sites, allowing for time-varying clusters using Gaussian processes and Dirichlet processes, demonstrated on ozone and PM10 data from Mexico City.

Contribution

It introduces a novel framework for time-varying clustering of multivariate functional data using Gaussian processes and Dirichlet processes, with a practical partitioning approach for continuous-time data.

Findings

Effective clustering of ozone and PM10 levels across sites.

Demonstrated the method on real environmental data from Mexico City.

Captured temporal changes in pollutant patterns.

Abstract

We consider the situation where multivariate functional data has been collected over time at each of a set of sites. Our illustrative setting is bivariate, monitoring ozone and PM$_{10}$ levels as a function of time over the course of a year at a set of monitoring sites. The data we work with is from 24 monitoring sites in Mexico City which record hourly ozone and PM$_{10}$ levels. We use the data for the year 2017. Hence, we have 48 functions to work with. Our objective is to implement model-based clustering of the functions across the sites. Using our example, such clustering can be considered for ozone and PM$_{10}$ individually or jointly. It may occur differentially for the two pollutants. More importantly for us, we allow that such clustering can vary with time. We model the multivariate functions across sites using a multivariate Gaussian process. With many sites and several functions at each site, we use dimension reduction to provide a stochastic process specification for the distribution of the collection of multivariate functions over the say $n$ sites. Furthermore, to cluster the functions, either individually by component or jointly with all components, we use the Dirichlet process which enables shared labeling of the functions across the sites. Specifically, we cluster functions based on their response to exogenous variables. Though the functions arise in continuous time, clustering in continuous time is extremely computationally demanding and not of practical interest. Therefore, we employ a partitioning of the time scale to capture time-varying clustering.

Figures

Tables

Table 1. Table 1 : Mean number of clusters, averaged over months and iterations

$\alpha =$	$\frac{1}{1000}$	$\frac{1}{10}$	$1$	$10$	$1000$
${\tau }^2\sim IG(1,1)$	5.33	5.42	5.75	5.83	6.39
${\tau }^2\sim IG(101,10)$	5.08	5.31	6.15	6.05	6.32

Equations

$$Y_i^{(k)}(t)=f_i^{(k)}(t)+{ϵ}_i^{(k)}(t)$$

$$f_i^{(k)}(t)={\mu }_i^{(k)}(t)+{\eta }_i^{(k)}(t).$$

$$\textbf{Y}_{i}(t)=\mathbf{B}_{i}\textbf{X}_{i}(t)+\mathbf{G}_{i}\textbf{Z}_{i}(t)+\boldsymbol{\eta}_{i}(t)+\mbox{\boldmath$\epsilon$}_{i}(t).$$

$$\textbf{Y}_{i}^{m}(t)=\mathbf{B}_{i}^{m}\textbf{X}_{i}(t)+\mathbf{G}_{i}^{m}\textbf{Z}_{i}(t)+\boldsymbol{\eta}_{i}(t)+\mbox{\boldmath$\epsilon$}_{i}(t),\hskip 5.69046ptt_{m-1}<t\leq t_{m}.$$

$$\mbox{\boldmath$\eta$}^{(k)}(t)=\left(\begin{array}[]{c}\eta_{1}^{(k)}(t)\\ \eta_{2}^{(k)}(t)\\ .\\ .\\ .\\ \eta_{n}^{(k)}(t)\\ \end{array}\right)=\mbox{\boldmath$\Lambda$}^{(k)}\mbox{\boldmath$\omega$}^{(k)}(t)$$

$$\mbox{\boldmath$\omega$}_{l}(t)=\left(\begin{array}[]{c}\omega_{l}^{(1)}(t)\\ \omega_{l}^{(2)}(t)\\ .\\ .\\ .\\ \omega_{l}^{(K)}(t)\\ \end{array}\right)=\textbf{A}\mbox{\boldmath$\nu$}_{l}(t)=\textbf{A}\left(\begin{array}[]{c}\nu_{l}^{(1)}(t)\\ \nu_{l}^{(2)}(t)\\ .\\ .\\ .\\ \nu_{l}^{(K)}(t)\\ \end{array}\right)$$

$$[{\nu }_l^{(k)}(t_1),{\nu }_l^{(k)}(t_2),...,{\nu }_l^{(k)}(t_{s-1}),{\nu }_l^{(k)}(t_s),...]$$

$$[{\nu }_l^{(k)}(t_1)][{\nu }_l^{(k)}(t_2)|{\nu }_l^{(k)}(t_1)]...[{\nu }_l^{(k)}(t_s)|{\nu }_l^{(k)}(t_{s-1})]...$$

$$\mbox{\boldmath$\Lambda$}^{(k)}=\left({\mbox{\boldmath$\Lambda$}^{(k)}_{1}},...,{\mbox{\boldmath$\Lambda$}^{(k)}_{n}}\right)^{T},$$

$$\text{cov}(\eta_{i}^{k}(t),\eta_{i^{\prime}}^{k^{\prime}}(t^{\prime})=(\mbox{\boldmath$\Lambda$}_{i}^{(k)})^{T}\mbox{\boldmath$\Lambda$}_{i^{\prime}}^{(k^{\prime})}\sum_{j=1}^{K}A_{kj}A_{k^{\prime}j}e^{-\phi_{(j)}(t^{\prime}-t)}.$$

$$p(\text{Y})$$

$$\displaystyle\pi(\mbox{\boldmath$\zeta$})$$

$$\pi ({\zeta }_i|{\zeta }_{-i})=\frac{\alpha }{\alpha +n-1}{G}_0({\zeta }_i)+\frac{1}{a+n-1}\sum _{l<i}{\delta }_{{\zeta }_l}\left({\zeta }_i\right).$$

$$Y_{i}^{(k)}(t)=\textbf{X}_{i}(t)^{T}\mbox{\boldmath$\beta$}_{i,m}^{(k)}+\textbf{Z}_{i}(t)^{T}\mbox{\boldmath$\gamma$}_{i,m}^{(k)}+\eta_{i}^{(k)}(t)+\epsilon_{i}^{(k)}(t),\hskip 5.69046ptt_{m-1}<t\leq t_{m}.$$

$$\pi(\mbox{\boldmath$\beta$}_{im}^{(1)},...,\mbox{\boldmath$\beta$}_{im}^{(K)})=\prod^{n_{s}}_{i=1}\left[\frac{\alpha}{\alpha+i-1}\prod^{K}_{k=1}\mathcal{N}\left(\mbox{\boldmath$\beta$}_{i,m}^{k}|m^{\beta^{(k)}},V^{\beta^{(k)}}\right)+\frac{1}{\alpha+i-1}\sum_{l<i}\prod^{K}_{k=1}\delta_{\mbox{\boldmath$\beta$}_{l,m}^{(k)}}\left(\mbox{\boldmath$\beta$}_{i,m}^{(k)}\right)\right],\\$$

$$\sqrt{Y_i^O(t)}$$

$$\log \left(Y_i^{PM}(t)\right)$$

$$\begin{gathered} \text{cov}({\eta }_i^k(t),{\eta }_{i^{\prime }}^{k^{\prime }}(t^{\prime }))=\text{cov}(\sum _{j=1}^K{A}_{kj}{\varphi }_{ij}^{(k)}(t),\sum _{j=1}^K{A}_{k^{\prime }j}{\varphi }_{i^{\prime }j}^{(k^{\prime })}(t^{\prime }))= \\ \sum _{j=1}^K\sum _{j^{\prime }=1}^K{A}_{kj}{A}_{k^{\prime }{j}^{\prime }}\text{cov}({\varphi }_{ij}^{(k)}(t),{\varphi }_{i^{\prime }{j}^{\prime }}^{(k^{\prime })}(t^{\prime })) \end{gathered}$$

$$\textbf{Y}(t)\sim MVN(\mbox{\boldmath$\mu$}(t),\mbox{\boldmath$\Sigma$}+\textbf{I}\otimes\textbf{D})$$

$$[\textbf{Y}(t^{\prime})|\textbf{Y}(t)]\sim MVN(\mbox{\boldmath$\mu$}(t^{\prime})+\mbox{\boldmath$\Sigma$}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})}(\mbox{\boldmath$\Sigma$}+\textbf{I}\otimes\textbf{D})^{-1}(\textbf{Y}(t)-\mbox{\boldmath$\mu$}(t)),$$

$$(\mbox{\boldmath$\Sigma$}+\textbf{I}\otimes\textbf{D})-\mbox{\boldmath$\Sigma$}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})}(\mbox{\boldmath$\Sigma$}+\textbf{I}\otimes\textbf{D})^{-1}\mbox{\boldmath$\Sigma$}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})}).$$

$$\textbf{Y}^{m}_{i}(t)=\textbf{G}^{m}_{i}\textbf{Z}_{i}(t)+\textbf{B}^{m}_{i}\textbf{X}_{i}(t)+\mbox{\boldmath$\eta$}_{i}(t)+\mbox{\boldmath$\epsilon$}_{i}(t),$$

$$\displaystyle\mbox{\boldmath$\gamma$}_{im}^{(k)}|\cdots$$

$$\displaystyle\mbox{\boldmath$\beta$}_{cm}^{(k)}|\cdots$$

$${\Lambda }_i^{(k)}$$

$$A_{kl}$$

$${\tau }_k^2$$

$${\xi }_{im}$$

$${\nu }^{(k)}(t)$$

$${\Sigma }_{{\gamma }_{im}^{(k)}}$$

$${\mu }_{{\gamma }_{im}^{(k)}}$$

$${\Sigma }_{{\beta }_{cm}^{(k)}}$$

$${\mu }_{{\beta }_{cm}^{(k)}}$$

$${\Sigma }_{{\Lambda }_i^{(k)}}$$

$${\mu }_{{\Lambda }_i^{(k)}}$$

$$v_{A_{kl}}$$

$${\mu }_{A_{kl}}$$

$$a_{{\tau }_k^2}^{\ast }$$

Full text

We model the multivariate functions across sites using a multivariate Gaussian process. With many sites and several functions at each site, we use dimension reduction to provide a stochastic process specification for the distribution of the collection of multivariate functions over the say $n$ sites. Furthermore, to cluster the functions, either individually by component or jointly with all components, we use the Dirichlet process which enables shared labeling of the functions across the sites. Specifically, we cluster functions based on their response to exogenous variables. Though the functions arise in continuous time, clustering in continuous time is extremely computationally demanding and not of practical interest. Therefore, we employ a partitioning of the time scale to capture time-varying clustering.

The data we work with is from 24 monitoring sites in Mexico City which record hourly ozone and PM10 levels. We use the data for the year 2017. Hence, we have 48 functions to work with. We provide a Gaussian process model for each function using continuous time meteorological variables as regressors along with adjustment for daily periodicity. We interpret similarity of functions in terms of their shape, captured through site-specific coefficients and use these coefficients to develop the clustering.

Keywords: Dimension reduction; Dirichlet process; hierarchical model; latent factor models; multivariate Gaussian process; ozone; PM10

1 Introduction

The beginnings of functional data analysis date to Grenander, (1950) and Rao, (1958), who applied the theory of stochastic processes to curve data. The term “functional data analysis” came into common usage with Ramsay, (1982) and Ramsay and Dalzell, (1991). Since then the field has undergone rapid growth, spurred by the re-framing of existing problems in the context of functional data analysis along with scientific advances giving rise to larger amounts and new types of functional data. Early work focused mainly on curves as the object of study. Ramsay, (2005) and Ramsay and Silverman, (2007) provide a thorough introduction to the basic concepts of functional data analysis. Applications have been found in MRI brain images (and, in fact, for imaging in general), finance, climatic variation, spectrometry data, and time-course gene expression data, to name a few. For a comprehensive overview of applications, see Ullah and Finch, (2013).

Functional data is viewed explicitly as living on $\mathbb{R}^{d}$, where $d=1,2,3$ in most applications. The response is indexed by arguments that are viewed as continuous. Of course, for any function, samples are only ever collected at a finite set of arguments. In fact, finite-dimensional multivariate approaches were used originally for FDA.

Explicit modeling of functional data is usually carried out in one of two ways: (i) as linear combinations of some set of basis functions, or (ii) as realizations of some stochastic process. In the former case, weights are chosen for the basis functions so that the fit to the data is optimized (Ramsay and Silverman, (2007); Morris, (2015) ). For instance, splines implement this by imposing a fixed partition on the domain of the function and allowing for the weights to vary across the partition elements. Observed data are then viewed as noisy realizations of the underlying splines, with implementation requiring finite truncation of the basis function expansion. To address the high dimensionality of functional data, we find extensions of principal component analysis to the functional domain. Formally, functional principal components analysis (FPCA) specifies each function as an infinite linear combination of orthonormal basis functions using the Karhunen-Loéve Theorem. Practically, a finite linear combination is adopted. Arguably, FPCA is among the most popular techniques for analyzing functional data.

The alternative approach which we pursue here is to view the data for each individual as noisy observations from a realization of a stochastic process. That is, the true functions are realizations of a stochastic process over the subset of $\mathbb{R}^{d}$ of interest. Gaussian processes are often a convenient choice since they require only the specification of a mean function and a valid covariance function. They are flexible and provide convenient standard multivariate normal distribution theory as well as straightforward interpolation. Extensions of this approach, based upon Dirichlet process modeling ideas, developed below, enable our functional modeling to exhibit clustering which is the primary contribution of this paper.

With multivariate functional data, methods employed for clustering of multivariate data vectors can be adapted to the task of clustering functional data, with the usual goal being to group observations such that within-cluster distance under some criterion is minimized. For example, K-means clustering is a widely applied algorithm for this problem (Hartigan and Wong,, 1979; Abraham et al.,, 2003). Here, we avoid an algorithmic approach and, instead, pursue model-based clustering which can be placed under the heading of mixture modeling. At a “big picture” level, the different mixture components act as the cluster centers, and optimization is equivalent to the model fitting process (Sugar and James,, 2003; Jacques and Preda,, 2014). We note that this approach addresses clustering for univariate functional data; see Wang et al., (2016) for a comprehensive review.

Multivariate functional data typically comprise several simultaneously recorded time course measurements for a sample of subjects or experimental units. They are viewed as realizations sampled from multivariate random functions. Analysis of multivariate functional data has received relatively less attention than univariate functional data methods though chemometrics is an area where multivariate FDA has been discussed (see, e.g., Wang and Chen,, 2015).. Again, functional principal component analysis (FPCA) serves as a fundamental tool. Methods for analyzing multivariate functional data that utilize univariate FPCA include dynamical correlation analysis for multivariate longitudinal observations (Dubin and Müller,, 2005), modeling the relationship of paired longitudinal observations (Zhou et al.,, 2008), and regularized FPCA for multidimensional functional data (Kayano and Konishi,, 2009). When the components of multivariate functional data are measured in the same units and have similar variation, classical multivariate FPCA that concatenates the multiple functions into one to perform univariate FPCA can work well.

The approach we adopt here employs multivariate Gaussian process specification for the functions. Shi and Choi, (2011) presents discussion in this vein. In particular, we consider the situation where multivariate functional data has been collected over time (so $d=1$) at each of a set of sites. Our illustrative setting is bivariate, monitoring ozone and PM10 levels as a function of time over the course of a year, at a set of monitoring sites. Our objective is to implement model-based clustering of the functions across the sites. Using our example, such clustering can be considered for ozone and PM10 individually or jointly. It may occur differentially for the two pollutants. More importantly for us, we allow that such clustering can vary with time.

Specifically, we consider the challenge of modeling collections of multivariate functional data incorporating their relationships to exogenous variables. We consider the setting where we have multiple functional outcomes that are related by regression to predictors that are themselves functional, up to pure error. In addition, mean [math] random effects are introduced using Gaussian processes, one for each function at each site. These random effects provide local (in time) adjustment to the regression functions.

In this setting, we interpret clustering to signify a common functional form for the functions, ignoring the random effects. If our functions are specified as Gaussian processes with a mean that is a regression in exogenous variables, then we cluster two functions if they share the same regression coefficients. In this regard, we can distinguish between the “shape” of the function and the “level” of the function. In our simple setting, sharing shape is captured through common slope coefficients; sharing level is captured by additionally sharing the intercept. Most importantly, we don’t cluster functions in terms of response. That is, the various functions may be observed at quite different levels of the exogenous variables resulting in quite different observed response levels. Again, we seek to cluster functions which share shape.

We are fortunate that all of our functions run for the same length and period of time. So, we can avoid some of the well-discussed challenges in the literature regarding registration and time warping in order to “align” the functions (see Gervini, (2014) for a general discussion and see Telesca and Inoue, (2008) within the Bayesian framework).

In our context, we have $n=24$ sites with two functions at each site. Each function is observed hourly over a year, resulting in $8760$ observations for each. We seek to capture dependence between the functions within a site as well as dependence between the functions across the sites. We use dimension reduction to specify the joint distribution of the collection of multivariate functions over the $24$ sites, expressing the $2\times 24$ functions in terms of $2\times r$ functions, where, illustratively, $r$ might be $3$ or $5$. We model these $2r$ functions using multivariate Gaussian processes which induce a Gaussian process specification for each function at each site.

With regard to our clustering objective, it is clear that two Gaussian process realizations over a subset say $(0,T]$ of $R^{1}$ will agree with probability [math]. So, in order to achieve model-based clustering, we imagine a finite or countably infinite set of functional realizations over $(0,T]$. Then, two functions are clustered if they share a common realization from this set. That is, we need to add a further level of hierarchical modeling to provide the set of functional realizations which are available to share. A convenient way to do this is through the Dirichlet process (which is discussed and referenced below). Such processes provide random discrete distributions, i.e., a random set of atoms with random probabilities assigned to the atoms. In our setting the random atoms are functions. When functional data is modeled in this fashion, two functions are clustered if they share the same atom, i.e., each atom is assigned a label and the two functions receive the same label.

In fact, there is more subtlety here according to whether or not the functions are explained through the introduction of regressors. If no regressors are considered, then all of the atoms are realizations of a constant mean Gaussian process over $(0,T]$ whence all of the functions we are modeling are constant mean Gaussian processes. If regressors are introduced into the means of the Gaussian processes, then we can consider two functions as clustered if they share the same mean function, ignoring the mean [math] residual (local adjustment) process. Specifically, with a regression function of the form say \textbf{X}(t)\mbox{\boldmath\beta}_{i}^{(k)} for the $k$th function at site $i$, clustering involves sharing of coefficient vectors for the \mbox{\boldmath\beta}_{i}^{(k)}. Now, the Dirichlet process would introduce atoms which are vectors of regression coefficients. We adopt this latter approach here since we seek to understand the nature of differential pollutant response to meteorological predictors across the sites.

We add a further wrinkle. We imagine that the relationship between response and predictors is changing over time. In particular, we use a partitioning of the time scale to capture time-varying change in these response/predictor relationships. In turn, this leads to an investigation of time-varying clustering of the functions. For instance, what can we say about how clustering of the sites arises say in summer vs. in winter? Though the functions arise in continuous time, continuous time clustering is not of practical interest and, additionally, is extremely computationally demanding. Instead, illustratively, we work at monthly scale.

Finally, though the sites are spatially referenced, we intentionally do not model them spatially. We are not interested in interpolation of the functions. Though there may be a conceptual realization of a multivariate function at every spatial location, we don’t imagine clustering functions which we don’t observe. In fact, we seek to cluster functions through their functional similarity which is not a spatial notion. As a post model fitting exercise, we can assess how strong clustering is in terms of inter-site distance.

The data we work with is from 24 monitoring sites in Mexico City which record hourly ozone and PM10 levels. We use the data for the year 2017. We have 48 functions to work with and we provide a Gaussian process model for each function using continuous time meteorological variables - relative humidity and temperature - as regressors along with daily periodicity to supply the mean of each of the Gaussian processes. The model fitting provides site-specific coefficients and we use these coefficients to develop the clustering. Again, since the data at each site is collected over the same window of time, no alignment techniques are needed. Furthermore, the fact that these stations measure the same quantities and are close in proximity suggests that the data can be effectively modeled using a representation of lower rank than the 24 pairs of curves we observe.

The format of the paper is as follows. We continue with a review of the modeling components we utilize in Section 2, focusing on multivariate Gaussian processes (Section 2.1) and Dirichlet processes (Section 2.2). We present our modeling framework in Section 3, including a discussion of prior distributions and model fitting. In Section 4, we present our analysis of ozone and $\text{PM}_{10}$ in Mexico City in year 2017. We discuss data attributes in Section 4.1, specific modeling details in 4.2, and results of our analysis in Section 4.3. Lastly, in Section 5, we summarize our approach and discuss possible extensions.

2 The Model Components

The general form of the modeling is as follows. At site $i$ for function $k$ at time $t$, we assume

[TABLE]

where $Y$ denotes the observed response, $f$ denotes the “true” function and $\epsilon$ denotes pure error. Further, we express the function in terms of a “fixed” effects term, $\mu_{i}^{(k)}(t)$ and a “random” effects term, $\eta_{i}^{(k)}(t)$, i.e.,

[TABLE]

Here, $\mu_{i}^{(k)}(t)$ is a regression term which, illustratively, we specify as \textbf{X}_{i}(t)^{T}\mbox{\boldmath\beta}_{i}^{(k)}+\textbf{Z}_{i}(t)^{T}\mbox{\boldmath\gamma}_{i}^{(k)}, i.e., a component and site specific vector of regression coefficients. Here, we let \mbox{\boldmath\beta}_{i}^{(k)} be the regression coefficients that are involved in clustering and \mbox{\boldmath\gamma}_{i}^{(k)} are the regression coefficients not involved in clustering allowing for settings where clustering may be of interest for only a subset of covariates. $\eta_{i}^{(k)}(t)$ is a mean 0 Gaussian process and these processes are dependent across $i$ and across $k$. We can assemble the model at site-level to

[TABLE]

Here, with $n$ sites, $i=1,2,...,n$. With $K$ functions at each site and $p$ predictors (including intercept), $\mathbf{B}_{i}$ and $\mathbf{G}_{i}$ are $K\times p$ with $k$th row, (\mbox{\boldmath\beta}_{i}^{(k)})^{T} and (\mbox{\boldmath\gamma}_{i}^{(k)})^{T}, respectively. $\boldsymbol{\eta}_{i}(t)$ is a $K$-dimensional Gaussian process and the $\boldsymbol{\eta}_{i}(t)$ are dependent across $i$. So, in total, we have $nK$ functions specified as an $nK$ dimensional Gaussian process. Finally, the \mbox{\boldmath\epsilon}_{i}(t) are pure error vectors, independent across components, and normally distributed with mean [math] and $k$th component having variance $\tau^{2(k)}$.

Equation (3) looks like a familiar mixed effects model for continuous time data. However, we introduce the following novelties. First, we partition the time window into say $M$ segments, i.e., $(0,t_{1},t_{2},...,t_{M-1},T]$. On the $m$th segment we revise (3) to

[TABLE]

As a result, we have time-varying coefficients. However, as residuals, we assume that the random effects Gaussian processes and the pure error processes are not partition-dependent.

Second, we implement a dimension reduction strategy to model the $\boldsymbol{\eta}_{i}(t)$ using an $rK$ dimensional Gaussian process where $r$ is much smaller than $n$, typically, say $3$ to $5$. Since the number of temporal observations for each function is large (hourly over an entire year in our case), even with independent Gaussian process, we will wind up working with very high dimensional covariance matrices for each process. For likelihood evaluation, these matrices require inverse and determinant calculation. There are strategies for handling such dimensionality challenges such as the predictive process (Banerjee et al.,, 2008) or the nearest neighbor Gaussian process (Datta et al.,, 2016). However, since our processes are over time, a convenient way to avoid this computational challenge is to adopt an exponential covariance function. This corresponds to a continuous time AR(1) process (see, e.g., Brockwell et al.,, 2007). With equally spaced time points, this enables writing of the joint distribution sequentially over time without resorting to calculating the joint covariance matrix.

Third, we recall our objective of clustering. With time-varying coefficients for each component function, we can investigate clustering of the components across the sites within each partition. We can investigate whether and, if so, how clustering behavior changes with partition, i.e., with time. As discussed in the Introduction, the Dirichlet process provides a convenient distributional mechanism to enable ties, i.e., to enable, within partition $m$, \mbox{\boldmath\beta}_{i}^{m(k)}=\mbox{\boldmath\beta}_{i^{\prime}}^{m(k)} for one or more $k$’s.

2.1 Multivariate Gaussian Process Factor Models

The multivariate Gaussian process (GP) supplies both a data model for multivariate data collected at continuously indexed sites as well as a multivariate random effects process model when dependence between the processes is sought. For examples of the first setting see, e.g., Wackernagel, (1994); Gneiting et al., (2010); Vandenberg-Rodes and Shahbaba, (2015). For the second, see, e.g., Banerjee et al., (2014), Chap 9. Regardless, specification of a multivariate Gaussian process requires a valid cross-covariance function. In the above notation, indexing over time, we seek a site specific $K\times K$ matrix valued function, $C_{i}(t,t^{\prime})$ with $(C_{i}(t,t^{\prime})_{k,k^{\prime}}=\text{cov}(\eta_{ik}(t),\eta_{ik^{\prime}}(t^{\prime}))$. There are numerous constructions for $C$ including coregionalization, kernel convolution, and convolving covariance functions (Banerjee et al.,, 2014, Chap 9). Direct development has been pursued for the Matérn covariance class in, e.g., Gneiting et al., (2010); Apanasovich et al., (2012). Development through stochastic partial differential equations is presented in, e.g., Vandenberg-Rodes and Shahbaba, (2015); Singh et al., (2018). For convenience, particularly since, in our application $K=2$, we employ coregionalization. However, any of the foregoing choices could be used equally well.

As above, we seek to have the \mbox{\boldmath\eta}_{i}(t) dependent across $i$. We want to specify an $nK$ dimensional GP through a dimension reduction which requires only an $rK$ dimensional GP. That is, for each $i$, we want to specify \eta_{i}^{(k)}(t)=\Lambda_{i}^{(k)T}\mbox{\boldmath\omega}^{(k)}(t) where $\Lambda_{i}^{(k)}$ is an $r\times 1$ vector and \mbox{\boldmath\omega}^{(k)}(t) is an $r$-dimensional GP with mean $\mathbf{0}$ and correlation matrix R. Assembling across $n$, we obtain the $n\times 1$ vector, \mbox{\boldmath\eta}^{(k)}(t) such that

[TABLE]

where \mbox{\boldmath\Lambda}^{(k)} is an $n\times r$ matrix. What we have here is a latent factor analysis (Lopes and West,, 2004) where \mbox{\boldmath\Lambda}^{(k)} provides the matrix of factor loadings for the $k$th component of the \mbox{\boldmath\eta}_{i}(t) vectors. We use correlation functions (covariance functions with unit variance) in the \mbox{\boldmath\omega}^{(k)}(t) above in order to identify the factor loadings. We emphasize that the dimension reduction is across sites, not across components of \mbox{\boldmath\eta}_{i}(t). Specifically, we provide this dimension reduction without introducing the spatial locations of the sites. In fact, \mbox{\boldmath\Lambda}^{(k)}\mbox{\boldmath\Lambda}^{(k)T} provides the $n\times n$ covariance matrix between the sites. In order to see how dependence for the $k$th component associates with distance between sites, we will investigate the (\mbox{\boldmath\Lambda}^{(k)}\mbox{\boldmath\Lambda}^{(k)T})_{i,i^{\prime}} vs. the corresponding inter-site distances, $d_{i,i^{\prime}}$.

To complete the specification we need to provide the modeling for the set of $rK$ processes, $\omega_{j}^{(k)}(t)$, $j=1,2,...,r;k=1,2,...,K$. Here, we bring in coregionalization, writing the $K\times 1$ vector, \mbox{\boldmath\omega}_{l}(t) as

[TABLE]

where A is a $K\times K$ coregionalization matrix and the $\nu_{l}^{(k)}(t)$ are i.i.d mean [math], variance $1$ GPs. In particular, by adopting an exponential correlation function, we assume they are continuous time first order Markov processes. Hence, the joint distribution across time can be peeled off sequentially, i.e.,

[TABLE]

becomes

[TABLE]

where $[\nu_{l}^{(k)}(t_{s})|\nu_{l}^{(k)}(t_{s-1})]=N(e^{-\phi_{(k)}(t_{s}-t_{s-1})}\nu_{l}^{(k)}(t_{s-1}),1-e^{-2\phi_{(k)}(t_{s}-t_{s-1})})$ with $[\nu_{l}^{(k)}(0)]=N(0,1)$.

Again, the factor loadings for the $k^{\text{th}}$ component are captured in the $n\times r$ matrix

[TABLE]

which brings in the dependence across the sites. The coregionalization matrix A induces dependence between the random effects at a given site. Since coregionalization is applied to random effects, for convenience we can take A to be lower triangular.

We recall that the factor model is not identifiable and thus constraints to $\Lambda$ are applied to remedy this. It is common to regain identifiability of $\Lambda$ by either (i) constraining $\Lambda$ to be orthogonal (e.g. Seber,, 2009) or lower triangular (see Geweke and Zhou,, 1996; Aguilar and West,, 2000; Lopes and West,, 2004) or (ii) using prior distributions that induce sparseness or that provide shrinkage (see Bernardo et al.,, 2003; Bhattacharya and Dunson,, 2011, for example). Identifiability can also be recovered by constraining parameters of the GPs. For example, the range parameters of the GPs could be ordered to make the factor loading identifiable. In this case, the factors would be interpretable in the way they each account for different ranges of autocorrelation. We discuss this approach in our model fitting discussion in Section 3.1 and in our analysis of hourly pollution levels presented in Section 4.

We conclude this subsection by noting that our specification enables convenient separation of the factor loadings and the coregionalization entries in the dependence structure for the $\eta_{i}^{k}(t)$. Specifically, for say $t^{\prime}\geq t$, we have

[TABLE]

We see the separation of the factor loadings at $i$ and $i^{\prime}$ for $k$ and $k^{\prime}$ from the coregionalization at $k$ and $k^{\prime}$. The proof of result (9) is provided in the Appendix. Following from this, we can directly obtain all of the marginal and conditional distribution theory needed for the model fitting. These results are also supplied in the Appendix.

2.2 Dirichlet Processes

Bayesian hierarchical specifications provide a convenient way to develop mixture models. Specifically, the Dirichlet process has been used to avoid explicitly specifying the number of clusters present in the data. That is, we need not struggle with the identifiability and slow-mixing challenges of finite mixture models. The Dirichlet process provides a countable mixture model where, again, the intent is clustering rather than learning about the number of mixture components and the parameters associated with each component. In this regard, the Dirichlet process model serves as a prior assigning random distributions which have countably infinite support; it dates at least to (Ferguson, 1973). The stick-breaking representation presented by (Sethuraman, 1994) consequentially expanded its applicability. Model fitting using Markov chain Monte Carlo (MCMC) was initially presented in (Escobar and West, 1995). In particular, Dirichlet process mixtures of stochastic processes are possible, allowing for extensions to accommodate functional data.

The random discrete distributions can place mass on any random objects of interest. Hence, they can adopt as functional atoms, realizations of Gaussian processes or realizations of spline-specified functions (see, e.g., MacEachern,, 2000; Gelfand et al.,, 2005; Duan et al.,, 2007; Petrone et al.,, 2009). However, as we noted in the Introduction, here, clustering corresponds to sharing regression coefficient vectors since, across sites, we seek to cluster pollutant response to meteorological predictors.

Following Ferguson, (1973), Dirichlet processes are defined in terms of a base distribution $G_{0}$ along with a concentration parameter $\alpha$; we say that a realization $G$ from a Dirichlet process is a random, almost surely discrete distribution that is close, in some sense, to the base measure $G_{0}$ with the closeness determined by the concentration parameter $\alpha$ and we write $G\sim DP(\alpha,G_{0})$.

The constructive representation introduced by Sethuraman (1994), showed that countably infinite mixtures of atoms with suitable stick-breaking weights are equivalent to draws from a Dirichlet process. The base measure of the Dirichlet process is then the distribution that provides the mixture components in the infinite mixture. That is, $G\sim DP(\alpha,G_{0})\Longleftrightarrow G=\sum_{l=1}^{\infty}p_{l}\delta_{\theta_{l}^{*}}$ where the $\theta_{l}^{*}$ are i.i.d. realizations from $G_{0}$ and the $p_{l}$ arise as stick-breaking probabilities, i.e., $p_{1}=w_{1}$ and $p_{l}=w_{l}\Pi_{j=1}^{l-1}(1-w_{j})$ with the $w_{l}\sim\text{Beta}(1,\alpha)$.

The stick-breaking representation provides intuition for sampling schemes. The mixture representation suggests introducing “latent classes” that individual observations are assigned to, resulting in distributions that are easier to sample. This enables simplified Markov chain Monte Carlo (MCMC) model fitting. In practice it is customary to truncate the number of $p_{l}$’s such that $\sum_{l=1}^{L}p_{l}$ is suitably close to $1$ or, equivalently, $\Pi_{l=1}^{L}(1-w_{l})$ is sufficiently close to [math]. Alternatively, the Dirichlet process can be represented be an urn scheme (Blackwell et al.,, 1973). The Gibbs sampler for this approach is straightforward and presented in Escobar, (1994); Escobar and West, (1995). Using the Pólya-urn representation of the Dirichlet process, the model for data $\textbf{Y}=\left(y_{1},...,y_{n}\right)$ is represented

[TABLE]

This representation reduces the infinite mixture to (at most) an $n$-dimensional mixture, where \mbox{\boldmath\zeta}=(\zeta_{1},...,\zeta_{n})^{T} assigns membership to specific mixture components. The full conditional for each $\zeta_{i}|\zeta_{-i}$ is

[TABLE]

The representations in (10) and (11) indicate that with more $\zeta_{i}$ assigned to the same mixture component, it becomes more likely that other observations will be assigned to that component.

In our setting, we use a DP mixture prior distribution on function-specific regression coefficients to identify unique responses to exogenous variables or predictors. Heuristically, this prior distribution groups functional data that have similar responses to predictors by inducing clusters or ties between regression coefficients. The DP-mixing accrues to the functions $Y_{i}^{(k)}(t)$’s. That is, each component at each site has a random coefficient vector, \mbox{\boldmath\beta}_{i}^{(k)}.

Using the Dirichlet process, these \mbox{\boldmath\beta}_{i}^{(k)} can be modeled jointly across $k$ or individually for each $k$ according to whether we seek to cluster by site or by components within site. In either case, under the Dirichlet process, each \mbox{\boldmath\beta}_{i}^{(k)} receives a random label. For joint clustering, we denote it as $L_{i}$ and if $L_{i}=l$ then we set $B_{i}=B_{l*}$, the $l$th atom in the distribution of $G$. If it is by component, then we denote it as $L_{i}^{(k)}$ and if $L_{i}^{(k)}=l$, \mbox{\boldmath\beta}_{i}^{(k)}=\mbox{\boldmath\beta}_{l}^{(k)*}, now, the $l$th atom in the distribution say, $G^{(k)}$ for the $k$th component. The DP-mixing occurs because Y_{i}^{(k)}(t)|\mbox{\boldmath\beta}_{i}^{(k)} is normally distributed, i.e., the $Y_{i}^{(k)}(t)$ arise from a mixture distribution given the \mbox{\boldmath\beta}_{l}^{(k)*} atoms. Articulating the specification in this fashion also reveals that we primarily care about when two sites share labels with little interest in the values of the associated atoms. Atoms will change across iterations of MCMC model fitting and so will labels.

Lastly, we allow this labeling to vary over time using exchangeable DP mixtures as prior distributions over coarse time scales. This requires introducing a superscript $m$ in each of the quantities in the previous paragraph. In some settings, incorporating dependence between the base measures for the time-varying DP mixtures may be warranted (e.g., the mean of the base measure may depend on the mean over the previous time period). Fitting such models is very challenging and is beyond the scope of our work here.

3 Methods and Models

We now return to the model in (4) which allows for partitioning of the time scale. Our data is still denoted by $Y_{i}^{(k)}(t)$. This model, at the component scale, becomes

[TABLE]

As a result, the only modification to the foregoing development is that we need to implement a Dirichlet process for the coefficients for each segment $m$ in order to obtain clustering for each partition. We assume that the Dirichlet processes are exchangeable across partitions, i.e., they share a common base measure $G$ and concentration parameter $\alpha$.

3.1 Prior Distributions and Model Fitting

Here, we present model fitting for specifications presented in (4)-(8). In our example below, we cluster on the coefficients associated with the exogenous variables. That is, we seek to cluster with regard to the functional response to environmental variables. So, we do not cluster on the variables capturing daily periodic behavior. As in (3), we include an additional term \textbf{z}_{i}(t)\mbox{\boldmath\gamma}_{i,m}^{(k)} in the mean to capture coefficients for which we are not interested in clustering. These coefficients are modeled with customary vague normal priors so there are no ties in the \mbox{\boldmath\gamma}_{im}^{(k)}).

Model fitting requires prior distributions for \mbox{\boldmath\gamma}_{im}^{(k)}, \mbox{\boldmath\beta}_{im}^{(k)}, $A_{kl}$, $\Lambda_{i}^{(k)}$, $\tau^{2(k)}$, covariance parameters ($\phi_{j}^{(k)}$) for the exponential covariance function), as well the joint DP for \mbox{\boldmath\beta}^{(k)}_{im}. Multivariate normal prior distributions for \mbox{\boldmath\gamma}_{im}^{(k)} are adopted, as is customary because of amenability to Gibbs sampling.

Components of the coregionalization matrix A and \mbox{\boldmath\Lambda}^{(k)} are not uniquely identifiable without constraints. To remedy this, we constrain all components on the diagonal of A to be 1 and all components above the diagonal to be 0. Then, lower diagonal elements are assumed a priori to be Gaussian with high variance. Prior specifications for $\Lambda_{i}^{(k)}$ depends on the approach used (e.g. Lopes and West,, 2004; Bhattacharya and Dunson,, 2011). Our approach relies on ordering $\phi_{j}^{(k)}$ for each GP factor indexed by $j$. Thus, the factor loadings represent how much each site loads onto GPs with slower or faster temporal decays in autocorrelation. With this approach, each $\Lambda_{i}^{(k)}$ is given a weakly informative multivariate normal prior distribution. Gamma, log-normal, or uniform prior distributions for $\phi_{j}^{(k)}$ are common. Alternatively, it is common to fix $\phi_{j}^{(k)}$ for identifiability reasons (Zhang,, 2004) using maximum likelihood or variogram fitting. In fact, we fixed the $\phi$’s (see below) and conducted some sensitivity analysis. Lastly, we use inverse-gamma prior distributions for $\tau^{2}_{k}$.

As in Section 2.2, we specify \mbox{\boldmath\beta}_{im}^{(k)} through a joint Dirichlet process, where the clustering varies month-to-month. For each time window $m$, the joint prior distribution for \mbox{\boldmath\beta}_{im}^{(1)},...,\mbox{\boldmath\beta}_{im}^{(K)} is

[TABLE]

where $\delta_{a}(x)$ is a Dirac-delta function of $x$ that has a point mass at $a$. This prior distribution assumes that the regression coefficients for all variables cluster or group together. This assumption could be relaxed to have independent Dirichlet process prior distributions for the regression coefficients for each outcome. To complete this specification, the concentration parameter $\alpha$ could have a prior distribution with a strictly positive support or could be fixed to control clustering behavior. We conducted sensitivity analysis using fixed values of $\alpha$. Using the prior distributions given here, closed-form full conditional distributions are available for all parameters except $\phi_{j}^{k}$. These are presented in Appendix B.

4 Analysis of Hourly Pollution Levels in Mexico City

4.1 Data

Here, we consider hourly ozone and $\text{PM}_{10}$ measurements from 24 stations in Mexico City during 2017. The locations of the 24 stations and their abbreviated names are plotted in Figure 1. The data are available under Datos/Horarios/Contaminante at http://www.aire.cdmx.gob.mx/default.php. In total, each station has 8760 ozone and $\text{PM}_{10}$ measurements.111Missing hourly measurements were imputed using the corresponding measurements at the nearest station within the same region. If no stations in that region recorded a measurement at that time, then the nearest station in a different region provided the missing value. This was done prior to our receiving the data for analysis.

For both ozone and $\text{PM}_{10}$ , we examine patterns in the data that motivate modeling decisions. For each station, we plot the average of each day, averaged over all hours in the day, and the hourly average, averaged over all hours of the day. In addition, we examine how the relationships between covariates (relative humidity and temperature) and outcomes (ozone and $\text{PM}_{10}$ ) vary month-by-month for each station. To explore this, we fit a linear model for both outcomes at each station for each month using relative humidity and temperature as covariates. Exploring these patterns allows us to assess how relationships between covariates and outcomes differ station-to-station and how they vary over 2017. These patterns for ozone are given in Figure 2 and for $\text{PM}_{10}$ in Figure 3.

Figure 2 shows significant variability in ozone levels over the year and as a function of time-of-day. Ozone levels peak in the Spring (April and May are the peak ozone season). As a function of the time-of-day, ozone levels are highest in the early afternoon, dip in the evening, and reach a minimum in the morning. While these trends are common across stations, there is variability in ozone levels and patterns between stations. In addition to the variability in ozone levels, the results from the monthly regression models show significant variability in the intercept and regression coefficients for humidity and temperature month-to-month, suggesting potential benefit of time-varying effects. However, regression coefficients are negative for relative humidity and positive for temperature.

The annual trends in $\text{PM}_{10}$ are less clear than they are for ozone; however, the highest levels in $\text{PM}_{10}$ occur in the winter. The hourly averages of $\text{PM}_{10}$ in Figure 3 exhibit clear peaks around 8 AM and 7 or 8 PM. Like the regression coefficients for ozone, the effects of humidity and temperature on $\text{PM}_{10}$ vary greatly over 2017 and station-to-station. Again, these summaries indicate that a hierarchical model with dynamic regression coefficients may be appropriate.

4.2 Model

For the Mexico City data, let $Y^{O}_{i}(t)$ be ozone levels in parts per billion (ppb), $Y^{PM}_{i}(t)$ be $\text{PM}_{10}$ concentration in $kg/m^{3}$ at location $i$ and hour $t$. In addition, let $\textbf{x}_{i}(t)$ be temperature in degrees Celsius at station $i$ and hour $t$. We model ozone on the square root scale, following Gelfand et al., (2003); Sahu et al., (2007); Berrocal et al., (2010); White et al., (2019); White and Porcu, (2019). We model $\text{PM}_{10}$ on the log scale as in Cocchi et al., (2007); Huang et al., (2018); White et al., (2019). In addition, we center and scale the transformed outcomes so that the likelihoods for each outcomes are similar. In both cases, these transformations stabilize variance between model mean and variance of residuals and improve predictive performance. We induce clustering hierarchically through an exchangeable monthly Dirichlet process on the regression coefficients. We index months by $m$ such that the models in (14) and (15) apply to $t\in m$.

Explicitly, we have

[TABLE]

where $\epsilon_{i}^{O}(t)$ and $\epsilon_{i}^{PM}(t)$ are assumed to be Gaussian with a mean of 0 and variance of $\tau^{2}_{1}$ and $\tau^{2}_{2}$, respectively. Here, we let coefficients (\mbox{\boldmath\gamma}^{O}_{im} and \mbox{\boldmath\gamma}^{PM}_{im}) for $\textbf{z}(t)^{T}$ (an intercept and sine and cosine terms with periods of 24 hours) be at monthly scale. As noted above, we do not include these parameters as a part of the Dirichlet process because the goal is clustering stations on with regard to response to meteorological variables (humidity and temperature). We fix the number of factors to $r=3$. For each factor, we fix the decay parameters $\phi^{(k)}_{j}$ such that $\phi^{(k)}_{1}=1/24$, $\phi^{(k)}_{2}=1/3$, and $\phi^{(k)}_{3}=1$ for all $k$. The first factor corresponds to long-range autocorrelation (with an effective range of 72 hours), the second factor captures middle-range autocorrelation (an effective range of nine hours), and the third factor accounts for short-term autocorrelation (an effective range of three hours).

We fit our model using the Gibbs sampler presented in Section 3 and Appendix B. Specifically, we ran the model for 150,000 iterations, discarding the first 50,000 samples. Then, we thinned the remaining 100,000 sample to 10,000 samples, keeping every 10th sample.

In Appendix C, we provide a discussion of the sensitivity of the model. To summarize this discussion, the number of clusters in the model is somewhat sensitive to selection of $\alpha$, specifications of the base measure of the Dirichlet process prior, and prior distributions on $\tau^{2}_{k}$. In addition, we found that joint clustering was driven exclusively by ozone. If we used independent Dirichlet processes for ozone and $\text{PM}_{10}$ , then clustering was limited to ozone.

4.3 Results

Our primary goal is to explore how the clustering of stations evolves over time. However, we first address possible issues with regard to our modeling choices. To examine the benefit of modeling ozone and $\text{PM}_{10}$ jointly, we note the significance of the coregionalization parameter $a_{12}$; it has a posterior mean of $0.158$ with 95% credible interval $(0.141,0.180)$. To address not including space in our model, because $\Lambda^{(k)}{\Lambda^{(k)}}^{T}$ represents the between-site covariance in our model, we plot the components of these matrices against the distance between sites in Figure 4. We can assess whether the factor loadings are capturing spatial behavior. If there was consequential spatial dependence, we would expect a negative relationship between the between-site covariance and distance. However, the plots in Figure 4 do not suggest any strong spatial patterns. We also plot the proportion of times that sites share a label, averaged over months, against the distance between sites in Figure 5. Again, there does not appear to be a very strong relationship between clustering and space.

Although Figures 4 and 5 suggest that the data do not have a strong spatial component after accounting for temporal autocorrelation, looking at clustering for individual months, we find some spatial patterns. For illustration, we provide the station location and cluster label (indexed by color) for May and December in Figure 6. In May, there appears to be a clear cluster in the north (in red) and three clusters in central Mexico City (teal, blue, and gold). In December, we see three clusters with strong spatial patters: red is concentrated in the east, purple in the southeast, and gold in central Mexico City. Within each month, these maps seem to indicate that the clustered regression coefficients capture spatial similarity that our Gaussian process model does not.

Again, we use $3$ GP factors, with decay parameters fixed, as above. To examine the relative importance of these factors, we compute and plot the relative Euclidean norms of the columns of \mbox{\boldmath\Lambda}^{(k)} in Figure 7. Interestingly, the long and short range factors account for 93.4% and 92.2% of the factor model for ozone and $\text{PM}_{10}$ , respectively. Therefore, three factors appear to be sufficient in this application.

In this analysis, the mean number of clusters varies somewhat month-to-month. To illustrate this, we plot the number of clusters as a function of month in Figure 8. Averaging over all months and iterations of our sampler, the mean number of clusters is 5.75. We note that clustering between sites is driven by ozone in this example. In fact, when we use independent DPs for ozone and $\text{PM}_{10}$ , there is no clustering for $\text{PM}_{10}$ . We also found that the number of clusters was relatively robust to the prior distribution on error terms $\tau^{2}_{k}$. However, we did see that the number of clusters was somewhat sensitive to the selection of the concentration parameter $\alpha$ and the base measure for regression coefficients. We discuss this in more detail in Appendix C.

Finally, we examine the regression coefficients for ozone and $\text{PM}_{10}$ . These should be interpreted in the context of our multivariate GP model with time-varying intercepts and seasonal terms. In Figure 9, we plot the cluster-specific posterior means for the regression coefficients for the effect of temperature and relative humidity on temperature. These show significant variability over the course of the year and within the month (i.e., between clusters). For January to May, temperature has a positive relationship with ozone for most clusters, but there is significant variability during the rest of the year. Humidity appears to have a negative relationship for most clusters from March to October. Temperature and $\text{PM}_{10}$ generally have a positive relationship from March to November. Relative humidity has a positive relationship with $\text{PM}_{10}$ during most of the year.

5 Conclusions and Future Work

In this paper, we have introduced a modeling framework for multivariate functional data over time and across sites with the goal of time-varying clustering on regression coefficients. We presented a continuous time scalable multivariate GP model that uses factor analysis to capture dependence between sites and coregionalization to capture dependence within sites. leverages a factor model to capture inter-site, co-regionalization to model covariance between variables. We used the exponential correlation function for temporal autocorrelation, resulting in a continuous time first-order Markov process. We used an exchangeable time-varying Dirichlet process prior distributions for regression coefficients to cluster sites that have a similar joint response to exogenous variables (covariates). The resulting hierarchical model was fitted within a Bayesian framework whence we presented appropriate prior distribution and associated model fitting details. Lastly, we analyzed our motivating dataset: hourly ozone and $\text{PM}_{10}$ data from Mexico City in 2017.

Turning to future work, we could allow component specific partitions. In fact, within the Dirichlet process mixing framework, we could attempt atom, i.e., label selection in continuous time following ideas in Petrone et al., (2009). In a second direction, rather than assuming that DPs were exchangeable month-to-month, we could incorporate dependence between base measures through a dynamic linear model (West and Harrison,, 1997). As a third possibility, though spatial autocorrelation was not a significant component in our motivating example, in some datasets sites may exhibit spatial autocorrelation after accounting for covariates and temporal autocorrelation. In these cases, we could incorporate spatial autocorrelation between factor loadings (see Christensen and Amemiya,, 2002; Hogan and Tchernis,, 2004, for early discussion on this).

Appendix A Distributional Properties

First, we derive the result in (9). Following the notation in Section 2.1, we seek $\text{cov}(\eta_{i}^{k}(t),\eta_{i^{\prime}}^{k^{\prime}}(t^{\prime})$ for say $t^{\prime}\geq t$. We have

[TABLE]

where $\phi_{ij}^{(k)}(t)=\sum_{l=1}^{r}\Delta_{il}^{(k)}\nu_{l}^{(j)}(t)$ and $\phi_{i^{\prime}j^{\prime}}^{(k^{\prime})}(t^{\prime})=\sum_{l=1}^{r}\Delta_{i^{\prime}l}^{(k^{\prime})}\nu_{l}^{(j^{\prime})}(t^{\prime})$. Since the $\nu_{l}^{(j)}(t)$ are independent across $l$ and $j$, $\text{cov}(\phi_{ij}^{(k)}(t),\phi_{i^{\prime}j^{\prime}}^{(k^{\prime})}(t^{\prime}))=0$ if $l\neq l^{\prime}$ or $j\neq j^{\prime}$. However, with $j=j^{\prime}$, we have $\text{cov}(\phi_{ij}^{(k)}(t),\phi_{i^{\prime}j}^{(k^{\prime})}(t^{\prime})=\sum_{l=1}^{r}\Delta_{il}^{(k)}\Delta_{i^{\prime}l}^{(k^{\prime})}e^{-\phi_{(j)}(t^{\prime}-t)}$.

So, $\sum_{j=1}^{K}\sum_{j^{\prime}=1}^{K}A_{kj}A_{k^{\prime}j^{\prime}}\text{cov}(\phi_{ij}^{(k)}(t),\phi_{i^{\prime}j^{\prime}}^{(k^{\prime})}(t^{\prime})=\sum_{j=1}^{K}A_{kj}A_{k^{\prime}j}\sum_{l=1}^{r}\Delta_{il}^{(k)}\Delta_{i^{\prime}l}^{(k^{\prime})}e^{-\phi_{(j)}(t^{\prime}-t)}$. After a little algebra this becomes (\mbox{\boldmath\Lambda}_{i}^{(k)})^{T}\mbox{\boldmath\Lambda}_{i^{\prime}}^{(k^{\prime})}\sum_{j=1}^{K}A_{kj}A_{k^{\prime}j}e^{-\phi_{(j)}(t^{\prime}-t)}.

Special cases of (9) used below are the following:

(i) when $t^{\prime}=t$, \text{cov}(\eta_{i}^{k}(t),\eta_{i^{\prime}}^{k^{\prime}}(t^{\prime})=(\mbox{\boldmath\Lambda}_{i}^{(k)})^{T}\mbox{\boldmath\Lambda}_{i^{\prime}}^{(k^{\prime})}\sum_{j=1}^{K}A_{kj}A_{k^{\prime}j},

(ii) when $i^{\prime}=i$, \text{cov}(\eta_{i}^{k}(t),\eta_{i^{\prime}}^{k^{\prime}}(t^{\prime})=(\mbox{\boldmath\Lambda}_{i}^{(k)})^{T}\mbox{\boldmath\Lambda}_{i}^{(k^{\prime})}\sum_{j=1}^{K}A_{kj}^{2}e^{-\phi_{(j)}(t^{\prime}-t)},

(iii) when $k^{\prime}=k$, \text{cov}(\eta_{i}^{k}(t),\eta_{i^{\prime}}^{k}(t^{\prime})=(\mbox{\boldmath\Lambda}_{i}^{(k)})^{T}\mbox{\boldmath\Lambda}_{i^{\prime}}^{(k)}\sum_{j=1}^{K}A_{kj}A_{k^{\prime}j}e^{-\phi_{(j)}(t^{\prime}-t)}.

So, the marginal distribution of $Y_{i}^{(k)}(t)$ is $Y_{i}^{(k)}(t)\sim N(\mu_{i}^{(k)}(t),(\Delta_{i}^{(k)})^{T}\Delta_{i}^{(k)}\sum_{j=1}^{K}A_{kj}^{2})$. Therefore, with \mu_{i}^{(k)}(t)=\textbf{X}_{i}(t)^{T}\mbox{\boldmath\beta}_{i}^{(k)}, for clustering across $i$, the levels of the $\textbf{X}_{i}(t)$ do not matter. We only ask whether \mbox{\boldmath\beta}_{i}^{(k)}=\mbox{\boldmath\beta}_{i^{\prime}}^{(k)}?

The marginal distribution for $\textbf{Y}_{i}(t)$ is also directly available, i.e., \textbf{Y}_{i}(t)\sim MVN(\mbox{\boldmath\mu}_{i}(t),\mbox{\boldmath\Sigma}_{ii}+\textbf{D}) where D is a $K\times K$ diagonal matrix with $(\textbf{D})_{kk}=\tau^{2(k)}$ and \mbox{\boldmath\Sigma}_{ii} also $K\times K$ with entries (\mbox{\boldmath\Sigma}_{ii})_{kk^{\prime}} equal to (ii) above with $t^{\prime}=t$.

The marginal distribution for the entire $nK\times 1$ vector \textbf{Y}(t)=\left(\begin{array}[]{c}\textbf{Y}_{1}(t)\\ \textbf{Y}_{2}(t)\\ .\\ .\\ .\\ \textbf{Y}_{n}(t)\\ \end{array}\right) is

[TABLE]

where I is the $n\times n$ identity matrix, D is as above, and $\Sigma$ is $nK\times nK$ with $K\times K$ blocks such that block \mbox{\boldmath\Sigma}_{ii^{\prime}} has entries, (\mbox{\boldmath\Sigma}_{ii^{\prime}})_{kk^{\prime}} equal to (i) above.

Finally, the conditional distribution, $[\textbf{Y}(t^{\prime})|\textbf{Y}(t)]$ can be obtained from the joint distribution of $(\textbf{Y}(t),\textbf{Y}(t^{\prime}))$. That is, \left(\begin{array}[]{c}\textbf{Y}(t)\\ \textbf{Y}(t^{\prime})\\ \end{array}\right)\sim MVN (\left(\begin{array}[]{c}\mbox{\boldmath\mu}(t)\\ \textbf{m}(t^{\prime})\\ \end{array}\right),\left(\begin{array}[]{cc}\mbox{\boldmath\Sigma}+\textbf{I}\otimes\textbf{D}&\mbox{\boldmath\Sigma}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})}\\ \mbox{\boldmath\Sigma}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})}&\mbox{\boldmath\Sigma}+\textbf{I}\otimes\textbf{D}\\ \end{array}\right)) where \mbox{\boldmath\Sigma}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})} is $nK\times nK$ symmetric with $K\times K$ blocks and block (\mbox{\boldmath\Sigma}_{\textbf{Y}(t),\textbf{Y}(t^{\prime})})_{ii^{\prime}} has entry $k,k^{\prime}$ as (9). Therefore,

[TABLE]

[TABLE]

Appendix B Gibbs Sampling

Here, we provide the Gibbs sampling for a slightly more general model than that presented between (4)-(8). In addition, we use the prior parameters denoted in Section 3.1. As we do with the Mexico City data, we allow limiting cluster based on a subset of the explanatory variables

[TABLE]

where the $k$th row of $\textbf{G}^{m}_{i}$ is \mbox{\boldmath\gamma}_{im}^{(k)}. Let $\theta|\cdots$ denote the posterior conditional distribution for $\theta$. With the exception of decay parameters of the exponential covariance functions, $\phi_{(j)}$, the posterior conditional distributions of all model parameters are conveniently available and presented below. Let $c$ indicate cluster and $n_{cm}$ be the number of clusters in time segment $m$. In addition, we constrain the coregionalization matrix A to be lower diagonal and index matrix elements $A_{kl}$ for $k<l$.

[TABLE]

where

[TABLE]

To update the Dirichlet process parameters, we use a collapsed sampler (Escobar and West,, 1995). To do this, we find the marginal distribution of the data for the at location $i$ and time block $m$, integrating over \mbox{\boldmath\beta}^{k}_{im} with respect to its prior distribution. This allows to compute the probability of creating a new group or cluster. Here, we let \mu^{(k)}_{i,-\beta}(t)=\textbf{z}_{i}(t)^{T}\mbox{\boldmath\gamma}_{im}^{(k)}+\eta_{i}^{(k)}(t) and $n_{c,-i}$ be the number stations in the $c$th if we exclude the $i$th location. Then, the relative probabilities of each cluster is

[TABLE]

where

[TABLE]

Using the Woodbury matrix formula, we can efficiently computate the inverse of the covariance matrix and the determinant.

For the GP parameters, we define a few terms that allow us to make an $r\times 1$ update. Each functional factor is defined using $K$ independent GPs. First, we define two diagonal matrices $D_{v}^{(k)}=\text{diag}\left(\frac{1}{1-e^{-2\phi_{1}^{(k)}}},...,\frac{1}{1-e^{-2\phi_{r}^{(k)}}}\right)$ and $D_{\phi}^{(k)}=\text{diag}\left(e^{-\phi_{1}^{(k)}},...,e^{-\phi_{r}^{(k)}}\right)$. In our example, $D_{v}^{(k)}=D_{v}^{(k^{\prime})}$ for all $k$ and $k^{\prime}$. Then, for $1<t<T$,

[TABLE]

When, $t=1$, $D_{v}^{(k)}+D_{\phi}^{(k)}D_{v}^{(k)}D_{\phi}^{(k)}$ becomes $\textbf{I}+D_{v}^{(k)}$ and there are no past terms in the mean. For $t=T$, $\Sigma_{\nu^{(k)}(t)}$ excludes $D_{\phi}^{(k)}D_{v}^{(k)}D_{\phi}^{(k)}$, as well as future terms in the mean.

If desired, coefficients, \mbox{\boldmath\gamma}_{im}^{(k)} and \mbox{\boldmath\beta}_{cm}^{(k)}, can be specified hierarchically or dynamically. In this case, the Gibbs sampler would include prior distributions and updates for $m_{\gamma^{(k)}_{m}}$, $m_{\beta^{(k)}_{m}}$,$V_{\gamma^{(k)}_{m}}$, and $V_{\beta^{(k)}_{m}}$.

Appendix C Sensitivity Analysis

Here, we take up the sensitivity of our model to the concentration parameter $\alpha$ and the prior distributions on the noise parameters $\tau^{2}_{k}$ and \mbox{\boldmath\beta}^{(k)}_{im}. First, we present how sensitive clustering is to prior specifications of noise parameters $\tau^{2}$ and concentration parameters. To do this, we supply the mean and median number of clusters as a function of month for many combinations of prior distributions for the error variance and the fixed concentration parameter $\alpha$ (See Figure 10). We provide the average number of clusters averaged over month and over iterations of the Gibbs sampler in Table 1. In general, with larger concentration parameters, the number of clusters increases. With more diffuse prior distributions on $\tau^{2}$ parameters, we generally had weaker clustering (more clusters); however, this depended somewhat upon the concentration parameter.

When creating a new cluster, the marginal likelihood, integrated over the base measure, must be high (See Appendix B). With a diffuse base measure, the marginal likelihood is lower and fewer clusters are expected. Because the variance of the base measure affects clustering, we explore how different base measures change the observed number of clusters. We plot the average number of clusters over month of the year in Figure 11. For the highest variance base measure, $\mathcal{N}\left(\textbf{0},10^{3}\textbf{I}\right)$, we observe 4.75 clusters, averaged over month and iteration. For base measures of $\mathcal{N}\left(\textbf{0},\textbf{I}\right)$ and $\mathcal{N}\left(\textbf{0},10^{-3}\textbf{I}\right)$, we observe the average number of clusters to be 5.29 and 6.90, respectively. As expected, when the variance of the base measure is small, we have more clusters, on average. On the other hand, a very diffuse base measure leads to fewer clusters. This results jibes with our intuition. Given that we have centered and scaled the outcomes and covariates, a very diffuse base measure for the regression coefficients seems unreasonable.

We also briefly discuss our decision to jointly cluster regression coefficients for ozone and $\text{PM}_{10}$ . In preliminary modeling, we noticed that regression coefficients for $\text{PM}_{10}$ would not cluster with reasonable values of $\alpha$. Ozone, on the other hand, would cluster. When the concentration parameter for both ozone and $\text{PM}_{10}$ was 1, the mean number of clusters averaging over month and iterations of the Gibbs sampler were 3.98 and 24 for ozone and $\text{PM}_{10}$ , respectively. Therefore, we argue that the joint clustering we have discussed in Section 4.3 is driven by clustering in ozone coefficients. In Figure 12, we plot the mean and median number of ozone clusters as a function of month.