Simultaneous clustering and estimation of additive shape invariant models for recurrent event data

Consider a set of recurrent events from repeated measurements of $n$ subjects $\mathcal{O}\equiv\big{\{}\left\{t_{i,r,j}\right\}_{j=1,\ldots,N_{i,r}(T)}:0<t_{i,r,1}<\cdots<t_{i,r,N_{i,r}(T)}<T,i=1,\ldots,n,r=1,\ldots,R\big{\}},$ where $t_{i,r,j}$ denotes the time of the $j$-th event of the $i$-th subject in the $r$-th observation, $T$ denotes the duration of each observation, $N_{i,r}(T)$ denotes the total number of events associated with the $i$-th subject and the $r$-th observation. For $i\in[n]$ and $r\in[R]$, we adopt the definition of counting process (see, e.g. Daley and Vere-Jones, 2003) as $N_{i,r}(t)\equiv\sum_{j=1}^{N_{i,r}(T)}\textbf{1}(t_{i,r,j}\leq t)$ for $t\in[0,T]$. We use the intensity function, i.e., $\lambda_{i,r}(t)\equiv\mathbb{E}\{\mathrm{d}N_{i,r}(t)/\mathrm{d}t\}$, to characterize each counting process.

We start with a model that tackles the first two of the aforementioned tasks: decomposition and alignment. To be specific, we assume that the true underlying intensities are formed through the superposition of multiple components, and these intensity components exhibit uniformity across subjects, differing only in temporal shifts. We arrive at the following model

where $a\in[0,\infty)$ represents the baseline intensity, $M\in\mathbb{N}^{+}$ is the number of components, $v_{i,m}\in[0,V]$ represents the subject-specific time shift associated with the $i$-th subject and the $m$-th component, $w^{*}_{r,m}\in[0,W]$ represents the known observation-specific time shift of the $m$-th component in the $r$-th observation, $S^{v}$ is the shift operator defined as $S^{v}x(t)\equiv x(t-v)$, and $f_{m}(t)\in\mathcal{F}$ represents the $m$-th intensity component. Here $V,W\in(0,T)$, and $\mathcal{F}\equiv\{f\in L^{2}(\mathbb{R}):f(t)=0\text{ for }t\in\mathbb{R}\setminus(0,T_{0})\}$ with $T_{0}\in(0,T)$. For further clarity, a graphical representation of the model in (1) is provided in Figure 2.

To provide a concrete example for Model 1, consider the experiment in Steinmetz et al. (2019) where $N_{i,r}(t)$ represents the recorded spike train of the $i$-th neuron in the $r$-th trial, and $\lambda_{i,r}(t)$ denotes the corresponding firing rate (see for instance Figure 1). The terms $a$ and $\{f_{m}(\cdot):m=1,2\}$ can be interpreted as the spontaneous firing rate and the neural response elicited by the visual gratings and auditory cue, respectively. The subject-specific time shift $v_{i,m}$ corresponds to the response latency of the $i$-th neuron in response to the $m$-th stimulus. The observation-specific time shift $w^{*}_{r,m}$ corresponds to the occurrence time of the $m$-th stimulus in the $r$-th trial, which is randomly generated by the experimenters.

We denote the collection of model parameters as $\boldsymbol{\theta}_{0}\equiv({a},\mathbf{f},\mathbf{v})\in\Theta_{0}$, where $\mathbf{f}\equiv(f_{m})_{m\in[M]}$, $\mathbf{v}\equiv(v_{i,m})_{i\in[n],m\in[M]}$, and $\Theta_{0}\equiv\{({a},\mathbf{f},\mathbf{v}):{a}\in[0,\infty),\mathbf{f}\in\mathcal{F}^{M},\mathbf{v}\in[0,V]^{n\times M}\}$. In addition, for any $N_{i,r}(t)$, we denote its conditional intensity given the observation-specific time shifts as $\lambda_{\boldsymbol{\theta}_{0},i}(t,\mathbf{w}^{*})\equiv\mathbb{E}_{\boldsymbol{\theta}_{0}}\{\mathrm{d}N_{i,r}(t)/\mathrm{d}t\mid\mathbf{w}^{*}_{r}=\mathbf{w}^{*}\}$, where $\mathbf{w}^{*}_{r}\equiv(w^{*}_{r,m})_{m\in[M]}$. We further denote $\boldsymbol{\lambda}_{\boldsymbol{\theta}_{0}}(t,\mathbf{w}^{*})\equiv(\lambda_{\boldsymbol{\theta}_{0},i}(t,\mathbf{w}^{*}))_{i\in[n]}$. In this context, identifiability pertains to the injectivity of the mapping $\boldsymbol{\theta}_{0}\mapsto\boldsymbol{\lambda}_{\boldsymbol{\theta}_{0}}(\cdot,\cdot)$ for $\boldsymbol{\theta}_{0}\in\Theta_{0}$.

With these notations, we can formally present the identifiability results in Proposition 1.

Assumption A1 posits that the observation duration sufficiently extends to avoid censorship. Given A1, the model described in (1) can be formulated in the form of a linear regression model in the frequency domain, where the response variables are determined by $N_{i,r}$’s, the explanatory variables depend on $w^{*}_{r,m}$’s, and the regression coefficients are functions of $f_{m}$’s and $v_{i,m}$’s. Within the context of linear regression, the significance of Assumption A2 becomes straightforward. Specifically, Assumption A2 effectively assumes no collinearity among explanatory variables. We note that Assumption A2 is satisfied if either (i) the observation-specific time shifts are randomized, i.e., $\{w^{*}_{m}\in\mathbb{C}:m=1,2,\ldots M\}$ or (ii) the gaps between the observation-specific time shifts are randomized, i.e., $w^{*}_{m}=w^{*}_{m-1}+\delta_{m-1}$ for $m=2,\ldots,M$, where $\{\delta_{m}\in\mathbb{C}:m\in[M-1]\}$ are independent random variables with non-zero variance. Detailed proof of Proposition 1 is presented in Section A of the Supplementary Materials.

In special scenarios, the model presented in (1) has connections with three distinct branches of research. Firstly, considering a special case where $M=1$ and $R=1$, the model in (1) can be simplified as follows:

where $g(t)\equiv a+S^{w^{*}_{1,1}}f_{1}(t)$. This model has been investigated in the realm of shape invariant models (see, e.g., Beath, 2007; Bigot and Gadat, 2010; Vimond, 2010; Bigot et al., 2013; Bigot and Gendre, 2013; Bontemps and Gadat, 2014). Consequently, the proposed model can be seen as a generalization of shape invariant models to incorporate additive components. Various assumptions have been proposed to ensure the identifiability of shape invariant models, including availability of repeated observations (Bigot and Gendre, 2013), partial knowledge of time shifts (Bigot and Gadat, 2010; Bigot et al., 2013), or the knowledge of the noise distribution (Bontemps and Gadat, 2014). For the proposed additive shape invariant model, where $M$ is allowed to exceed 2, a combination of repeated observations and partial knowledge of time shifts suffices to ensure the model identifiability, as elucidated in Proposition 1.

Secondly, considering the scenario where the intensity components $\{S^{v_{i,m}+w^{*}_{r,m}}f_{m}(t):m\in[M]\}$ have non-overlapping supports for all $i\in[n]$ and $r\in[R]$, model (1) can be expressed as:

where $g(s)\equiv a+\sum_{m\in[M]}f_{m}(s)$ and $h_{i,r}(t)\equiv\sum_{m\in[M]}\{t-(v_{i,m}+w^{*}_{r,m})\}\times\mathbf{1}[f_{m}\{t-(v_{i,m}+w^{*}_{r,m})\}\neq 0]$. The reformulated model in (3) means that $\lambda_{i,r}$’s are identical subject to unknown time warping functions (i.e., $h_{i,r}(t)$’s). The task of estimating time warping functions has been explored in the domain of curve registration (Kneip and Gasser, 1992; Ramsay and Li, 1998; James, 2007; Telesca and Inoue, 2008; Cheng et al., 2016).

Lastly, considering the scenario where the variances of $v_{i,m}$’s and $w^{*}_{r,m}$’s are both close to zero, model (1) can be approximated using Taylor expansion as:

where $\mu(t)\equiv a+\sum_{m\in[M]}f_{m}(t-\mathbb{E}u_{i,r,m})$, $u_{i,r,m}\equiv v_{i,m}+w^{*}_{r,m}$, $\zeta_{i,r,m}\equiv-(u_{i,r,m}-\mathbb{E}u_{i,r,m})\|Df_{m}(t-\mathbb{E}u_{i,r,m})\|_{t}$, $Df_{m}$ denotes the first order derivative of $f_{m}$, and $\psi_{m}(t)\equiv Df_{m}(t-\mathbb{E}u_{i,r,m})\|Df_{m}(t-\mathbb{E}u_{i,r,m})\|_{t}^{-1}$. We defer the derivation of this approximation to Section C of the Supplementary Materials. When functions $\{S^{\mathbb{E}u_{i,r,m}}f_{m}(t):m\in[M]\}$ exhibit non-overlapping supports, the approximate model in (4) corresponds to the models of functional principal component analysis (FPCA) (Yao et al., 2005; Morris and Carroll, 2006; Di et al., 2009; Crainiceanu et al., 2009; Xu et al., 2020).

We consider the case that $N_{i,r}(t)$’s are Poisson processes. To effectively estimate the parameters in (1), we exploit the two sources of stochasticity of the Poisson processes. Firstly, the number of events over $[0,T]$ for any Poisson process is a random variable that follows a Poisson distribution, in other words, $N_{i,r}(T)\sim\mathrm{Poisson}(\Lambda_{i,r})$, where $\Lambda_{i,r}\equiv\int_{0}^{T}\lambda_{i,r}(t)\mathrm{d}t$ represents the expected event count. When Assumption A1 holds, we can derive from the model in (1) that $\Lambda_{i,r}=aT+\sum_{m\in[M]}\int_{0}^{T}f_{m}(t)\mathrm{d}t$, denoted as $\Lambda$ for simplicity. Secondly, conditioning $N_{i,r}(T)$, the event times $\{t_{i,r,j}:j=1,\ldots,N_{i,r}(T)\}$ can be regarded as independent and identically distributed random variables. The probability density function of these event times, or event time distribution, is characterized by $\lambda_{i,r}(t)\Lambda^{-1}=a\Lambda^{-1}+\sum_{m\in[M]}S^{v_{i,m}+w^{*}_{r,m}}f_{m}(t)\Lambda^{-1}$.

Accordingly, we estimate $a,\mathbf{f},\mathbf{v}$ through the following reparameterization. Letting $a^{\prime}\equiv a\Lambda^{-1}$and $\mathbf{f}^{\prime}\equiv(f_{m}\Lambda^{-1})_{m\in[M]}$, we have the following optimization problem:

where $\|\cdot\|_{t}$ denotes the $L^{2}$-norm with respect to $t$, $\beta_{i,r}\equiv N_{i,r}(T)$, $y_{i,r}(t)\equiv\{N_{i,r}(t+\Delta t)-N_{i,r}(t)\}\Delta t^{-1}$, and $\Delta t$ represents a infinitesimally small value. In (5), the objective function $L_{1}({a}^{\prime},\mathbf{f}^{\prime},\mathbf{v})$ measures the discrepancy between the empirical and estimated distributions of event timings, where $y_{i,r}(t)N_{i,r}(T)^{-1}$ and $\{a^{\prime}+\sum_{m\in[M]}S^{v_{i,m}+w^{*}_{r,m}}f^{\prime}_{m}(t)\}$ serve as the empirical distribution and the estimated distribution of $\{t_{i,r,j}:j=1,\ldots,N_{i,r}(T)\}$, respectively. The term $\beta_{i,r}$ serves as the weight of the counting process $N_{i,r}(t)$. For instance, setting $\beta_{i,r}=N_{i,r}(T)$ means equal weights for all events. Finally, we estimate $\Lambda$ using the empirical mean

The parameters $a,\mathbf{f}$ can be estimated from $\hat{\Lambda},\hat{a}^{\prime},\hat{\mathbf{f}}^{\prime}$ by $\hat{a}\equiv\hat{a}^{\prime}\hat{\Lambda}$ and $\hat{\mathbf{f}}\equiv\hat{\mathbf{f}}^{\prime}\hat{\Lambda}$.

We extend model (1) to simultaneously perform decomposition, alignment, and clustering. Assume that the $n$ subjects can be classified into $K$ distinct clusters, that is, $[n]=\cup_{k=1}^{K}\mathcal{C}_{k}$, where $\mathcal{C}_{1},\ldots,\mathcal{C}_{K}$ represent mutually exclusive subsets of $[n]$. These clusters are delineated based on the similarity of intensity components across subjects. Specifically, we introduce the following model:

where $z_{i}\in[K]$ represents the cluster membership of subject $i$ such that $z_{i}=k$ if $i\in\mathcal{C}_{k}$. We refer to the model in (7) as the additive shape invariant mixture model, or ASIMM for short. Conditioning on each cluster, the additive shape invariant mixture model in (7) simplifies to the additive shape invariant model in (1). Similar to the connection between model in (1) and FPCA, the additive shape invariant mixture model in (7) has a close connection with clustering methods based on FPCA (Chiou and Li, 2007; Bouveyron and Jacques, 2011; Jacques and Preda, 2013; Yin et al., 2021).

In the context of neural data (Steinmetz et al., 2019), $\mathcal{C}_{k}$’s can be interpreted as functional groups of neurons, wherein neurons exhibit similar firing patterns. The additive shape invariant mixture model in (7) enables us to simultaneously identify functional groups of neurons (i.e., $\mathcal{C}_{k}$’s), discern representative neural firing patterns (i.e., $f_{k,m}$’s), and estimate individual neural response latencies (i.e., $v_{i,m}$’s). It is worthwhile to emphasize that the applicability of the proposed method extends beyond neural data analysis. For instance, the additive shape invariant mixture model in (7) can be employed in analyzing recurrent consumer activity in response to advertising (Xu et al., 2014; Zadeh and Sharda, 2014; Tanaka et al., 2016; Bues et al., 2017), or studying hospital admission rates following the implementation of societal disease prevention policies (Barone-Adesi et al., 2006; Sims et al., 2010; Klevens et al., 2016; Evans et al., 2021). Additionally, the proposed model has potential for application on diverse datasets (Tang and Li, 2023; Schoenberg, 2023; Dempsey, 2023; Ganggang Xu and Guan, 2024; Djorno and Crawford, 2024).

We denote the collection of unknown parameters in (7) as $\mathbf{z}\equiv(z_{i})_{i\in[n]}$, $\mathbf{a}\equiv(a_{k})_{k\in[K]}$, $\mathbf{f}\equiv(f_{k,m})_{k\in[K],m\in[M]}$, $\mathbf{v}\equiv(v_{i,m})_{i\in[n],m\in[M]}$. Let $(\mathbf{z}^{*},\boldsymbol{a}^{*},\mathbf{f}^{*},\mathbf{v}^{*})\in\Theta_{1}$ denote the true parameters, where $\Theta_{1}\equiv\{(\mathbf{z},\boldsymbol{a},\mathbf{f},\mathbf{v}):\mathbf{z}\in[K]^{n},\boldsymbol{a}\in[0,\infty)^{K},\mathbf{f}\in\mathcal{F}^{K\times M},\mathbf{v}\in[0,V]^{n\times M}\}$. When conditioning on $\mathbf{z}^{*}$, Proposition 1 establishes the identifiability of $\boldsymbol{a}^{*},\mathbf{f}^{*},\mathbf{v}^{*}$. However, to ensure the identifiability of $\mathbf{z}^{*}$, an additional assumption regarding the separability of clusters is required. The formal presentation of model identifiability is provided in Proposition 2.

Assumption A3 mandates that each cluster exhibits at least one signature intensity component that is unique to this cluster. Statement P4 directly stems from Assumption A3. Statements P5, P6, and P7 can be derived by applying Proposition 1 to each individual cluster. Detailed proof of Proposition 2 is presented in Section B of the Supplementary Materials.

We estimate the parameters in (7) by generalizing the optimization approach in Section 2.4 to incorporate the cluster structure. For any $k\in[K]$ and $m\in[M]$, let $\Lambda_{k}\equiv a_{k}T+\sum_{m\in[M]}\int_{0}^{T}f_{k,m}(t)\mathrm{d}t$, $a^{\prime}_{k}\equiv a_{k}\Lambda_{k}^{-1}$, and $f^{\prime}_{k,m}\equiv f_{k,m}\Lambda_{k}^{-1}$. Denoting $\mathbf{a}^{\prime}\equiv(a^{\prime}_{k})_{k\in[K]}$, $\mathbf{f}^{\prime}\equiv(f^{\prime}_{k,m})_{k\in[K],m\in[M]}$, and $\boldsymbol{\Lambda}\equiv(\Lambda_{k})_{k\in[K]}$, we propose the following optimization problem:

where

and $\gamma\in(0,\infty)$ is a tuning parameter. In essence, $L_{1}$ and $L_{2}$ assess the within-cluster variance of event time distributions and event counts, respectively. When the number of clusters is reduced to one (i.e., $K=1$), the definitions of $L_{1}$ in (9) is identical to the definition in (5). The tuning parameter $\gamma$ modulates the relative importance of $L_{2}$ compared to $L_{1}$ in the optimization with respect to $\mathbf{z}$. When $\gamma$ is sufficiently small, the estimator $\hat{\mathbf{z}}$ defined in (8) is predominantly determined by $L_{1}$, resulting in a potentially suboptimal value of $L_{2}$. Conversely, when $\gamma$ is sufficiently large, the dominance shifts towards $L_{2}$, resulting in a $\hat{\mathbf{z}}$ that achieves the minimum value of $L_{2}$, while $L_{1}$ may be relegated to suboptimal values. Subsequent to addressing the optimization problem in (8), the estimations of $\mathbf{a}$ and $\mathbf{f}$ can be established via $\hat{a}_{k}\equiv\hat{a}^{\prime}_{k}\hat{\Lambda}_{k}$ and $\hat{\mathbf{f}}_{k,m}\equiv\hat{\mathbf{f}}^{\prime}_{k,m}\hat{\Lambda}_{k}$ for $k\in[K]$, $m\in[M]$.

The optimization problem in (8) involves two tuning parameters $\gamma$ and $K$. To determine these tuning parameters, we employ a heuristic method. We first establish a preliminary estimate of $K$ using simple methods, such as applying the k-means algorithm on $N_{i,r}(T)$’s and selecting $K$ using the elbow method (Thorndike, 1953). Given this preliminary estimation of $K$, we choose the largest $\gamma$ before observing a significant increase in $L_{1}$. We provide simulation experiments to justify this heuristic method and demonstrate the robustness of selected $\gamma$ to the change of the preliminary selection of $K$ in Section F.1 in the Supplementary Material. Finally, conditioning on the selected $\gamma$, we refine the value of $K$ by identifying the elbow point on the curve of the overall objective function in (8) against $K$.

The centering step (11) involves two optimization problems. To solve the first optimization problem, we formulate it in the frequency domain as follows:

where $\boldsymbol{\phi}^{\prime}\equiv(\phi^{\prime}_{k,m,l})_{k\in[K],m\in[M],l\in\mathbb{Z}}$, $\{\phi^{\prime}_{k,m,l}:l\in\mathbb{Z}\}$ denotes the Fourier coefficients of $f^{\prime}_{k,m}(t)$, $\{\eta_{i,r,l}:l\in\mathbb{Z}\}$ denotes the Fourier coefficients of $y_{i,r}(t)$, and $\operatorname{j}$ denotes the imaginary unit. Notably, the multiplicative term $\exp(-\operatorname{j}2\pi l[\hat{v}_{i,m}+w^{*}_{r,m}]T^{-1})$ serves as the frequency domain counterpart of the shift operator $S^{\hat{v}_{i,m}+w^{*}_{r,m}}$. In essence, the Fourier transformation converts the shift operators into multiplication. As a result, the objective function becomes more tractable compared to its counterpart in the original domain. Indeed, an analytical solution for $\hat{\boldsymbol{\phi}}^{\prime}$ can be derived. Let $\boldsymbol{\phi}^{\prime}_{k,*,l}\equiv({\phi}^{\prime}_{k,m,l})_{m\in[M]}$ for any $k\in[K]$ and $l\in\mathbb{Z}$. For $l\neq 0$, the solution to (13) with respect to $\boldsymbol{\phi}^{\prime}_{k,*,l}$ can be expressed as follows:

Here, $\mathbf{E}_{k,l}$ is defined as

where $\hat{\mathcal{C}}_{k}\equiv\{i\in[n]:\hat{z}_{i}=k\}$. $\overline{\mathbf{E}_{k,l}}$ denotes the complex conjugate of $\mathbf{E}_{k,l}$, $\mathbf{B}_{k}$ is a diagonal matrix of $(\beta_{i,r})_{(i,r)\in\hat{\mathcal{C}}_{k}\times[R]}$, $\mathbf{h}_{k,l}\equiv(\eta_{i,r,l}N_{i,r}(T)^{-1})_{(i,r)\in\hat{\mathcal{C}}_{k}\times[R]}$, and $\ell_{0}\in\mathbb{N}$ is a truncation parameter introduced to ensure the numerical feasibility of computing $\hat{\boldsymbol{\phi}}^{\prime}_{k,*,0}$. For $l\neq 0$, the objective function concerning $\boldsymbol{\phi}^{\prime}_{k,*,l}$ in (13) is a weighted sum of squares, hence the estimate in (14) can be derived using the well-known least squares estimator. For $l=0$, the estimate in (15) can be obtained by exploiting the definition of $\mathcal{F}$. The detailed derivations of (14) and (15) can be found in Section D of the Supplementary Material.

Upon obtaining $\hat{\boldsymbol{\phi}}^{\prime}$, the solution to the first optimization problem in (11) can be derived as follows. For $k\in[K]$ and $m\in[M]$,

where (79) is obtained by substituting $\hat{\boldsymbol{\phi}}^{\prime}$ into (13), and (18) follows from the inverse Fourier transformation.

The second optimization problem in (11) aims to minimize the within-cluster variances of event counts given cluster memberships. The solution to this optimization problem is straightforward: for any $k\in[K]$,

where $|\hat{\mathcal{C}}_{k}|$ denotes the cardinality of the set $\hat{\mathcal{C}}_{k}$. In summary, the solution to the centering step is encapsulated by equations (79), (18), and (19).

The optimization problem in (12) can be scaled down to the subject level, allowing for the independent estimation of parameters associated with each subject. For any subject $i\in[n]$, let $\mathbf{v}_{i}\equiv(v_{i,m})_{m\in[M]}$ denote its associated time shifts. In (12), the parameters $z_{i}$ and $\mathbf{v}_{i}$ are estimated through the following sub-problem:

where $L_{1,i}$ and $L_{2,i}$ are defined as

The parameter dimension for the problem in (20) is significantly reduced compared to the original optimization problem in (12). Consequently, solving the problem in (20) is computationally more efficient than addressing the original problem stated in (12).

To solve the optimization problem in (20), we employ the following procedure:

In the first step (23), we determine the optimal time shift for each potential cluster membership $k\in[K]$. The optimization problem in this step can be solved in the frequency domain using the Newton’s method (see Section E of the Supplementary Material). In the second step (24), we evaluate the objective function for each possible cluster membership, leveraging the optimal time shift corresponding to that particular cluster. Subsequently, we designate the cluster associated with the minimal objective function value as the estimated cluster membership. In the last step (25), we choose the optimal time shift for the estimated cluster membership. In summary, the solution to the clustering step is encapsulated by (23), (24), and (25).

The objective function in (8) exhibits a non-convex nature. This characteristic poses a crucial need for an appropriate initialization scheme and convergence criterion. Our proposed initialization scheme is described in Remark 1. The overall estimation procedure is summarized in Algorithm 1.