Temporal Point Process Graphical Models

Graphical modeling of point processes is drawing increasing attention, as data in the form of point processes on a graph are emerging in scientific and business applications. Examples include biological neural networks with neural spike train data (see e.g. Paninski et al., 2007; Pillow et al., 2008; Bolding and Franks, 2018); social media data recording timestamps of each user’s actions (Perry and Wolfe, 2013; Zhou et al., 2013; Fox et al., 2016); high frequency financial data recording times of market orders (Bacry et al., 2013; Aït-Sahalia et al., 2015; Bacry et al., 2015). Understanding the connectivity structure of these point processes, i.e, the edges on the graphs, is a challenging problem of growing interest. Our work is motivated by a recent spike train data set from Bolding and Franks (2018) . This data set collected spike train data from two important brain regions for olfactory of mice, namely the olfactory blub and the piriform cortex. The neural activity in the olfactory bulb is highly concentration-dependent, while that in the piriform cortex is concentration-invariant. The study of Bolding and Franks (2018) suggests that there exists a complex excitation-inhabitation mechanism between the olfactory blub and the piriform cortex that allows the nervous system to produce a concentration-invariant representation of odors.

In this article, we propose a new temporal point process graphical model to address this problem. We model the spike train data of the neurons as a multivariate temporal point process, where each component of the point process represents the activity of a neuron. We allow future events be excited or inhabited by the past events in a non-additive manner. We organize the transfer coefficients into a matrix and impose a sparse structure, which is welcomed for interpretation in neuroscience and other applications (Zhao and Yu, 2006; Zou, 2006). Our loss function is convex, so many highly scalable algorithms, such as alternating direction method of multipliers and coordinate descent method, can be applied to parameter estimation. We establish a non-asymptotic error bound of the estimator and allow the number of nodes in the graph to diverge. We also analyze the approximation error introduced during the implementation phase.

There have been a great number of works focusing on graphical models, including the Gaussian graphical model (Yuan and Lin, 2007; Liu, 2013) and directed acyclic graph (Kalisch and Bühlman, 2007; Van de Geer et al., 2013; Zheng et al., 2018). In this models, each node of the graph is corresponding to a random variable and the data are realizations of the random variables. The graphs of classical graphical models are used to represent the structure of conditional independence among the random variables. In contrast, we aim at the case where each node of the graph represents a point process and the graph can be viewed as a multivariate point process. In our setting, the graph encodes the local independence relationship among the components of the multivariate point process (Didelez, 2008). The basic idea of local independence is that once we know about specific past events, the intensity of the future event is independent of other past events. In our model, local independence is fully characterised by the parameters (see Proposition 1).

There have been a family of models targeting temporal point processes. A Hawkes process (Hawkes, 1974) assumes that historical events can trigger future events and has been widely used in modelling neuronal spike train data and social network data. The model we proposed differs from existing research of the Hawkes process in many ways. Most prior art has focus on parameter estimation in the linear case, where the transfer function of the multivariate Hawkes process is the identity function. Hansen et al.(2015) considered estimating parameters in a linear Hawkes process of fixed number of nodes with least-square loss function and $\ell_{1}$ regularization. Bacry et al.(2020) organized the transferring coefficients as a matrix and impose a low-rank structure. A crucial weakness of linear Hawkes process is that it does not allow inhibitory relations, i.e, the past events can only trigger more future events, not silence future events. Our model extends to a nonlinear setting where inhibitory influences are allowed and the influence of past events from different nodes can accumulate in a non-additive manner.

Tang and Li(2021) also considered parameter estimation of a nonlinear Hawkes process in their work. Our work differs from Tang and Li(2021) in several aspects. First, we focus on the estimation of sparse graphical models while they focus on the estimation of a low-rank tensor of regression parameters. Second, we provide a stronger theoretical guarantee by proving the global minimum of our penalized loss function has a $\sqrt{\frac{s\log p}{T}}$ convergence rate while their theory only shows that there exists a local minimum of the penalized loss function in the neighborhood of the true parameter. Finally, our loss function is convex while theirs is not. A convex loss function with restricted strong convexity allows us to find global minimum efficiently.

The paper is organized as follows. We will finish this section with notation. We briefly introduce the temporal point process and our model in Section 2. In Section 3, we present the construction of our estimator and our main theorem. In Section 4, we analyze the approximation error introduced during the implementation phase. We show the results of simulations in Section 5 and an application on a spick train data set in Section 6. We conclude with a discussion in Section 7. Technical proofs are provided in the Appendix.

Some notations will be employed through out this paper. Let $[n],n\in\mathbb{N}$ denote the set $\{1,\ldots,n\}$. For a set $S\subset[p]$ and a vector $v\in\mathbb{R}^{p}$, we define $v_{S}\in\mathbb{R}^{p}$ as $(v_{S})_{i}=v_{i}\mathbb{I}(i\in S)$. Let $\|\cdot\|_{p},\ p\in[1,\infty)$ denote the $\ell_{p}$-norm of a vector and $\|\cdot\|_{0}$ denote the number of non-zero entries of a vector. For a matrix $A\in\mathbb{R}^{m\times n}$, let $\|A\|_{F}=(\sum_{i=1}^{m}\sum_{j=1}^{n}a_{ij}^{2})^{1/2}$ denote the Frobenius norm of a matrix. We use $S_{q}^{p}(r)$ for a $\|\cdot\|_{q}$ ball with radius $r$ in $\mathbb{R}^{p}$. We use $c_{1},c_{2},\ldots$ to denote positive numbers that do not depend on $p$ or $T$ but can still depend on some values that we consider as constants in this paper, and their value can change from line to line.

We will begin with a very brief review of temporal point processes. For a more systematic and comprehensive discussion of point processes, we refer to Brémaud (1981) and Daley and Vere-Jones (2007).

Let $N(t),t\geq 0$ be a simple point process. $N(t)$ is corresponding to a sequence of real-valued random variables $\{t_{k}\}_{k\in\mathbb{Z}}$, with $t_{k+1}>t_{k}$ and $t_{1}\geq 0$ almost surely. $\{t_{k}\}_{k\in\mathbb{Z}}$ are called the event times of $N(t)$. For a set $A$ in the Borel $\sigma$-field of the real line, $N(A)=\sum_{k}\mathbb{I}_{\{t_{k}\in A\}}$. We write $N([t,t+\mathrm{d}t))$ as $\mathrm{d}N(t)$, where $\mathrm{d}t$ denotes an arbitrarily small increment of $t$. Let $\mathcal{H}_{t}$ denote the $\sigma$-field generated by $N(\tau),\tau\in[0,t]$. The intensity of $N(t)$ is defined as

$$\lambda(t)\mathrm{d}t=\mathbb{P}(\mathrm{d}N(t)=1|\mathcal{H}_{t}).$$

Now suppose we have $p$ simple point processes $N_{1},\ldots,N_{p}$ with $t_{j,1}<t_{j,2}<\cdots<t_{j,n}<\cdots$ denoting the event times of the process $N_{j},j\in[p]$. Then we employ $\mathcal{N}=(N_{1},\ldots,N_{p})$ to denote a $p$-variate point process. Let $\mathcal{H}_{t}$ denote the $\sigma$-field generated by $\mathcal{N}(\tau),\tau\in[0,T]$. The intensity process $\bm{\lambda}(t)=(\lambda_{1}(t),\ldots,\lambda_{p}(t))^{\mathrm{T}}$ is a $p$-variate $\mathcal{H}_{t}$-predictable process defined as

$$\lambda_{j}(t)\text{d}t=\mathbb{P}\big{(}\text{d}N_{j}(t)=1|\mathcal{H}_{t}\big{)},\ j\in[p].$$

In our model, we assume this intensity process takes the form

$$\lambda_{j}(t)=h_{j}\Big{(}\mu_{j}+\sum_{k=1}^{p}\int_{0}^{t}\omega_{j,k}(\tau)\text{d}N_{k}(t-\tau)\Big{)},\ j\in[p],$$

(1)

where $h_{j}(\cdot),j\in[p]$ are link functions and can be nonlinear. $\mu_{j},j\in[p]$ are the background intensities of the components of the multivariate point process. $\omega_{j,k}(\cdot):\mathbb{R}^{+}\rightarrow\mathbb{R},j,k\in[p]$ are the transfer functions. $\bm{\lambda}(t)=\big{(}\lambda_{1}(t),\cdots,\lambda_{p}(t)\big{)}$ is called the intensity process. The convolutional expression $\int_{0}^{t}\omega_{j,k}(\tau)\text{d}N_{k}(t-\tau)$ summaries the historical information of the $k$-th component. In this paper, we consider parameterized transfer functions that can be written as

$$\omega_{j,k}(t)=\beta_{j,k}\kappa_{j,k}(t),\ j,k\in[p],$$

where $\kappa_{j,k}(t),\ j,k\in[p]$ are known kernel functions.

In graphical models, the edges encodes conditional independence structures among the random variables (Maathuis et al., 2018; Qiao et al., 2019; Engelke and Hitz, 2020). For multivariate point processes, the question of interest is that whether the intensity of one component depends on the history of another component, which can be represent by the local independence for multivariate point process (Didelez, 2008).

In temporal point process graphical models, the local independence structures are encoded in the transfer matrix $B=(\beta_{j,k})\in\mathbb{R}^{p\times p}$, because $\lambda_{j}(t)=h_{j}\big{(}\mu_{j}+\sum_{k=1}^{p}\int_{0}^{t}\beta_{j,k}\kappa_{j,k}(\tau)\text{d}N_{k}(t-\tau)\big{)}\in\mathcal{F}_{t}^{A}$ for index set $A$ as long as $\{k,\beta_{j,k}\neq 0\}\subset A$. We have the following proposition:

Thus by estimating the transfer matrix $B$, we can recover the local independence relationships in $\mathcal{N}(t)$.

In this section, we present an approach for estimating the temporal point process graphical models. We will begin with the construction of our estimator in Section 3.1, then provide a detailed comparison among our method with the widely employed method in prior arts, the least-square estimation and the maximum likelihood estimation in Section 3.2. In Section 3.3, we will present theoretical analysis that guarantees the performance of our method.

In this section, we investigate the performance of our estimation procedure in a simulation study. We propose two estimators, the naively weighted estimator with $W_{j}(t)=1,\ \forall j\in[p]$ and the iterative reweighted maximum likelihood estimator where $W_{j}(t)=\frac{h_{j}^{\prime}(t)}{h_{j}(t)}$ is iteratively updated.

We set the number of nodes to $p\in\{30,60\}$. We consider two popular graphical structures: block and chain. For the block structure, $B$ is a $60$ ($30\times 30$) blockwise diagonal matrix of 12 (6) identical symmetrical Toeplitz matrices whose first row is $(0,0.3,-0.3,0.3,-0.3)$ for $p=60$ ($p=30$). For the chain structure, $B$ is a symmetrical Toeplitz matrix whose first row is $(0,0.3,-0.3,0,0,\cdots,0)$. Their explicit forms are shown in Figure 1 and 2 for $p=30$ and $p=60$, where the excitatory connection are shown in blue with entry $\beta_{j,k}=0.3$ and the inhibitory connection are shown in red with entry $\beta_{j,k}=-0.3$. The background intensity is $\mu_{j}=0.5$. We consider two sets of link function and transfer kernel. For setting 1, we consider an $\arctan$ function as the link function and a restricted linear transfer kernel.

$\displaystyle\begin{cases}h_{j}(x)&=4+\frac{8}{\pi}\arctan(x),\ j\in[p],\\ \kappa_{j,k}(x)&=(1-x)\mathbb{I}_{[0,1]}(x),\ j,k\in[p].\end{cases}$

For setting 2, we consider a sigmoid function as our link function and a exponential transfer kernel.

$\displaystyle\begin{cases}h_{j}(x)&=1+\frac{4e^{x}}{1+e^{x}},\ j\in[p],\\ \kappa_{j,k}(x)&=e^{-x}\mathbb{I}_{[0,5]}(x),\ j,k\in[p].\end{cases}$

The length of observation is $T\in\{200,400\}$ and the number of sub-intervals $M=10T$. We use a 5-fold cross-validation to select the penalty level. We evaluate the estimation accuracy by the relative $\ell_{1}$ error $\frac{\sum_{i,j}|\hat{\beta}_{i,j}-\beta_{i,j}|}{\sum_{i,j}|\beta_{i,j}|}$ and the relative mean square error of the estimated transferring coefficient matrix $\frac{\|\hat{B}-B^{*}\|_{F}^{2}}{\|B^{*}\|_{F}^{2}}$. Three methods are compared: the $\ell_{1}$ penalized naively weighted estimator (Sparse-Naive) the $\ell_{1}$ penalized iteratively reweighted maximum likelihood estimator (Sparse-MLE) and the vanilla maximum likelihood estimator. The average relative $\ell_{1}$ error and $\ell_{2}$ error based on 50 replications with the standard errors in the parenthesis are reported.

Table 1 and 2 summarize the results based on $50$ data replications. It is seen that both of our proposed methods outperform the vanilla maximum likelihood estimator, which is known to be asymptotically efficient (Ogata,1978) in estimating parameters when the transfer function is linear. Overall, the iteratively reweighted penalized maximum likelihood estimator achieves better performance. But the loss surface for the $\ell_{1}$ penalized MLE is not convex, so there is no permission that the working algorithm will find a feasible solution. The naive weighted estimator provides a convex loss function that guarantees that the algorithm returns a solution in the neighborhood of a feasible solution.

To demonstrate the accuracy of variable selection for our proposed methods, we calculated true-positive rates and false-positive rates, i.e, the edges of the network that were correctly identified for each method and tuning parameter. Plotting true positive rate versus false positive rate over a fine grid of values of the tuning parameter produces a ROC curve, with curves near the top left corner indicating a method that performs well. Figure 3 shows the median of the ROC curves for our two proposed methods and table 3 provides the area under the ROC curve (average over the 50 simulation runs). It is seen that both methods are successful in recovering the graphical structure in all the settings we considered.

Finally, we report the computation time. The naively weighted estimator has a convex loss function and does not need to calculate the reweigh parameters and runs faster than the iteratively reweigh MLE. The simulation code is programmed in Python. For the simulation example in setting 1 with $p=60$, $T=400$ and $M=4000$, the average computing time was about 1.5 minutes for Naive estimator and 1.7 minutes for MLE. For the simulation example in setting 2 with $p=60$, $T=400$ and $M=4000$, the average computing time was $2.1$ minutes for Naive estimator and $2.5$ minutes for MLE. All computations were performed with a Intel(R) Core(TM) i9 10900K [email protected].

We consider the task of recovering the connectivity network among the neurons using the spike train data from the optogenetic study of Bolding and Franks(2018) discussed in the Section 1. In this experiment, the spike trains are recorded at 30 kHz on the olfactory bulb (OB) and the piriform cortex (PCx) of the mice, while a laser pulse is applied directly on the OB cells. For specially designed transgenic mice, light pulses will elicited an increase in OB spiking and remain elevated for the duration of the stimulus. We consider the spike train data collected at three intensity levels, $0\ mW/mm^{2}$, $5\ mW/mm^{2}$ and $10\ mW/mm^{2}$ from $15$ OB cells and $45$ PCx cells. For each intensity level, the observation lasts about $T=90$ seconds. We choose the number of sub-intervals $M=2000$ by pre-experiments.

Figure 4 shows the number of firings of OB cells and PCx cells. The firing numbers vary largely among cells. We use an adaptive link-function to better characterize the average intensity of different neurons:

$$h_{i}(t)=\hat{\lambda}_{i}\big{(}1+\frac{2\arctan(t)}{\pi}\big{)},\ i=1,\ldots,60,$$

where $\hat{\lambda}_{i}=\frac{N_{i}([0,T])}{T}$. Based on the research of the duration of influence among the neurons in Bolding and Franks(2018), we choose $\kappa_{j,k}(x)=\mathbb{I}_{[0,0.25]}(x),j,k,\in[p]$. Since our goal was to provide interpretable visualizations and investigate the influence of the laser on the neural connectivity, we computed sparse graphs with approximately $5\%$ connected edges. We performed a bootstrap procedure by randomly selecting $2000$ sub-intervals with replacement from the original data, finding a tuning parameters $\gamma_{n}$ to achieve $5\%$ sparsity level, fitting the point process graphical model and repeating 50 times. The final graph was constructed from the excitatory and inhibitory edges that occurred in at least $50\%$ of the bootstrap replications.

Figure 5 Shows the networks estimated using a point process graphical model for three groups. The bootstrapped graphical model estimated a sparser network with sparsity level $3.28\%$ for the $0\ mW/mm^{2}$ group, $4.08\%$ for the $5\ mW/mm^{2}$ group and $3.67\%$ for the $10\ mW/mm^{2}$ group. We observe a few apparent patterns. First, PCx cells are densely connected in all three groups but the $0mW/mm^{2}$ group has increased connectivity relative to other groups, suggesting that there may be a multi-layer excitatory-inhibitory neuron network in the piriform cortex to stabilize the odor signal. Second, when the laser was applied to OB cells, the number of edges from other cells to the OB cells increase and the sign of if several connections changed, indicating that the laser affected the connectivity structure of OB cells and PCx cells and this effect was heterogenous. Finally, we observed that a few neurons fired less frequently and had little effect on other cells (such as OB-9, PCx-1, PCx-15, PCx-42, PCx-44, etc). Their role in the odor-processing mechanism needs to be further explored.

In this paper, we have proposed the temporal point process graphical model, a new class of graphical models for learning the temporal dependencies among different components of a multivariate point process. We have shown that the locally independent structure of the process is fully encoded in the transfer matrix. We have provided a class of estimators with a non-asymptotic estimation error bound, which extends the classical estimators of temporal point process. Our model naturally includes the case of sparse temporal point process regression, where the predictors and the responses are not the same process. The performance of our estimator has been tested by simulations and a spike train data set.

Several challenges still remain in the analysis of the temporal point process graphical models. Using cross-validation for selecting tuning parameters is proved to be successful in the classical Lasso estimators (Chetverikov et al., 2021). However, a detailed theoretical analysis is needed to demonstrate its performance in our case. The criteria for optimal tuning parameter selection for estimation and variable selection are left for future research. The second challenge is to assess the uncertainty of the estimated parameters, which is less addressed in the existing work.