Nonparametric collective spectral density estimation with an application to clustering the brain signals

TL;DR

This paper introduces a nonparametric method for jointly estimating spectral density functions of multiple time series, enabling improved estimation, visualization, and clustering of brain signals based on their spectral features.

Contribution

It proposes a collective spectral density estimation technique using basis functions and a low-dimensional manifold, enhancing efficiency and enabling clustering of brain signals.

Findings

Improved spectral density estimation accuracy.

Effective clustering of EEG signals based on spectral features.

Visualization tools for spectral density analysis.

Abstract

In this paper, we develop a method for the simultaneous estimation of spectral density functions (SDFs) for a collection of stationary time series that share some common features. Due to the similarities among the SDFs, the log-SDF can be represented using a common set of basis functions. The basis shared by the collection of the log-SDFs is estimated as a low-dimensional manifold of a large space spanned by a pre-specified rich basis. A collective estimation approach pools information and borrows strength across the SDFs to achieve better estimation efficiency. Also, each estimated spectral density has a concise representation using the coefficients of the basis expansion, and these coefficients can be used for visualization, clustering, and classification purposes. The Whittle pseudo-maximum likelihood approach is used to fit the model and an alternating blockwise Newton-type algorithm is developed for the computation. A web-based shiny App found at "https://ncsde.shinyapps.io/NCSDE" is developed for visualization, training and learning the SDFs collectively using the proposed technique. Finally, we apply our method to cluster similar brain signals recorded by the electroencephalogram for identifying synchronized brain regions according to their spectral densities.

Figures

Tables

Table 1. Table 1 : Comparison based on 100 100 100 simulation runs. The mean and standard errors (in parentheses) of the angle are reported.

ts-$n$-$m$	Ps	tSVD.Ps	NSDE	tSVD.NSDE	NCSDE
ts-100-06	84.97(0.38)	78.47(0.77)	88.20(0.17)	68.51(1.74)	37.50(1.98)
ts-100-15	89.13(0.08)	80.50(0.71)	90.00(0.00)	79.10(0.80)	33.84(2.64)
ts-100-30	88.93(0.09)	80.66(0.73)	90.00(0.00)	79.10(0.70)	24.12(1.31)
\hdashlinets-200-06	84.93(0.43)	75.40(0.85)	89.65(0.04)	59.34(2.14)	20.36(0.50)
ts-200-15	88.29(0.14)	78.40(0.95)	90.00(0.00)	51.74(2.28)	15.50(0.89)
ts-200-30	89.38(0.05)	80.59(0.69)	90.00(0.00)	66.85(1.53)	11.80(0.32)
\hdashlinets-400-06	85.51(0.36)	66.16(1.30)	89.90(0.01)	43.12(1.72)	13.91(0.55)
ts-400-15	87.66(0.18)	71.88(1.16)	90.00(0.00)	38.41(1.53)	12.39(0.34)
ts-400-30	88.92(0.08)	73.07(1.18)	90.00(0.00)	34.02(1.45)	7.76(0.17)

Table 2. Table 2 : Comparison based on 100 100 100 simulation runs. The mean and standard errors (in parentheses) of the adjusted Rand similarity coefficient ( 𝖠𝖱𝖨 𝖠𝖱𝖨 \mathsf{ARI} ) are reported.

ts-$n$-$m$	Ps	tSVD.Ps	NSDE	tSVD.NSDE	NCSDE
ts-100-06	0.651(0.012)	0.631(0.010)	0.686(0.014)	0.691(0.014)	0.914(0.012)
ts-100-15	0.542(0.007)	0.559(0.005)	0.483(0.007)	0.483(0.007)	0.952(0.008)
ts-100-30	0.526(0.006)	0.550(0.003)	0.464(0.005)	0.464(0.005)	0.959(0.006)
\hdashlinets-200-06	0.639(0.010)	0.640(0.011)	0.742(0.016)	0.731(0.016)	0.996(0.003)
ts-200-15	0.559(0.008)	0.570(0.006)	0.784(0.020)	0.681(0.022)	0.994(0.003)
ts-200-30	0.526(0.005)	0.547(0.003)	0.636(0.018)	0.567(0.015)	0.998(0.001)
\hdashlinets-400-06	0.637(0.011)	0.689(0.012)	0.773(0.016)	0.773(0.016)	0.998(0.002)
ts-400-15	0.545(0.008)	0.598(0.006)	0.811(0.015)	0.822(0.014)	1.000(0.000)
ts-400-30	0.534(0.006)	0.581(0.004)	0.872(0.014)	0.884(0.012)	1.000(0.000)

Equations

$$\gamma (h)={\int }_{-1/2}^{1/2}f(\omega )e^{2\pi i\omega h}d\omega ,$$

$$f(\omega )=\sum _{h=-\infty }^{\infty }\gamma (h)e^{-2\pi i\omega h},-1/2\le \omega \le 1/2.$$

$$I_n({\omega }_j)=\sum _{h=-(n-1)}^{n+1}\hat{\gamma }(h)e^{-2\pi i{\omega }_{j}h},j=0,1,...,n-1,$$

$$L(\mathbf{\theta })=\frac{1}{(2\pi {)}^{n/2}|{\Gamma }_{n,\mathbf{\theta }}{|}^{1/2}}\exp \left(-\frac{1}{2}{\mathbf{X}}^{\top }{\Gamma }_{n,\mathbf{\theta }}^{-1}\mathbf{X}\right),$$

$$\ell (\mathbf{\theta })\propto log(|{\Gamma }_{n,\mathbf{\theta }}|)+\mathrm{tr}\left(\mathbf{X}{\mathbf{X}}^{\top }{\Gamma }_{n,\mathbf{\theta }}^{-1}\right).$$

$${\ell }_W(\mathbf{\theta })=n{\int }_{-1/2}^{1/2}\left\{log(f_{\mathbf{\theta }}(\omega ))+I_n(\omega ){f_{\mathbf{\theta }}(\omega )}^{-1}\right\}d\omega ,$$

$$\ell_{W}({\boldsymbol{\theta}})\stackrel{{\scriptstyle.}}{{=}}n\sum_{-1/2<\omega_{j}<1/2}\Bigl{\{}\hbox{log}({f_{\boldsymbol{\theta}}(\omega_{j})})+{I_{n}(\omega_{j})}{f_{\boldsymbol{\theta}}(\omega_{j})}^{-1}\Bigr{\}}.$$

$$u_i(\omega )=\sum _{k=1}^K{\varphi }_k(\omega ){\alpha }_{ik},i=1,\dots ,m.$$

$$f_{i}(\omega)=\exp{u_{i}(\omega)}=\exp\Biggl{\{}\sum_{k=1}^{K}\phi_{k}(\omega)\alpha_{ik}\Biggr{\}},\quad i=1,\dots,m.$$

$${\varphi }_k(\omega )=\sum _{\ell =1}^L{b}_{\ell }(\omega ){\theta }_{\ell k},k=1\dots ,K.$$

$$u_i(\omega )=\mathbf{\varphi }(\omega {)}^{\top }{\mathbf{\alpha }}_i=\mathbf{b}(\omega {)}^{\top }\mathbf{\Theta }{\mathbf{\alpha }}_i,i=1,\dots ,m.$$

$$u_i(\omega )=\mathbf{b}(\omega {)}^{\top }{\mathbf{\psi }}_i,i=1,\dots ,m.$$

$$I_{i,n}(\omega)={\Bigl{|}\displaystyle{{1}\over{n}}\sum_{t=1}^{n}X_{it}e^{-2\pi it\omega}\Bigr{|}^{2}}.$$

$$\ell_{W}({\boldsymbol{\Theta}},{\bf A})=\sum_{i=1}^{m}\sum_{{j}\in\bm{J}}\Biggl{\{}u_{i}(\omega_{j})+I_{i,n}(\omega_{j})\exp\bigl{[}-u_{i}(\omega_{j})\bigr{]}\Biggr{\}},$$

$$-2{\ell }_W(\mathbf{\Theta },\mathbf{A})+\lambda \sum _{k=1}^K\mathsf{PEN}({\varphi }_k),$$

$$\sum _{k=1}^K\mathsf{PEN}({\varphi }_k)=\sum _{i=k}^K{\mathbf{\theta }}_k^{\top }\mathbf{R}{\mathbf{\theta }}_k=\mathrm{tr}\{{\mathbf{\Theta }}^{\top }\mathbf{R}\mathbf{\Theta }\}.$$

$${\mathbf{\alpha }}_i^{new}$$

$$\displaystyle{\boldsymbol{\theta}}_{k}^{new}\ ={\boldsymbol{\theta}}_{k}^{old}-\tau\,\biggl{[}\frac{\partial^{2}}{\partial{\boldsymbol{\theta}}_{k}\partial{\boldsymbol{\theta}}_{k}^{\top}}\{\ell_{W}({\boldsymbol{\Theta}},{\bf A})\}-\lambda{\bf R}\biggr{]}^{-1}\times\qquad\qquad\qquad\qquad\qquad\qquad$$

$$\displaystyle\qquad\qquad\qquad\biggl{[}\frac{\partial}{\partial{\boldsymbol{\theta}}_{k}}\{\ell_{W}({\boldsymbol{\Theta}},{\bf A})\}-\lambda{\bf R}{\boldsymbol{\theta}}_{k}\biggr{]}\bigg{|}_{{\boldsymbol{\Theta}}={\boldsymbol{\Theta}}^{old},{\bf A}={\bf A}^{old}},$$

$$\displaystyle={\boldsymbol{\theta}}_{k}^{old}-\tau\biggl{[}\sum_{i=1}^{m}\alpha_{ik}^{2}\sum_{j}\Bigl{\{}{\bf b}(\omega_{j})I_{i,n}(\omega_{j})\exp\bigl{[}-u_{i}(\omega_{j})\bigr{]}{\bf b}(\omega_{j})^{\top}\Bigr{\}}-\lambda{\bf R}\biggr{]}^{-1}\times\qquad$$

$$\displaystyle\qquad\biggl{[}\sum_{i=1}^{m}\alpha_{ik}\sum_{j}\Bigl{\{}{\bf b}(\omega_{j})-{\bf b}(\omega_{j})I_{i,n}(\omega_{j})\exp\bigl{[}-u_{i}(\omega_{j})\bigr{]}\Bigr{\}}-\lambda{\bf R}{\boldsymbol{\theta}}_{k}\biggr{]}\bigg{|}_{{\boldsymbol{\Theta}}={\boldsymbol{\Theta}}^{old},{\bf A}={\bf A}^{old}}$$

$$\mathrm{AIC}(\lambda )=-2{\ell }_W(\hat{\mathbf{\Theta }},\hat{\mathbf{A}})+2\mathsf{df}(\lambda ),$$

$$\displaystyle{\sf df}(\lambda)=\sum_{k=1}^{K}{\rm trace}\biggl{\{}\biggl{[}\sum_{i=1}^{m}\alpha_{ik}^{2}\sum_{j}\Bigl{\{}{\bf b}(\omega_{j})I_{i,n}(\omega_{j})\exp\bigl{[}-u_{i}(\omega_{j})\bigr{]}{\bf b}(\omega_{j})^{\top}\Bigr{\}}-\lambda{\bf R}\biggr{]}^{-1}\times\biggr{.}$$

$$\displaystyle\biggl{.}\biggl{[}\sum_{i=1}^{m}\alpha_{ik}^{2}\sum_{j}\Bigl{\{}{\bf b}(\omega_{j})I_{i,n}(\omega_{j})\exp\bigl{[}-u_{i}(\omega_{j})\bigr{]}{\bf b}(\omega_{j})^{\top}\Bigr{\}}\biggr{]}\biggr{\}}.\qquad\quad$$

$${\hat{\lambda }}_{i+1}^{-1}=\frac{trace({\hat{\mathbf{\Theta }}}_i^{\top }\mathbf{R}{\hat{\mathbf{\Theta }}}_i)}{\mathsf{df}({\hat{\lambda }}_i)-(a-1)},$$

$${\Delta }_2{\mathbf{\theta }}_k:={\Delta }_1{\Delta }_1{\mathbf{\theta }}_k={\mathbf{\theta }}_k-2{\mathbf{\theta }}_{k-1}+{\mathbf{\theta }}_{k-2}.$$

$$\mbox{\boldmath$L$}_{1}={\left[{\begin{array}[]{*{20}{c}}1&{-1}&0&\ldots&0\\ 0&1&{-1}&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&0\\ 0&\cdots&0&1&{-1}\end{array}}\right]_{(M-1\times M)}}$$

$$\{\begin{array}{ccc}\mathrm{Model}\mathrm{I}& AR(3):\mathbf{\phi }=(0.1,0.5,0.1),& \mathrm{with}prob.p_1,\\ \mathrm{Model}\mathrm{II}& AR(3):\mathbf{\phi }=(0.1,0.1,0.5),& \mathrm{with}prob.p_2,\\ \mathrm{Model}\mathrm{III}& AR(3):\mathbf{\phi }=(0.5,0.1,0.1),& \mathrm{with}prob.p_3.\end{array}$$

$$f(\omega )={\sigma }_{\epsilon }^2{\left|1-\sum _{k=1}^3{\phi }_k{e}^{-2\pi ik\omega }\right|}^{-2},$$

$$I_{i,n}({\boldsymbol{\omega}})={\Bigl{|}\displaystyle{{1}\over{n}}\sum_{t=1}^{n}X_{it}e^{-2\pi it{\boldsymbol{\omega}}}\Bigr{|}^{2}},$$

$$\mathsf{ARI}=\frac{\sum_{i=0}^{1}\sum_{j=0}^{1}\binom{n_{ij}}{2}-\Bigl{[}\sum_{i}\binom{n_{i\cdot}}{2}+\sum_{j}\binom{n_{\cdot j}}{2}\Bigr{]}/\binom{m}{2}}{\frac{1}{2}\Bigl{[}\sum_{i}\binom{n_{i\cdot}}{2}+\sum_{j}\binom{n_{\cdot j}}{2}\Bigr{]}-\Bigl{[}\sum_{i}\binom{n_{i\cdot}}{2}+\sum_{j}\binom{n_{\cdot j}}{2}\Bigr{]}/\binom{m}{2}}.$$

Full text

Nonparametric collective spectral density estimation with an application to clustering the brain signals

Mehdi Maadooliat111Department of Mathematics, Statistics and Computer Science, Marquette University, Milwaukee, WI 53201-1881 USA. E-mail: mehdi@mscs.mu.edu, Ying Sun222CEMSE Division, King Abdullah University of Science and Technology, Thuwal 23955-6900 Saudi Arabia. E-mail: ying.sun@kaust.edu.sa and Tianbo Chen2

Abstract

In this paper, we develop a method for the simultaneous estimation of spectral density functions (SDFs) for a collection of stationary time series that share some common features. Due to the similarities among the SDFs, the log-SDF can be represented using a common set of basis functions. The basis shared by the collection of the log-SDFs is estimated as a low-dimensional manifold of a large space spanned by a pre-specified rich basis. A collective estimation approach pools information and borrows strength across the SDFs to achieve better estimation efficiency. Also, each estimated spectral density has a concise representation using the coefficients of the basis expansion, and these coefficients can be used for visualization, clustering, and classification purposes. The Whittle pseudo-maximum likelihood approach is used to fit the model and an alternating blockwise Newton-type algorithm is developed for the computation. A web-based shiny App found at “https://ncsde.shinyapps.io/NCSDE” is developed for visualization, training and learning the SDFs collectively using the proposed technique. Finally, we apply our method to cluster similar brain signals recorded by the electroencephalogram for identifying synchronized brain regions according to their spectral densities.

MSC 2010 subject classifications: Primary 62G05, 62M15; secondary 62P10.

Keywords and phrases: Collective estimation; nonparametric estimation; roughness penalty; time series clustering; spectral density function; Whittle likelihood

Short title: Nonparametric collective spectral density estimation

1 Introduction

The spectral density plays an important role in time series analysis in frequency domain. We suppose that $X_{t}$ represents a zero-mean, weakly stationary time series with the autocovariance function $\gamma(h)$ defined as $\gamma(h)=E(X_{t}X_{t+h}),$ $h=0,\pm 1,\pm 2,\ldots.$ If the autocovariance function is absolutely summable, i.e., $\sum_{h=-\infty}^{\infty}|\gamma(h)|<\infty,$ then the autocovariance sequence $\gamma(h)$ has the spectral representation

[TABLE]

where $f(\omega)$ is the spectral density of $X_{t}$, which has the inverse Fourier representation

[TABLE]

Given the time series $\{x_{t}$, $t=1,\ldots,n\}$, when ${\gamma}(h)$ is replaced by the sample covariance $\hat{\gamma}(h)$, the periodogram is defined as

[TABLE]

which can be calculated as $I_{n}(\omega_{j})=|d(\omega_{j})|^{2}$, where $\displaystyle d(\omega_{j})=n^{-1/2}\sum_{t=1}^{n}x_{t}e^{-2\pi i\omega_{j}t}$ is the discrete Fourier transform of $x_{t}$ at the fundamental frequencies $\omega_{j}=j/n$.

However, the raw periodogram is not a consistent estimator for the spectral density of a stationary random process. One classical method to obtain a consistent estimator is to smooth the periodogram across frequencies. Yuen (1979) analyzed the performance of three methods of periodogram smoothing for spectrum estimation. Wahba (1980) developed an objective optimum smoothing procedure for estimating the log-spectral density using the spline to smooth the log-periodogram. A discrete spectral average estimator and lag window estimators were introduced in Brockwell and Davis (2013). Both of the two methods are consistent. Brillinger (2001) introduced periodogram kernel smoothing. One critical issue in periodogram smoothing is span selection. Lee (1997) proposed a span selector based on unbiased risk estimation. Ombao et al. (2001) proposed using the Gamma-deviance generalized cross-validation (Gamma GCV), and Lee (2001) proposed a bandwidth selection based on coupling of the so-called plug-in and the unbiased risk estimation ideas.

Another popular method for spectral density estimation is based on likelihoods. For example, Capon (1983) used the maximum likelihood method to estimate the spectral density of signals with noise. In Güler et al. (2001), the spectral density of brain signal data was analyzed using the maximum likelihood. The well-known Whittle likelihood was developed for time series analysis in Whittle (1953), Whittle (1954b), and Whittle (1962). Pawitan and O’Sullivan (1994) proposed using a penalized Whittle likelihood to estimate spectral density functions. Chow and Grenander (1985) developed a penalized likelihood-type method for the nonparametric estimation of the spectral density of Gaussian processes. Fan and Kreutzberger (1998) used local polynomial techniques to fit the Whittle likelihood for spectral density estimation. We suppose there is a time series ${\bf X}$ of length $n$ from a mean-zero stationary Gaussian process with a parametric covariance function $\gamma(h;{\boldsymbol{\theta}})$, where ${\boldsymbol{\theta}}$ is the unknown parameter. If we let $\Gamma_{n,{\boldsymbol{\theta}}}$ be the covariance matrix of the random vector ${\bf X}$, then the likelihood function is

[TABLE]

and $\ell({\boldsymbol{\theta}})=-2\times\hbox{log}L({\boldsymbol{\theta}})$ is the log-likelihood,

[TABLE]

Whittle (1954a) proposed an approximation of the above log-likelihood function, known as the Whittle likelihood approximation:

[TABLE]

where ${f_{\boldsymbol{\theta}}(\omega)}$ is the spectral density function. For a discrete frequency range, the Whittle approximation (1) can be written as

[TABLE]

Using the fast Fourier transform (FFT), $\ell_{W}({\boldsymbol{\theta}})$ can be obtained efficiently with only $O\left(n\hbox{log}_{2}n\right)$ operations.

We apply our method to the neuroscience study of functional connectivity in the brain with electroencephalgram (EEG) data. The EEG is a method of monitoring spontaneous electrical activity in the brain over a period of time. EEGs are typically recorded from multiple electrodes, referred as EEG channels, placed on the scalp. The clustered EEG channels are useful for understanding how different brain regions communicate with each other. To identify synchronized EEG channels that show similar oscillations or waveforms, many time series clustering algorithms have been developed to cluster similar smoothed periodograms estimated from the EEG channels separately (see Hennig et al., 2015, for some examples).

In this paper, we propose collectively estimating multiple spectral densities that share common features. The collective estimation was first proposed by Maadooliat et al. (2016) for the estimation of probability density functions in protein structure analysis. For multiple stationary time series that may share similar spectral density functions, we propose a nonparametric spectral density estimation approach that minimizes a penalized Whittle likelihood. The dimension is reduced by representing the raw periodograms on a set of common basis functions while allowing each one to differ by introducing random effects into the model. The collectively estimated spectral densities can then be used for clustering purposes. Collective estimation is a novel approach for estimating multiple spectral densities that share common features in a statistically more efficient way.

The rest of the paper is organized as follows. Section 2 presents the core of the proposed method. Subsection 2.1 provides the blockwise Newton-Raphson algorithm, and subsections 2.2-2.5 provide the implementation details. Section 3 reports the simulation results to illustrate the proposed collective estimation approach and to compare it with competitive non-collective estimation approaches. The application to EEG data is given in Section 4. Section 5 concludes the paper.

2 Methodology

We consider a collection of stationary time series, ${\bf X}_{i}$s, where $i=1,\cdots,m$ and  ${\bf X}_{i}^{\top}=(X_{i1},X_{i2},\cdots,X_{in})$ with the associated spectral density $f_{i}$. We want to estimate the spectral density functions together. The rationale of this collective spectral density estimation approach is to improve the efficiency of the estimates by using a shared basis to obtain the SDFs.

We assume that each log-spectral density function can be represented by a linear combination of a common set of basis functions $\{\phi_{k}(\omega),k=1,\cdots,K\}$, and that each has its own set of coefficients. More specifically, we assume that $\hbox{log}\{f_{i}(\omega)\}=u_{i}(\omega)$, with

[TABLE]

Equivalently, the spectral density functions can be written as

[TABLE]

For identifiability, we require that $1,\phi_{k},\ k=1,\dots,K$, to be linearly independent. We would like $K$ to be a small number so that the number of parameters to be estimated remains manageable even when we estimate a large number of spectral densities (i.e., when $m$ is large).

In our setting, the basis functions are not prespecified and need to be determined from the data. To this end, we suppose that these basis functions fall in a low-dimensional subspace of a function space spanned by a rich family of fixed basis functions, $\{b_{\ell}(\omega),\ell=1,\cdots,L\}$ ($L\gg K$), such that

[TABLE]

For identifiability, we require that $1,\ b_{\ell},\ \ell=1,\dots,L$, to be linearly independent. A large enough $L$ ensures the necessary flexibility to represent the unknown spectral densities. For univariate cases, the fixed basis can be the monomials, the B-splines, or the Fourier basis. Bivariate splines can be used as the fixed basis functions for bivariate spectral densities.

To simplify the presentation, we now introduce some vectors and matrices to denote the quantities of interest: ${\boldsymbol{\phi}}(\omega)=(\phi_{1}(\omega),\cdots,\phi_{K}(\omega))^{\top}$, ${\boldsymbol{\alpha}}_{i}=(\alpha_{i1},\cdots,\alpha_{iK})^{\top}$, ${\bf b}(\omega)=(b_{1}(\omega),\cdots,b_{L}(\omega))^{\top}$, ${\boldsymbol{\theta}}_{k}=(\theta_{1k},\cdots,\theta_{Lk})^{\top}$, and ${\boldsymbol{\Theta}}=({\boldsymbol{\theta}}_{1},\cdots,{\boldsymbol{\theta}}_{K})$. Then, from (2) and (4), we can rewrite $u_{i}(\omega)$ in the vector-matrix form as

[TABLE]

One may combine the equations above for $i=1,\dots,m$ into a matrix form by evaluating the log-SDFs over the discrete frequencies ${\boldsymbol{\omega}}=(\omega_{1},\cdots,\omega_{\tilde{n}})^{\top}$. The $m$ equations given in (5) can be written as ${\bf U}={\bf B}{\boldsymbol{\Theta}}{\bf A}^{\top}$, where ${\bf U}=\left(u_{1}({\boldsymbol{\omega}}),\cdots,u_{m}({\boldsymbol{\omega}})\right)$ is an $\tilde{n}\times m$ matrix that represents the log-SDFs, ${\bf B}=\left({\bf b}(\omega_{1}),\cdots,{\bf b}(\omega_{\tilde{n}})\right)^{\top}$ is an $\tilde{n}\times L$ matrix that represents the rich basis functions at the discrete frequencies ${\boldsymbol{\omega}}$, and ${\bf A}=({\boldsymbol{\alpha}}_{1},\dots,{\boldsymbol{\alpha}}_{m})^{\top}$. The unknown parameters can then be collectively written as the pair $({\boldsymbol{\Theta}},{\bf A})$. There is an identifiability issue caused by the non-uniqueness of the parametrization of $({\boldsymbol{\Theta}},{\bf A})$. This issue can be resolved by introducing some restrictions on the parameterization; see subsection 2.5.

We could have used the fixed basis $\{b_{\ell}(\omega),\ell=1,\cdots,L\}$ in (2) and (3); however, that would be either too restrictive (if $L$ is small) or produce a large number of parameters (if $L$ is large). Alternatively, if we were to model the individual density functions separately using the fixed basis $\{b_{\ell}(\omega),\ell=1,\cdots,L\}$, then we would write

[TABLE]

We let $\mathbf{\Psi}=(\boldsymbol{\psi}_{1},\dots,\boldsymbol{\psi}_{m})$ be the $L\times m$ matrix of coefficients from the basis expansions given in (6). By comparing (5) and (6), we obtain $\mathbf{\Psi}={\boldsymbol{\Theta}}{\bf A}^{\top}$, which is a rank $K$ matrix. Thus, the collective modeling approach introduces a low-rank structure to the coefficient matrix in the basis expansion of the log-spectral densities. This dimensionality reduction allows us to significantly reduce the number of parameters to be estimated and, thus, gain estimation efficiency.

2.1 Estimation using penalized Whittle likelihood

The periodogram $I_{i,n}$ is a rough estimate of the spectral density, $f_{i}$, associated with the time series ${\bf X}_{i}$ observed over $n$ time points:

[TABLE]

We consider the periodogram at the Fourier frequencies $\bm{\omega}=2\pi\bm{J}/n$ for $\bm{J}=\{-\lfloor(n-1)/2\rfloor,\cdots,n-\lfloor n/2\rfloor\}.$

The Whittle likelihood for estimating $m$ spectral densities has the following form:

[TABLE]

where $u_{i}(\omega)$ are defined in (5). It is concave in ${\boldsymbol{\alpha}}_{i}$ when other parameters are fixed and also concave in ${\boldsymbol{\theta}}_{k}$ when other parameters are fixed. Applying the roughness penalty approach of function estimation (Green and Silverman, 1994), we estimate the model parameters by minimizing the penalized likelihood criterion

[TABLE]

where ${\sf PEN}(\phi_{k})$ is a roughness penalty function that regularizes the estimated basis function $\phi_{k}$ to ensure that it is a smooth function, and $\lambda>0$ is a penalty parameter. The penalty function can be written in a quadratic form as

[TABLE]

Two choices for the penalty matrices (${\bf R}_{1}$ and ${\bf R}_{2}$), based on the second derivative of the basis functions and the difference operator, are given in subsection 2.3.

We use an alternating blockwise Newton-Raphson algorithm to minimize the penalized Whittle likelihood approximation. Our algorithm cycles through updating ${\boldsymbol{\alpha}}_{i}$ for $i=1,\dots,m$ and ${\boldsymbol{\theta}}_{k}$ for $k=1,\dots,K$ until convergence. Following the usual step-halving strategy for the Newton-Raphson iteration, the updating formulas are

[TABLE]

and

[TABLE]

where $\tau$ is the first result from the sequence $\{(1/2)^{\delta},\delta=0,1,\dots\}$ such that the objective function in (8) is reduced. The initial values of the Newton-Raphson iteration can be obtained by projecting the raw spectral density estimates (e.g., peridograms) to the model space of (3).

2.2 Selecting the tuning parameter

We may select the penalty parameter by minimizing the AIC (Akaike, 1974),

[TABLE]

where $\ell_{W}({\boldsymbol{\Theta}},{\bf A})$ is the log-likelihood defined in (7), and the degrees of freedom ${\sf df}(\lambda)$ is defined as

[TABLE]

The parameters in these formulas are replaced by their estimated values. The AIC can be derived as an approximation of the leave-one-out cross-validation (O’Sullivan, 1988; Gu, 2002).

Selecting the tuning parameter that minimizes the AIC requires training the model for different values of $\lambda$s and then picking the one that minimizes the criterion function, which can be very expensive in time. Instead, we present an alternative procedure that updates the value of the tuning parameter within the Newton-Raphson iterations. This idea has been used in a generalized mixture model to iteratively update the smoothing parameter (Schall, 1991). Schellhase and Kauermann (2012) and Najibi et al. (2017) extended this approach for density estimation. We borrow their formulation, and use the parameter estimates in the $i^{th}$ step to update the tuning parameter, $\hat{\lambda}_{i+1}$, through

[TABLE]

where $a$ is the order of the differences (derivative) used in the penalty matrix ${\bf R}$ (see Section 2.3). From what we have seen in the implementation of the new procedure, updating the tuning parameter within the Newton-Raphson iterations, on average, does not increase the number of iterations required for convergence. Therefore, the new procedure obtains the final result $p$ times faster than the old procedure, where $p$ is the number of $\lambda$s used in the grid search to minimize the AIC.

2.3 Choices of penalties

2.3.1 Second derivative of the basis

For the univariate spectral density estimation, if we use the usual squared-second-derivative penalty ${\sf PEN}(\phi_{k})=\int\{\phi_{k}^{\prime\prime}(\omega)\}^{2}\,dx$ and $\phi_{k}(\omega)={\bf b}(\omega)^{\top}{\boldsymbol{\theta}}_{k}$, then ${\bf R}_{1}=\displaystyle\int\ddot{{\bf b}}(\omega)\ddot{{\bf b}}(\omega)^{\top}\,d\omega$ with $\ddot{{\bf b}}(\omega)=({b}_{1}^{\prime\prime}(\omega),\dots,b_{L}^{\prime\prime}(\omega))^{\top}$.

2.3.2 Difference operator

We can also control the roughness of the estimated functions by using the difference penalty (Eilers and Marx, 1996) to achieve the appropriate level of smoothness. The variability is controlled through a difference function of order $a$, $\Delta_{a}$, where $\Delta_{1}{\boldsymbol{\theta}}_{k}:={\boldsymbol{\theta}}_{k}-{\boldsymbol{\theta}}_{k-1}$, and $\Delta_{a}$ is obtained recursively. For example, the second order difference function, $\Delta_{2}$, has the following form:

[TABLE]

We can write the difference functions $\Delta_{a}$ into a matrix form, \mbox{\boldmathL}_{a}. For example, when $a=1$ we have

[TABLE]

The positive definite penalty matrix used to control the smoothness is denoted as ${\bf R}_{2}$, and it has the quadratic form {\bf R}_{2}=\mbox{\boldmathL}_{a}^{\top}\mbox{\boldmathL}_{a}.

In the web-application NCSDE, we incorporated both penalties and either one can be used by the user.

2.4 Number of basis functions (clusters)

A critical step in dealing with real data is identifying the number of common basis functions (choice of $K$ in the context of NCSDE), which is directly related to the number of clusters, $\tilde{k}$. We identify the number of basis functions (clusters), using the elbow method as typically used in clustering analysis and it can be traced back to the work by Thorndike (1953). First, we obtain $\hbox{log}(\mathsf{S.Ps})={\bf B}\left({\bf B}^{\top}{\bf B}\right)^{-1}{\bf B}^{\top}\hbox{log}({\bf I})$, which is then used as an input to the elbow method, implemented in the R package factoextra (Kassambara and Mundt, 2016), to determine the number of clusters and basis functions.

In the elbow method, we use hierarchical clustering for partitioning the data, and obtain the total within-cluster sum of squares (WSS) based on different numbers of clusters $\tilde{k}$ (e.g., by varying $\tilde{k}$ from $1$ to $10$). We note that WSS measures the concentration of the clusters and it is desired to be as small as possible. Therefore, the optimal number of clusters can be obtained by plotting WSS against the number of clusters $\tilde{k}$. The location of an elbow (turning point) in this plot is generally considered an indicator of the appropriate number of clusters (see Figures 1(a), 5(a) and 6(a) for illustration purposes). From what we have seen in the implementation of NCSDE, a clear-cut elbow is achievable in longer time series ($n\geq 400$).

2.5 Identifiability of $(\Theta,\mbox{A})$

The non-uniqueness of the parametrization of $({\boldsymbol{\Theta}},{\bf A})$ causes an identifiability issue. Specifically, if ${\bf U}$ is a $K\times K$ orthogonal matrix, then ${\boldsymbol{\Theta}}{\boldsymbol{\alpha}}_{i}=({\boldsymbol{\Theta}}{\bf U})({\bf U}^{\top}{\boldsymbol{\alpha}}_{i})$. Thus, $\tilde{\boldsymbol{\Theta}}={\boldsymbol{\Theta}}{\bf U}$ and $\tilde{\boldsymbol{\alpha}}_{i}={\bf U}^{\top}{\boldsymbol{\alpha}}_{i}$ give the same representation of (5). To gain identifiability, we require that $(i)$ ${\boldsymbol{\Theta}}^{\top}{\boldsymbol{\Theta}}={\bf I}$, $(ii)$ ${\bf A}^{\top}{\bf A}={\bf D}^{2}$ is a diagonal matrix, $(iii)$ the columns of ${\bf A}$ are ordered such that the diagonal elements of ${\bf D}^{2}$ are in strictly decreasing order, and $(iv)$ the first non-zero element of each column of ${\boldsymbol{\Theta}}$ is positive. With such ${\boldsymbol{\Theta}}$ and ${\bf A}$, if the diagonal elements of ${\bf D}$ are all different, and we set $\bar{{\bf A}}={\bf A}{\bf D}^{-1}$ so that $\bar{\bf A}^{\top}\bar{\bf A}={\bf I}$, then we have ${\boldsymbol{\Theta}}{\bf A}^{\top}={\boldsymbol{\Theta}}{\bf D}\bar{{\bf A}}^{\top}$, which is a uniquely defined singular value decomposition (SVD). The desired identifiability of $({\boldsymbol{\Theta}},{\bf A})$ then follows from the uniqueness of the SVD.

3 Simulation study

We conducted a simulation study to evaluate the proposed nonparametric collective spectral density estimation method and compare it with non-collective spectral density estimation approaches. First, we describe the simulation setup and, in the next subsection, we list the competitive non-collective alternatives.

3.1 Simulation setup

In the simulations, we considered autoregressive models of order three, $AR(3)$. In each simulation run, we generated an $n\times m$ matrix ${\bf X}=({\bf X}_{1},{\bf X}_{2},\cdots,{\bf X}_{m})$, such that ${\bf X}_{i}^{\top}=(X_{i1},X_{i2},\cdots,X_{in})$. Each ${\bf X}_{i}$ was a time series of length $n$ from one of the following models:

[TABLE]

The $AR(3)$ model is given by $X_{t}=\sum_{i=1}^{3}\varphi_{i}X_{t-i}+\varepsilon_{t}$, with the associated SDF

[TABLE]

where $\sigma_{\varepsilon}^{2}=Var(\varepsilon_{t})$ (see von Storch and Zwiers, 2001, for details). We set $n=(100,200,400)$, $m=(6,15,30)$, and $p_{k}=\displaystyle{{1}\over{3}}$, for $k=1,2,3.$ Therefore, we had nine different pairs of $(n,m)$ and, for each pair we ran the simulation $N=100$ times.

3.2 Competitive approaches

In each simulation run we used the $n\times m$ data matrix ${\bf X}$ as an input to obtain estimates of $m$ SDFs ($\hat{\bf f}$: $\tilde{n}\times m$ matrix) from the following six methods:

• •

Periodograms (Ps):

We used the periodogram ${\bf I}$ as a rough estimate. For each run, the $i^{th}$ column of the $\tilde{n}\times m$ matrix ${\bf I}$ was a vector of size $\tilde{n}$ obtained via

[TABLE]

evaluated at the discrete frequencies ${\boldsymbol{\omega}}=\left(\omega_{1},\cdots,\omega_{\tilde{n}}\right).$

• •

Smoothed Periodograms (S.Ps):

The next estimate was obtained by smoothing ${\bf I}$ using the rich set of basis functions, \mathsf{S.Ps}=\exp\Bigl{[}{\bf B}\bigl{(}{\bf B}^{\top}{\bf B}\bigr{)}^{-1}{\bf B}^{\top}\hbox{log}({\bf I})\Bigr{]}.

• •

tSVD Periodograms (tSVD.Ps):

We used the truncated SVD to obtain the rank $K$ approximation of the smoothed periodograms $\left({\bf B}^{\top}{\bf B}\right)^{-1}{\bf B}^{\top}\hbox{log}({\bf I})$, called ${\boldsymbol{\Theta}}_{1}{\bf A}_{1}^{\top}$. This was referred to in Stewart (1993) as the approximation theorem (a.k.a. the Eckart-Young theorem). Therefore, the third estimate, “tSVD Periodograms”, is $\mathsf{tSVD.Ps}=\exp\left[{\bf B}{\boldsymbol{\Theta}}_{1}{\bf A}_{1}^{\top}\right]$.

• •

Separate Estimations (NSDE):

We avoided collective estimation and obtained the nonparametric spectral density estimates (NSDE) by maximizing the Whittle likelihood using the rich family of basis functions ${\bf B}$ separately. From (6), $u_{i}({\boldsymbol{\omega}})={\bf B}\boldsymbol{\psi}_{i}$. Therefore, our fourth estimate is $\mathsf{NSDE}=\exp\left[{\bf B}\mathbf{\Psi}\right]$, where $\mathbf{\Psi}=(\boldsymbol{\psi}_{1},\dots,\boldsymbol{\psi}_{m})$.

• •

tSVD Separate Estimations (tSVD.NSDE):

Having obtained $\mathbf{\Psi}$ by maximizing the Whittle likelihood separately, we used the truncated SVD to obtain the rank $K$ approximation of $\mathbf{\Psi}$ called ${\boldsymbol{\Theta}}_{s}{\bf A}_{s}^{\top}$. The fifth estimate is $\mathsf{tSVD.NSDE}=\exp\left[{\bf B}{\boldsymbol{\Theta}}_{s}{\bf A}_{s}^{\top}\right]$.

• •

Collective Estimations (NCSDE):

Finally, we used the proposed method from Section 2 to obtain $\mathsf{NCSDE}=\exp\left[{\bf B}\widehat{{\boldsymbol{\Theta}}}\widehat{{\bf A}}^{\top}\right]$.

3.3 Measures of quality

3.3.1 Adjusted Rand similarity coefficient

The adjusted Rand similarity coefficient (adjusted Rand index) was used to compare the (dis-)similarity between two clustering results. It is defined as follows (Vinh et al., 2009):

[TABLE]

To calculate the adjusted Rand index, we computed the $2\times 2$ contingency table, which consists of the following four cells:

• •

$n_{11}$: the number of observation pairs where both observations are comembers in both clusterings.

• •

$n_{10}$: the number of observation pairs where the observations are comembers in the first clustering but not the second.

• •

$n_{01}$: the number of observation pairs where the observations are comembers in the second clustering but not the first.

• •

$n_{00}$: the number of observation pairs where no pairs are comembers in either clustering result.

Furthermore $n_{i\cdot}$ and $n_{\cdot j}$ are defined as $n_{i0}+n_{i1}$ and $n_{0j}+n_{1j}$, respectively. Note that the $\mathsf{ARI}$ ranges from [math] to $1$, with [math] indicating that the two clusters do not agree on any pairs and $1$ indicating that the clusters are exactly the same. In subsection 3.4, we present and compare the clustering results from the approaches given in subsection 3.2, considering the original labels from the simulation as the gold standard. We did not use the class labels when applying the clustering algorithms; we only use the class labels to evaluate the clustering results.

3.3.2 Canonical angle

The canonical angle between the column spaces of ${\bf u}=\hbox{log}(\bf f)$ and $\hat{{\bf u}}=\hbox{log}(\hat{\bf f})$ is defined as the maximum angle between any two vectors from the two spaces. Mathematically, it can be computed as $\mathsf{angle}=\cos^{-1}(\rho)\times 180/\pi$, where $\rho$ is the minimum singular value of the matrix $\mathbf{Q}_{\hat{{\bf u}}}^{\top}\mathbf{Q}_{{\bf u}}$, and $\mathbf{Q}_{\hat{{\bf u}}}$ and $\mathbf{Q}_{{\bf u}}$ are orthonormal matrices obtained by the QR decomposition of matrices $\hat{{\bf u}}$ and ${\bf u}$, respectively (Golub and Van Loan, 2013). Smaller values of the angle indicate better estimations of the SDF weight functions.

3.4 Simulation results

In each simulation run, we used the approaches given in subsection 3.2 for nine different pairs of ($m,n$) to obtain matrices of size $\tilde{n}\times m$ that represent the associated estimates of $m$ SDFs for each approach. Figure 1(a) illustrates how to use the elbow method to obtain the number of clusters ($\tilde{k}=3$), and the remaining subfigures (Figure 1(g)-(f)) provide the results of the five approaches (S.Ps, tSVD.Ps, NSDE, tSVD.NSDE, and NCSDE) when estimating the SDFs for a randomly selected simulation run with $m=30$ and $n=400$. NCSDE obtained the results that are the smoothest and the closest to the true SDFs. In order to further verify the effectiveness of NCSDE technique, we used the canonical angle (details are given in the previous subsection) as a measure of closeness. Table 1 shows that $(i)$ the dimension reduction technique clearly improved the efficiency of the results (tSVD.Ps is better than Ps, and tSVD.NSDE is better than NSDE), and $(ii)$ incorporating the dimension reduction technique within the iterative procedure of estimating the SDFs (NCSDE) outperformed the sequential approach of estimating the SDFs separately and then proceeding with the dimension reduction technique (tSVD.NSDE).

The distance matrices obtained from the approaches given in subsection 3.2 can be used as an input to the hierarchical clustering algorithm. In our shiny App, we implemented the hclust function with option {method=”ward.D2”} in the R package stats to obtain dendrograms (R Core Team, 2016) (see Figure 2 for a comparison of four of the approaches, Ps, tSVD.Ps, NSDE, and NCSDE, on a randomly selected simulation run with $m=30$ and $n=400$). Since we knew the number of clusters and the cluster labels (gold standard) in advance, we cut the dendrograms into the original number of clusters.

To evaluate the performance of the proposed method in discovering the correct label (gold standard), in a more systematic framework, we used the popular external measure, adjusted Rand index, which is commonly used in the clustering evaluation literature and is discussed in subsection 3.3. The associated results are given in Table 2. NCSDE clearly outperformed the competing approaches in this clustering task, based on the adjusted Rand index similarity ($\mathsf{ARI}$).

We proceeded with a two-step procedure to visualize the level of similarity among the SDFs:

• i)

Obtain the $\tilde{n}\times m$ matrix ($\hat{\bf f}$) from any of the approaches given; then

• ii)

Use a multidimensional scaling (MDS) technique (Borg and Groenen, 2005).

However, another advantage of the proposed method is that the coefficients of the basis expansion for the fitted densities (${\bf A}$ matrix) provide a low-dimensional representation that can be directly used to visualize and cluster the spectral densities. Figure 3 shows the scatter plot of the first three coefficients (${\bf A}_{.1}$, ${\bf A}_{.2}$, and ${\bf A}_{.3}$) for the last two approaches in subsection 3.2 (tSVD.NSDE and NCSDE). We observe a clear segregation into three classes using NCSDE, indicating that these coefficients are useful for clustering purposes.

4 EEG data application

4.1 EEG data introduction

The EEG data was collected during the resting state from a single subject, a male student randomly selected from the 17 subjects in the experiment, which is described in detail in Wu et al. (2014). During the experiment, the EEG signals were recorded from 256 channels on the scalp surface, with a millisecond resolution (1000 recordings per second). From the 256 channels, 62 channels were eliminated due to issues with the quality of the data (Wu et al., 2014). The locations of the channels are shown in Figure 4. In this application, we cluster the remaining 194 channels using the time series data from the first and the last minutes. The total length of the time series is 60,000. However, we only consider the low-frequency bands (i.e., we truncate the perdiodograms at length $3000$).

4.2 Clustering results

The elbow method (Figures 5(a) and 6(a) suggested the presence of 4 different clusters in both cases, the first and the last minutes. Therefore, we set $\tilde{k}=4$ for both cases and proceeded with the NCSDE iterative procedure to obtain the Whittle maximum likelihood estimator.

For the data from the first minute, the algorithm converged after $33$ iterations. Figure 5(b) provides the result of the hierarchical clustering algorithm based on the estimated coefficients of the basis expansion, $\widehat{{\bf A}}$ matrix, where the associated tree has been cut into $\tilde{k}=4$ clusters. Out of the $194$ channels, $57,\ 42,\ 46$, and $49$ SDFs were assigned respectively to the four clusters. Figure 5(c) shows the $194$ estimated SDFs associated with the four different clusters. Figure 5(d) presents the clustering results in the 2D brain map.

The results from the last minute of data are presented in Figure 6. The algorithm converged after $170$ iterations, and Figure 6(b) provides the hierarchical clustering result based on the estimated $\widehat{{\bf A}}$ matrix for the last minute of data. By cutting the tree into $\tilde{k}=4$ clusters, we determined that $22,\ 75,\ 54$, and $43$ SDFs should be assigned to the first through fourth clusters, respectively. Figure 6(c) shows the functional means of the estimated SDFs in each cluster to illustrate the difference between the four groups. Figure 6(d) presents the clustering results in the 2D brain map.

As shown in Figure 5 and Figure 6, we have the following findings:

For both datasets (first and last minutes), the elbow method suggested using four clusters, and both of the brain maps provide relatively well-separated, symmetric regions from the clustering results. 2. 2.

All of the brain maps give three regions: the front, middle , and back regions. From the functional mean of each cluster, we find that the channels in the middle region of the brain have a lower density in the low-frequency band than the front and back regions. 3. 3.

Both data are associated with the resting state; when we compare the results from the first and last minutes, we find that there are no significant differences between them, except in the front region where we obtain two clusters in first minute (Figure 5(d)), but just one cluster in last minute (Figure 6(d)). 4. 4.

The most left-front channels in Figure 6(d) have a different pattern from the other channels.

5 Conclusion

A novel approach for collectively estimating multiple SDFs was developed in this paper. By pooling data from different time series in the frequency domain and using a shared basis to represent the SDFs, the collective estimation approach is statistically more efficient than non-collective estimation approaches. The proposed method uses the penalized Whittle likelihood approximation to yield a flexible family of spectral densities. As an output of the new method, each estimated log-spectral density is expressed in a basis expansion where the basis is estimated from the data, assuming that the SDFs lie in a low-dimensional manifold of the large space spanned by a pre-specified rich basis. The collective spectral density estimation approach is widely applicable when there is a need to estimate multiple SDFs from different populations. Moreover, the coefficients of the basis expansion for the fitted spectral densities provide a concise, low-dimensional representation that could be useful for visualization and clustering. Another advantage of the new procedure is that it speeds up the process by updating the smoothness parameter within the Newton-Raphson iterations and avoids a grid search over the space of the smoothing parameter, $\lambda$, which could be very time consuming.

A web application that can be used by the research community to reproduce the results in this paper or to estimate SDFs collectively based on NCSDE for any other related applications is available at “https://ncsde.shinyapps.io/NCSDE.”