---
abstract: 'We propose a novel class of network models for temporal dyadic interaction data. Our objective is to capture important features often observed in social interactions: sparsity, degree heterogeneity, community structure and reciprocity. We use mutually-exciting Hawkes processes to model the interactions between each (directed) pair of individuals. The intensity of each process allows interactions to arise as responses to opposite interactions (reciprocity), or due to shared interests between individuals (community structure). For sparsity and degree heterogeneity, we build the non time dependent part of the intensity function on compound random measures following [@Todeschini2016]. We conduct experiments on real-world temporal interaction data and show that the proposed model outperforms competing approaches for link prediction, and leads to interpretable parameters.'
author:
- |
    Xenia  Miscouridou$^{\mathbf{1}}$, François Caron$^{\mathbf{1}}$, Yee Whye Teh$^{\mathbf{1,2}}$\
    $^{1}$Department of Statistics, University of Oxford\
    $^{2}$DeepMind\
    `{miscouri, caron, y.w.teh}@stats.ox.ac.uk`\
bibliography:
- '../biblio.bib'
title: 'Modelling sparsity, heterogeneity, reciprocity and community structure in temporal interaction data'
---

Introduction
============

There is a growing interest in modelling and understanding temporal dyadic interaction data. Temporal interaction data take the form of time-stamped triples $(t,i,j)$ indicating that an interaction occurred between individuals $i$ and $j$ at time $t$. Interactions may be directed or undirected. Examples of such interaction data include commenting a post on an online social network, exchanging an email, or meeting in a coffee shop. An important challenge is to understand the underlying structure that underpins these interactions. To do so, it is important to develop statistical network models with interpretable parameters, that capture the properties which are observed in real social interaction data.

One important aspect to capture is the *community structure* of the interactions. Individuals are often affiliated to some latent communities (e.g. work, sport, etc.), and their affiliations determine their interactions: they are more likely to interact with individuals sharing the same interests than to individuals affiliated with different communities. An other important aspect is *reciprocity*. Many events are responses to recent events of the opposite direction. For example, if Helen sends an email to Mary, then Mary is more likely to send an email to Helen shortly afterwards. A number of papers have proposed statistical models to capture both community structure and reciprocity in temporal interaction data [@Blundell2012; @Dubois2013; @Linderman2014]. They use models based on Hawkes processes for capturing reciprocity and stochastic block-models or latent feature models for capturing community structure.

In addition to the above two properties, it is important to capture the global properties of the interaction data. Interaction data are often *sparse*: only a small fraction of the pairs of nodes actually interact. Additionally, they typically exhibit high degree (number of interactions per node) *heterogeneity*: some individuals have a large number of interactions, whereas most individuals have very few, therefore resulting in empirical degree distributions being heavy-tailed. As shown by @Karrer2011, @Gopalan2013 and @Todeschini2016, failing to account explicitly for degree heterogeneity in the model can have devastating consequences on the estimation of the latent structure.

Recently, two classes of statistical models, based on random measures, have been proposed to capture sparsity and power-law degree distribution in network data. The first one is the class of models based on exchangeable random measures [@Caron2017; @Veitch2015; @Herlau2015; @Borgs2016; @Todeschini2016; @Palla2016; @Janson2017]. The second one is the class of edge-exchangeable models [@Crane2015; @Crane2017; @Cai2016; @Williamson2016; @Janson2017a; @Ng2017]. Both classes of models can handle both sparse and dense networks and, although the two constructions are different, connections have been highlighted between the two approaches [@Cai2016; @Janson2017a].

The objective of this paper is to propose a class of statistical models for temporal dyadic interaction data that can capture all the desired properties mentioned above, which are often found in real world interactions. These are *sparsity*, *degree heterogeneity*, *community structure* and *reciprocity*. Combining all the properties in a single model is non trivial and there is no such construction to our knowledge. The proposed model generalises existing reciprocating relationships models [@Blundell2012] to the sparse and power-law regime. Our model can also be seen as a natural extension of the classes of models based on exchangeable random measures and edge-exchangeable models and it shares properties of both families. The approach is shown to outperform alternative models for link prediction on a variety of temporal network datasets.

The construction is based on Hawkes processes and the (static) model of @Todeschini2016 for sparse and modular graphs with overlapping community structure. In Section \[sec:background\], we present Hawkes processes and compound completely random measures which form the basis of our model’s construction. The statistical model for temporal dyadic data is presented in Section \[sec:model\] and its properties derived in Section \[sec:properties\]. The inference algorithm is described in Section \[sec:inference\]. Section \[sec:experiments\] presents experiments on four real-world temporal interaction datasets.

Background material {#sec:background}
===================

Hawkes processes
----------------

Let $(t_k)_{k\geq 1}$ be a sequence of event times with $t_k\geq 0$, and let $\mathcal H_t=(t_k|t_k\leq t )$ the subset of event times between time $0$ and time $t$. Let $N(t)=\sum_{k\geq 1}1_{t_k\leq t}$ denote the number of events between time $0$ and time $t$, where $1_{A}=1$ if $A$ is true, and 0 otherwise. Assume that $N(t)$ is a counting process with conditional intensity function $\lambda(t)$, that is for any $t\geq 0$ and any infinitesimal interval $dt$ $$\Pr(N(t+dt)-N(t)=1|\mathcal H_t)=\lambda(t)dt.\label{eq:counting}$$ Consider another counting process $\tilde{N}(t)$ with the corresponding $(\tilde{t}_k)_{k\geq 1}, \mathcal{\tilde{H}}_t, \tilde{\lambda}(t)$. Then, $N(t),\tilde{N}(t)$ are mutually-exciting Hawkes processes  [@self_exc_HP] if the conditional intensity functions $\lambda(t)$ and $\tilde{\lambda}(t)$ take the form $$\begin{aligned}
\lambda(t)=\mu + \int_0^t g_\phi(t-u)\, d\tilde{N}(u)\label{eq:hawkesintensity}\qquad
\tilde{\lambda}(t)=\tilde{\mu} + \int_0^t {g}_{\tilde{\phi}}(t-u)\, d{N}(u)\end{aligned}$$ where $\mu=\lambda(0)>0, \tilde{\mu}=\tilde{\lambda}(0)>0$ are the base intensities and $g_\phi,g_{\tilde{\phi}}$ non-negative kernels parameterised by $\phi$ and $\tilde{\phi}$. This defines a pair of processes in which the current rate of events of each process depends on the occurrence of past events of the opposite process.

Assume that $\mu=\tilde{\mu},\, \phi= \tilde{\phi}$ and $g_\phi(t) \geq 0$ for $t > 0$, $g_\phi(t)=0$ for $t<0$. If $g_\phi$ admits a form of fast decay then this results in strong local effects. However, if it prescribes a peak away from the origin then longer term effects are likely to occur. We consider here an exponential kernel $$g_\phi(t-u)= \eta e^{-\delta (t-u)}, t>u
\label{expo_kernel}$$ where $\phi=(\eta,\delta)$. $\eta\geq 0$ determines the sizes of the self-excited jumps and $\delta>0$ is the constant rate of exponential decay. The stationarity condition for the processes is $\eta < \delta$. Figure \[fig:process\_and\_inten\] gives an illustration of two mutually-exciting Hawkes processes with exponential kernel and their conditional intensities.

Compound completely random measures
-----------------------------------

A homogeneous completely random measure (CRM) [@Kingman1967; @Kingman1993] on $\mathbb R_+$ without fixed atoms nor deterministic component takes the form $$\begin{aligned}
W=\sum_{i\geq 1}w_i \delta_{\theta_i}\end{aligned}$$ where $(w_i,\theta_i)_{i\geq 1}$ are the points of a Poisson process on $(0,\infty)\times\mathbb R_+$ with mean measure $\rho(dw)H(d\theta)$ where $\rho$ is a L' evy measure, $H$ is a locally bounded measure and $\delta_x$ is the dirac delta mass at $x$. The homogeneous CRM is completely characterized by $\rho$ and $H$, and we write $W\sim \operatorname{CRM}(\rho,H)$, or simply $W\sim \operatorname{CRM}(\rho)$ when $H$ is taken to be the Lebesgue measure. @Griffin2017 proposed a multivariate generalisation of CRMs, called compound CRM (CCRM). A compound CRM $(W_1,\ldots,W_p)$ with independent scores is defined as $$\begin{aligned}
W_k=\sum_{i\geq 1}{w_{ik}}\delta_{\theta_i},~~k=1,\ldots,p\end{aligned}$$ where $w_{ik}=\beta_{ik}w_{i0}$ and the scores $\beta_{ik}\geq 0$ are independently distributed from some probability distribution $F_k$ and $W_0=\sum_{i\geq 1}w_{i0}\delta_{\theta_i}$ is a CRM with mean measure $\rho_0(dw_0)H_0(d\theta)$. In the rest of this paper, we assume that $F_k$ is a gamma distribution with parameters $(a_k,b_k)$, $H_0(d\theta)=d\theta$ is the Lebesgue measure and $\rho_0$ is the Lévy measure of a generalized gamma process $$\begin{aligned}
\rho_0(dw)=\frac{1}{\Gamma(1-\sigma)}w^{-1-\sigma}e^{-\tau w}dw\label{eq:levyGGP}\end{aligned}$$ where $\sigma\in(-\infty, 1)$ and $\tau>0$.

Statistical model for temporal interaction data {#sec:model}
===============================================

Consider temporal interaction data of the form $\mathcal D=(t_k,i_k,j_k)_{k\geq 1}$ where $(t_k,i_k,j_k)\in \mathbb R_+\times \mathbb N_*^2$ represents a directed interaction at time $t_k$ from node/individual $i_k$ to node/individual $j_k$. For example, the data may correspond to the exchange of messages between students on an online social network.

We use a point process $(t_k,U_k,V_k)_{k\geq 1}$ on $\mathbb R^3_+$, and consider that each node $i$ is assigned some continuous label $\theta_i\geq 0$. the labels are only used for the model construction, similarly to [@Caron2017; @Todeschini2016], and are not observed nor inferred from the data. A point at location $(t_k,U_k,V_k)$ indicates that there is a directed interaction between the nodes $U_k$ and $V_k$ at time $t_k$. See Figure \[fig:model\] for an illustration.

For a pair of nodes $i$ and $j$, with labels $\theta_i$ and $\theta_j$, let $N_{ij}(t)$ be the counting process $$\begin{aligned}
N_{ij}(t)=\sum_{k|(U_k,V_k)=(\theta_i,\theta_j)}1_{t_k\leq t}
\label{eq:PointProcess}\end{aligned}$$ for the number of interactions between $i$ and $j$ in the time interval $[0,t]$. For each pair $i,j$, the counting processes $N_{ij}, N_{ji}$ are mutually-exciting Hawkes processes with conditional intensities $$\begin{aligned}
\lambda_{ij}(t)=\mu_{ij} + \int_0^t g_\phi(t-u)\, dN_{ji}(u),\qquad \lambda_{ji}(t)=\mu_{ji} + \int_0^t g_\phi(t-u)\, dN_{ij}(u)
\label{eq:intensity}\end{aligned}$$ where $g_\phi$ is the exponential kernel defined in Equation . Interactions from individual $i$ to individual $j$ may arise as a response to past interactions from $j$ to $i$ through the kernel $g_\phi$, or via the base intensity $\mu_{ij}$. We also model assortativity so that individuals with similar interests are more likely to interact than individuals with different interests. For this, assume that each node $i$ has a set of positive latent parameters $(w_{i1},\ldots,w_{ip})\in \mathbb R_+^p$, where $w_{ik}$ is the level of its affiliation to each latent community $k=1,\ldots,p$. The number of communities $p$ is assumed known. We model the base rate $$\mu_{ij}=\mu_{ji}=\sum_{k=1}^p w_{ik}w_{jk}.\label{eq:muij}$$Two nodes with high levels of affiliation to the same communities will be more likely to interact than nodes with affiliation to different communities, favouring assortativity.

In order to capture sparsity and power-law properties and as in @Todeschini2016, the set of affiliation parameters $(w_{i1},\ldots,w_{ip})$ and node labels $\theta_i$ is modelled via a compound CRM with gamma scores, that is $W _0 =\sum_{i=1}^{\infty}w_{i0} \delta_{\theta_{i}} \sim\text{CRM} (\rho_{0})$ where the Lévy measure $\rho_0$ is defined by Equation , and for each node $i\geq 1$ and community $k=1,\ldots,p$ $$\begin{aligned}
w_{ik}= w_{i0}  \beta_{ik}\text{, where }\beta_{ik} \overset{\text{ind}}{\sim}\operatorname{Gamma}(a_k,b_k).
\label{weights}\end{aligned}$$ The parameter $w_{i0}\geq 0$ is a degree correction for node $i$ and can be interpreted as measuring the overall popularity/sociability of a given node $i$ irrespective of its level of affiliation to the different communities. An individual $i$ with a high sociability parameter $w_{i0}$ will be more likely to have interactions overall than individuals with low sociability parameters. The scores $\beta_{ik}$ tune the level of affiliation of individual $i$ to the community $k$. The model is defined on $\mathbb R_+^3$. We assume that we observe interactions over a subset $[0,T]\times[0,\alpha]^2\subseteq \mathbb R_+^3$ where $\alpha$ and $T$ tune both the number of nodes and number of interactions. The whole model is illustrated in Figure \[fig:model\].

![Degree distribution for a Hawkes graph with different values of $\sigma$ (left) and $\tau$ (right). The degree of a node $i$ is defined as the number of nodes with whom it has at least one interaction. The value $\sigma$ tunes the slope of the degree distribution, larger values corresponding to a higher slope. The parameter $\tau$ tunes the exponential cut-off in the degree distribution.[]{data-label="fig:degree"}](figures/tau_and_sigma_degreea "fig:"){width=".48\textwidth"}![Degree distribution for a Hawkes graph with different values of $\sigma$ (left) and $\tau$ (right). The degree of a node $i$ is defined as the number of nodes with whom it has at least one interaction. The value $\sigma$ tunes the slope of the degree distribution, larger values corresponding to a higher slope. The parameter $\tau$ tunes the exponential cut-off in the degree distribution.[]{data-label="fig:degree"}](figures/tau_and_sigma_degreeb "fig:"){width=".48\textwidth"}

The model admits the following set of hyperparameters, with the following interpretation:\
$\bullet$ The hyperparameters $\phi=(\eta,\delta)$ where $\eta\geq 0$ and $\delta\geq 0$ of the kernel $g_\phi$ tune the *reciprocity*.\
$\bullet$ The hyperparameters $(a_k,b_k)$ tune the *community structure* of the interactions. $a_k/b_k=\mathbb E[\beta_{ik}]$ tunes the size of community $k$ while $a_k/b_k^2=\text{var}(\beta_{ik})$ tunes the variability of the level of affiliation to this community; larger values imply more separated communities.\
$\bullet$ The hyperparameter $\sigma$ tunes the *sparsity* and the *degree heterogeneity*: larger values imply higher sparsity and heterogeneity. It also tunes the slope of the degree distribution. Parameter $\tau$ tunes the exponential cut-off in the degree distribution. This is illustrated in Figure \[fig:degree\].\
$\bullet$ Finally, the hyperparameters $\alpha$ and $T$ tune the overall number of interactions and nodes.

We follow [@Rasmussen2013] and use vague Exponential$(0.01)$ priors on $\eta$ and $\delta$. Following @Todeschini2016 we set vague Gamma$(0.01, 0.01)$ priors on $\alpha$, $1-\sigma$, $\tau$, $a_k$ and $b_k$. The right limit for time window, $T$ is considered known.

Properties {#sec:properties}
==========

Connection to sparse vertex-exchangeable and edge-exchangeable models {#sec:properties2}
---------------------------------------------------------------------

The model is a natural extension of sparse vertex-exchangeable and edge-exchangeable graph models. Let $z_{ij}(t)=1_{N_{ij}(t)+N_{ji}(t)>0}$ be a binary variable indicating if there is at least one interaction in $[0,t]$ between nodes $i$ and $j$ in either direction. We assume $$\Pr(z_{ij}(t)=1|(w_{ik},w_{jk})_{k=1,\ldots,p})=1-e^{-2t\sum_{k=1}^p w_{ik}w_{jk}}$$ which corresponds to the probability of a connection in the static simple graph model proposed by @Todeschini2016. Additionally, for fixed $\alpha>0$ and $\eta=0$ (no reciprocal relationships), the model corresponds to a rank-$p$ extension of the rank-1 Poissonized version of edge-exchangeable models considered by @Cai2016 and @Janson2017. The sparsity properties of our model follow from the sparsity properties of these two classes of models.

Sparsity
--------

The size of the dataset is tuned by both $\alpha$ and $T$. Given these quantities, both the number of interactions and the number of nodes with at least one interaction are random variables. We now study the behaviour of these quantities, showing that the model exhibits sparsity. Let $I_{\alpha,T}, E_{\alpha,T}, V_{\alpha,T}$ be the overall number of interactions between nodes with label $\theta_i\leq \alpha$ until time $T$, the total number of pairs of nodes with label $\theta_i\leq \alpha$ who had at least one interaction before time $T$, and the number of nodes with label $\theta_i\leq \alpha$ who had at least one interaction before time $T$ respectively.

$$\begin{aligned}
\begin{split}
V_{\alpha,T}&=\sum_{i}1_{\sum_{j\neq i} N_{ij}(T)1_{\theta_j\leq \alpha} >0}1_{\theta_i\leq \alpha}
\qquad I_{\alpha,T}=\sum_{i\neq j}N_{ij}(T)1_{\theta_i\leq \alpha}1_{\theta_j\leq \alpha} \\
E_{\alpha,T}&=\sum_{i<j}1_{N_{ij}(T)+N_{ji}(T)>0}1_{\theta_i\leq \alpha}1_{\theta_j\leq \alpha}
\end{split}\end{aligned}$$

We provide in the supplementary material a theorem for the exact expectations of $I_{\alpha,T}$, $E_{\alpha,T}$ $V_{\alpha,T}$. Now consider the asymptotic behaviour of the expectations of $ V_{\alpha,T}$, $E_{\alpha,T}$ and $I_{\alpha,T}$ , as $\alpha$ and $T$ go to infinity.[^1] Consider fixed $T>0$ and $\alpha$ that tends to infinity. Then, $$\begin{aligned}
\mathbb E[ V_{\alpha,T}]= \left\{
                                 \begin{array}{ll}
                                  \Theta( \alpha) & \text{if }\sigma<0 \\
                                   \omega(\alpha) & \text{if }\sigma\geq 0
                                 \end{array}
                               \right.,\qquad
\mathbb E[ E_{\alpha,T}]=\Theta(\alpha^2),\qquad
\mathbb E[ I_{\alpha,T}]=\Theta(\alpha^2)\end{aligned}$$ as $\alpha$ tends to infinity. For $\sigma<0$, the number of edges and interactions grows quadratically with the number of nodes, and we are in the dense regime. When $\sigma\geq 0$, the number of edges and interaction grows subquadratically, and we are in the sparse regime. Higher values of $\sigma$ lead to higher sparsity. For fixed $\alpha$, $$\begin{aligned}
\mathbb E[ V_{\alpha,T}]= \left\{
                                 \begin{array}{ll}
                                  \Theta( 1 ) & \text{if }\sigma<0 \\
                                   \Theta(\log T) & \text{if }\sigma=0 \\
                                  \Theta( T^{\sigma}) & \text{if }\sigma>0 \\
                                 \end{array}
                               \right.,\qquad
\mathbb E[ E_{\alpha,T}]= \left\{
                                 \begin{array}{ll}
                                  \Theta( 1 ) & \text{if }\sigma<0 \\
                                   O(\log T) & \text{if }\sigma=0 \\
                                  O( T^{3\sigma/2} ) & \text{if }\sigma\in(0,1/2] \\
                                  O( T^{(1+\sigma)/2}) & \text{if }\sigma\in(1/2,1) \\
                                 \end{array}
                               \right.\end{aligned}$$ and $\mathbb E[ I_{\alpha,T}]=\Theta(T )$ as $T$ tends to infinity. Sparsity in $T$ arises when $\sigma\geq 0$ for the number of edges and when $\sigma> 1/2$ for the number of interactions. The derivation of the asymptotic behaviour of expectations of $ V_{\alpha,T}$, $E_{\alpha,T}$ and $I_{\alpha,T}$ follows the lines of the proofs of Theorems 3 and 5.3 in [@Todeschini2016] $(\alpha\rightarrow\infty)$ and Lemma D.6 in the supplementary material of [@Cai2016] $(T\rightarrow\infty)$, and is omitted here.

Approximate Posterior Inference {#sec:inference}
===============================

Assume a set of observed interactions $\mathcal D=(t_k,i_k,j_k)_{k\geq 1}$ between $V$ individuals over a period of time $T$. The objective is to approximate the posterior distribution $\pi(\phi,\xi|\mathcal D)$ where $\phi$ are the kernel parameters and $\xi=((w_{ik})_{i=1,\ldots,V,k=1,\ldots,p},(a_k,b_k)_{k=1,\ldots,p},\alpha,\sigma,\tau)$, the parameters and hyperparameters of the compound CRM. One possible approach is to follow a similar approach to that taken in [@Rasmussen2013]; derive a Gibbs sampler using a data augmentation scheme which associates a latent variable to each interaction. However, such an algorithm is unlikely to scale well with the number of interactions. Additionally, we can make use of existing code for posterior inference with Hawkes processes and graphs based on compound CRMs, and therefore propose a two-step approximate inference procedure, motivated by modular Bayesian inference [@Jacob2017].

Let $\mathcal Z=(z_{ij}(T))_{1\leq i,j\leq V}$ be the adjacency matrix defined by $z_{ij}(T)=1$ if there is at least one interaction between $i$ and $j$ in the interval $[0,T]$, and 0 otherwise. We have $$\begin{aligned}
\pi(\phi,\xi|\mathcal D)&=\pi(\phi,\xi|\mathcal D,\mathcal Z)=\pi(\xi|\mathcal D,\mathcal Z)\pi(\phi|\xi ,\mathcal D).\end{aligned}$$ The idea of the two-step procedure is to (i) Approximate $\pi(\xi|\mathcal D,\mathcal Z)$ by $\pi(\xi|\mathcal Z)$ and obtain a Bayesian point estimate $\widehat\xi$ then (ii) Approximate $\pi(\phi|\xi ,\mathcal D)$ by $\pi(\phi|\widehat \xi ,\mathcal D)$.

The full posterior is thus approximated by $\widetilde\pi(\phi,\xi)=\pi(\xi|\mathcal Z)\pi(\phi|\widehat \xi ,\mathcal D).$ As mentioned in Section \[sec:properties2\], the statistical model for the binary adjacency matrix $\mathcal Z$ is the same as in [@Todeschini2016]. We use the MCMC scheme of [@Todeschini2016] and the accompanying software package SNetOC[^2] to perform inference. The MCMC sampler is a Gibbs sampler which uses a Metropolis-Hastings (MH) step to update the hyperparameters and a Hamiltonian Monte Carlo (HMC) step for the parameters. From the posterior samples we compute a point estimate $(\widehat w_{i1},\ldots,\widehat w_{ip})$ of the weight vector for each node. We follow the approach of @Todeschini2016 and compute a minimum Bayes risk point estimate using a permutation-invariant cost function. Given these point estimates we obtain estimates of the base intensities $\widehat \mu_{ij}$. Posterior inference on the parameters $\phi$ of the Hawkes kernel is then performed using Metropolis-Hastings, as in [@Rasmussen2013]. Details of the two-stage inference procedure are given in the supplementary material.

**Empirical investigation of posterior consistency.** To validate the two-step approximation to the posterior distribution, we study empirically the convergence of our approximate inference scheme using synthetic data. Experiments suggest that the posterior concentrates around the true parameter value. More details are given in the supplementary material.

Experiments {#sec:experiments}
===========

We perform experiments on four temporal interaction datasets from the Stanford Large Network Dataset Collection[^3] [@SNAP2014]:\
$\bullet$ The [Email]{} dataset consists of emails sent within a large European research institution over 803 days. There are 986 nodes, 24929 edges and 332334 interactions. A separate interaction is created for every recipient of an email.\
$\bullet$ The [College]{} dataset consists of private messages sent over a period of 193 days on an online social network (Facebook-like platform) at the University of California, Irvine. There are 1899 nodes, 20296 edges and 59835 interactions. An interaction $(t,i,j)$ corresponds to a user $i$ sending a private message to another user $j$ at time $t$.\
$\bullet$ The [Math]{} overflow dataset is a temporal network of interactions on the stack exchange website Math Overflow over 2350 days. There are 24818 nodes, 239978 edges and 506550 interactions.An interaction $(t,i,j)$ means that a user $i$ answered another user’s $j$ question at time $t$, or commented on another user’s $j$ question/response.\
$\bullet$ The [Ubuntu]{} dataset is a temporal network of interactions on the stack exchange website Ask Ubuntu over 2613 days. There are 159316 nodes, 596933 edges and 964437 interactions. An interaction $(t,i,j)$ means that a user $i$ answered another user’s $j$ question at time $t$, or commented on another user’s $j$ question/response. **Comparison on link prediction.** We compare our model (Hawkes-CCRM) to five other benchmark methods: (i) our model, without the Hawkes component (obtained by setting $\eta=0$), (ii) the Hawkes-IRM model of @Blundell2012 which uses an infinite relational model (IRM) to capture the community structure together with a Hawkes process to capture reciprocity, (iii) the same model, called Poisson-IRM, without the Hawkes component, (iv) a simple Hawkes model where the conditional intensity given by Equation is assumed to be same for each pair of individuals, with unknown parameters $\delta$ and $\eta$, (v) a simple Poisson process model, which assumes that interactions between two individuals arise at an unknown constant rate. Each of these competing models capture a subset of the features we aim to capture in the data: sparsity/heterogeneity, community structure and reciprocity, as summarized in Table \[tab:results\_and\_prop\]. The models are given in the supplementary material. The only model to account for all the features is the proposed Hawkes-CCRM model.

  ------------- ------------- ------------- ------------- -------------
                                                          
                    email        college        math         ubuntu
  Hawkes-CCRM     **10.95**     **1.88**      **20.07**     **29.1**
  CCRM             $12.08$       $ 2.90$       $89.0$        $36.5$
  Hawkes-IRM       $14.2$        $ 3.56$        96.9          59.5
  Poisson-IRM      $31.7$        $15.7$         204.7         79.3
  Hawkes           $154.8$      $153.29$       220.10        191.39
  Poisson        $\sim10^ 3$   $ \sim10^4$   $ \sim10^4$   $\sim10^4 $
  ------------- ------------- ------------- ------------- -------------

  : (Left) Performance in link prediction. (Right) Properties captured by the different models.

\[tab:results\_and\_prop\]

[c@c@c@]{} sparsity/ & community & reciprocity\
heterogeneity & structure\
& &\
& &\
&&\
&&\
& &\
\

We perform posterior inference using a Markov chain Monte Carlo algorithm. For our Hawkes-CCRM model, we follow the two-step procedure described in Section \[sec:inference\]. For each dataset, there is some background information in order to guide the choice of the number $p$ of communities. The number of communities $p$ is set to $p=4$ for the [Email]{} dataset, as there are 4 departments at the institution, $p=2$ for the [College]{} dataset corresponding to the two genders, and $p=3$ for the [Math]{} and [Ubuntu]{} datasets, corresponding to the three different types of possible interactions. We follow @Todeschini2016 regarding the choice of the MCMC tuning parameters and initialise the MCMC algorithm with the estimates obtained by running a MCMC algorithm with $p = 1$ feature with fewer iterations. For all experiments we run 2 chains in parallel for each stage of the inference. We use 100000 iterations for the first stage and 10000 for the second one. For the Hawkes-IRM model, we also use a similar two-step procedure, which first obtains a point estimate of the parameters and hyperparameters of the IRM, then estimates the parameters of the Hawkes process given this point estimate. This allows us to scale this approach to the large datasets considered. We use the same number of MCMC samples as for our model for each step.

We compare the different algorithms on link prediction. For each dataset, we make a train-test split in time so that the training datasets contains $85 \%$ of the total temporal interactions. We use the training data for parameter learning and then use the estimated parameters to perform link prediction on the held out test data. We report the root mean square error between the predicted and true number of interactions for each directed pair in the test set . The results are reported in Table \[tab:results\_and\_prop\]. On all the datasets, the proposed Hawkes-CCRM approach outperforms other methods. Interestingly, the addition of the Hawkes component brings improvement for both the IRM-based model and the CCRM-based model.

**Community structure, degree distribution and sparsity.** Our model also estimates the latent structure in the data through the weights $w_{ik}$, representing the level of affiliation of a node $i$ to a community $k$. For each dataset, we order the nodes by their highest estimated feature weight, obtaining a clustering of the nodes. We represent the ordered matrix $\left(z_{ij}\left(T\right)\right)$ of binary interactions in Figure \[fig:adjacency\] (a)-(d). This shows that the method can uncover the latent community structure in the different datasets. Within each community, nodes still exhibit degree heterogeneity as shown in Figure \[fig:adjacency\] (e)-(h). where the nodes are sorted within each block according to their estimated sociability $\widehat w_{i0}$. The ability of the approach to uncover latent structure was illustrated by @Todeschini2016, who demonstrate that models which do not account for degree heterogeneity, cannot capture latent community estimation but they rather cluster the nodes based on their degree. We also look at the posterior predictive degree distribution based on the estimated hyperparameters, and compare it to the empirical degree distribution in the test set. The results are reported in Figure \[fig:adjacency\] (i)-(l) showing a reasonable fit to the degree distribution. Finally we report the $95\%$ posterior credible intervals (PCI) for the sparsity parameter $\sigma$ for all datasets. Each PCI is $(-0.69, -0.50), (-0.35, -0.20 ),(0.15, 0.18), (0.51,  0.57)$ respectively. The range of $\sigma$ is $ (-\infty,1)$. [Email]{} and [college]{} give negative values corresponding to denser networks whereas [Math]{} and [Ubuntu]{} datasets are sparser.

Conclusion
==========

We have presented a novel statistical model for temporal interaction data which captures multiple important features observed in such datasets, and shown that our approach outperforms competing models in link prediction. The model could be extended in several directions. One could consider asymmetry in the base intensities $\mu_{ij}\neq \mu_{ji}$ and/or a bilinear form as in [@Zhou2015]. Another important extension would be the estimation of the number of latent commmunities/features $p$.

**Acknowledgments.** The authors thank the reviewers and area chair for their constructive comments. XM, FC and YWT acknowledge funding from the ERC under the European Union’s 7th Framework programme (FP7/2007-2013) ERC grant agreement no. 617071. FC acknowledges support from EPSRC under grant EP/P026753/1 and from the Alan Turing Institute under EPSRC grant EP/N510129/1. XM acknowledges support from the A. G. Leventis Foundation.

Supplementary Material {#supplementary-material .unnumbered}
======================

A: Background on compound completely random measures {#a-background-on-compound-completely-random-measures .unnumbered}
====================================================

We give the necessary background on compound completely random measures (CCRM). An extensive account of this class of models is given in [@Griffin2017]. In this article, we consider a CCRM $(W_1,\ldots,W_p)$ on $\mathbb R_+$ characterized, for any $(t_1,\ldots,t_p)\in\mathbb R_+^p$ and measurable set $A\subset \mathbb R_+$, by $$\mathbb E[e^{-\sum_{k=1}^p t_k W_k(A)}] =\exp (-H_0(A)\psi(t_1,\ldots,t_p))$$ where $H_0$ is the Lebesgue measure and $\psi$ is the multivariate Laplace exponent defined by $$\psi(t_1,\ldots,t_p)  = \int_{\mathbb{R}^p_+}(1-e^{-\sum_{k=1}^p t_k w_k})\rho(dw_{1},\ldots,dw_p).
\label{eq:laplaceexponent}$$ The multivariate Lévy measure $\rho$ takes the form $$\rho(dw_{1},\ldots,dw_{p})=\int_{0}^{\infty
}w_{0}^{-p}\prod_{k=1}^p F_k\left(  \frac{dw_{k}}{w_{0}}\right)
\rho_{0}(dw_{0})
\label{eq:levymeasure}$$ where $F_k$ is the distribution of a Gamma random variable with parameters $a_k$ and $b_k$ and $\rho_0$ is the Lévy measure on $(0,\infty)$ of a generalized gamma process $$\rho_0(dw_0)=\frac{1}{\Gamma(1-\sigma)}w_0^{-1-\sigma}\exp(-w_0\tau)dw_0$$ where $\sigma\in(-\infty, 1)$ and $\tau>0$.

Denote ${\boldsymbol{w}}_i=(w_{i1},\dots,w_{ip})^T$, ${\boldsymbol{\beta}}_i=\left(\beta_{i1},\dots,\beta_{ip} \right)^T$ and $\rho(d{\boldsymbol{w}})=\rho(dw_1,\dots,dw_p)$. $w_0$ always refers to the scalar weight corresponding to the measure $\rho_0$.

B: Expected number of interactions, edges and nodes {#b-expected-number-of-interactions-edges-and-nodes .unnumbered}
===================================================

Recall that $I_{\alpha,T}, E_{\alpha,T}$ and $V_{\alpha,T}$ are respectively the overall number of interactions between nodes with label $\theta_i\leq \alpha$ until time $T$, the total number of pairs of nodes with label $\theta_i\leq \alpha$ who had at least one interaction before time $T$, and the number of nodes with label $\theta_i\leq \alpha$ who had at least one interaction before time $T$ respectively, and are defined as $$\begin{aligned}
I_{\alpha,T}&=\sum_{i\neq j}N_{ij}(T)1_{\theta_i\leq \alpha}1_{\theta_j\leq \alpha} \\
E_{\alpha,T}&=\sum_{i<j}1_{N_{ij}(T)+N_{ji}(T)>0}1_{\theta_i\leq \alpha}1_{\theta_j\leq \alpha}\\
V_{\alpha,T}&=\sum_{i}1_{\sum_{j\neq i} N_{ij}(T)1_{\theta_j\leq \alpha} >0}1_{\theta_i\leq \alpha}\end{aligned}$$

\[expected\_num\] The expected number of interactions $I_{\alpha,T}$, edges $E_{\alpha,T}$ and nodes $V_{\alpha,T}$ are given as follows: $$\begin{aligned}
\mathbb{E}[I_{\alpha,T}]{} & = \alpha^2 {\boldsymbol{\mu}}_w^T{\boldsymbol{\mu}}_w \left ( \frac{\delta}{\delta - \eta }T  -
      \frac{\eta}{(\eta -\delta)^2}  \left( 1 - e^{-T(\delta - \eta)}\right)\right )\\
\mathbb{E}[E_{\alpha,T}]&=
   \frac{\alpha^2}{2} \int_{\mathbb{R}^p_+} \psi(2Tw_1,\ldots,2Tw_p)  \rho(dw_{1},\ldots,dw_p)\\
\mathbb{E}[V_{\alpha,T}]&=
\alpha  \int_{\mathbb{R}^p_+} \left (1-e^{- \alpha \psi(2Tw_1,\ldots,2Tw_p)}\right )\rho(dw_1,\ldots,dw_p)\end{aligned}$$ where $\mu_w =\int_{\mathbb{R}^p_+} {\boldsymbol{w}}\rho(dw_1,\ldots,dw_p)$.

The proof of Theorem \[expected\_num\] is given below and follows the lines of Theorem 3 in [@Todeschini2016].

Mean number of nodes $\mathbb{E}[V_{\alpha,T}]$ {#mean-number-of-nodes-mathbbev_alphat .unnumbered}
-----------------------------------------------

We have $$\begin{split}
\mathbb{E} \left[ V_{\alpha,T} \right] &=  \mathbb{E} \left[ \sum_i (1-1_{N_{ij}(T)=0, \forall j\neq i|\theta_j\leq \alpha)} \right]1_{\theta_i\leq \alpha} \nonumber\\
&= \sum_i \left\{ 1- \mathbb{P}(N_{ij}(T)=0, \forall j\neq i|\theta_j\leq \alpha) \right\}1_{\theta_i\leq \alpha}
\end{split}$$ Using the Palm/Slivnyak-Mecke formula and Campbell’s theorem, see e.g. [@Moller2003 Theorem 3.2] and [@Kingman1993], we obtain $$ \begin{split}
\mathbb{E}[V_{\alpha,T}]&=  \mathbb{E} \left[ \mathbb{E} \left[ V_{\alpha,T} |W_1,\ldots,W_p \right] \right] \\
&=  \mathbb{E} \left[\sum_i \left( 1-   e^{-2T{\boldsymbol{w}}_i^T \sum_{j\neq i} {\boldsymbol{w}}_j 1_{\theta_j \leq \alpha}  } \right) 1_{\theta_i \leq \alpha} \right]  \\
&= \alpha \int  \mathbb{E} \left( 1-   e^{-2T{\boldsymbol{w}}^T \sum_j {\boldsymbol{w}}_j 1_{\theta_j \leq \alpha}  } \right) \rho(d{\boldsymbol{w}}) \\
&= \alpha
\int_{\mathbb{R}_+^p} \left(1-e^{- \alpha \psi(2Tw_1,\ldots,2Tw_p)}\right) \rho \left(dw_1,\ldots,dw_p \right)
\end{split}
\label{nodes}$$

Mean number of edges $\mathbb{E}[E_{\alpha,T}]$ {#mean-number-of-edges-mathbbee_alphat .unnumbered}
-----------------------------------------------

Using the extended Slivnyak-Mecke formula, see e.g. [@Moller2003 Theorem 3.3] $$\begin{split}
 \mathbb{E}\left[E_{\alpha,T}\right]   &=\mathbb{E}\left[ \mathbb{E}\left[E_{\alpha,T}| W_1,\ldots,W_p \right] \right]  \nonumber\\
 & = \mathbb{E} \left[  \sum_i 1_{\theta_i<\alpha}\frac{1}{2} \sum_{j\neq i} 1_{\theta_j \leq \alpha}(1-e^{-2T{\boldsymbol{w}}_i^T{\boldsymbol{w}}_j})\right] \nonumber\\
   &= \frac{\alpha^2 }{2}\int_{\mathbb{R}^p_+} \psi(2Tw_1,\dots,2Tw_p)  \rho(dw_1,\ldots,dw_p)
\end{split}$$

Mean number of interactions $ \mathbb{E}[I_{\alpha,T}] $  {#mean-number-of-interactions-mathbbei_alphat .unnumbered}
---------------------------------------------------------

We have $$\begin{aligned}
 \mathbb{E}[I_{\alpha,T}]
& =\mathbb{E}\left[\mathbb{E}\left[I_{\alpha,T} | W_1,\ldots,W_p\right]\right] \nonumber\\
& = \mathbb{E}\left[ \sum_{i\neq j}\left(\mathbb{E}\left[\int_0^T \lambda_{ij}(t)dt\mid  W_1,\ldots,W_p \right] \right) \right]\nonumber\\
& = \mathbb{E}\left[ \sum_{i\neq j}\left( \int_0^T \mathbb{E}\left[   \mu_{ij} \frac{\delta}{\delta - \eta}  - \mu_{ij}\frac{\eta}{\delta - \eta}e^{-{t(\delta - \eta)}}    \mid W_1,\ldots,W_p \right]  dt \right) \right]\nonumber\\
& = \mathbb{E}\left[ \sum_{i\neq j}{\boldsymbol{w}}_i^T {\boldsymbol{w}}_j  \right]    \int_0^T \left[    \frac{\delta}{\delta - \eta}  - \frac{\eta}{\delta - \eta}e^{-{t(\delta - \eta)}} \right]  dt \nonumber\\
& =  \alpha^2 {\boldsymbol{\mu}}_w^T{\boldsymbol{\mu}}_w \left (\frac{\delta}{\delta - \eta} T - \frac{\eta}{(\delta - \eta)^2}(1 - e^{-T(\delta - \eta)})\right )
\label{dir_edges}\end{aligned}$$ where the third line follows from [@Dassios_Zhao], the last line follows from another application of the extended Slivnyak-Mecke formula for Poisson point processes and ${\boldsymbol{\mu}}_w =\int_{\mathbb{R}^p_+} {\boldsymbol{w}}\rho(dw_1,\ldots,dw_p)$.

C: Details of the approximate inference algorithm {#c-details-of-the-approximate-inference-algorithm .unnumbered}
=================================================

Here we provide additional details on the two-stage procedure for approximate posterior inference. The code is publicly available at <https://github.com/OxCSML-BayesNP/HawkesNetOC>.

Given a set of observed interactions $\mathcal D=(t_k,i_k,j_k)_{k\geq 1}$ between $V$ individuals over a period of time $T$, the objective is to approximate the posterior distribution $\pi(\phi,\xi|\mathcal D)$ where $\phi=(\eta,\delta)$ are the kernel parameters and $\xi=((w_{ik})_{i=1,\ldots,V,k=1,\ldots,p},(a_k,b_k)_{k=1,\ldots,p},\alpha,\sigma,\tau)$, the parameters and hyperparameters of the compound CRM. Given data $\mathcal D$, let $\mathcal Z=(z_{ij}(T))_{1\leq i,j\leq V}$ be the adjacency matrix defined by $z_{ij}(T)=1$ if there is at least one interaction between $i$ and $j$ in the interval $[0,T]$, and 0 otherwise.

For posterior inference, we employ an approximate procedure, which is formulated in two steps and is motivated by modular Bayesian inference (Jacob et al., 2017). It also gives another way to see the two natures of this type of temporal network data. Firstly we focus on the static graph i.e. the adjacency matrix of the pairs of interactions. Secondly, given the node pairs that have at least one interaction, we learn the rate for the appearance of those interactions assuming they appear in a reciprocal manner by mutual excitation.

We have $$\begin{aligned}
\pi(\phi,\xi|\mathcal D)&=\pi(\phi,\xi|\mathcal D,\mathcal Z)=\pi(\xi|\mathcal D,\mathcal Z)\pi(\phi|\xi ,\mathcal D).\end{aligned}$$ The idea of the two-step procedure is to

1.  Approximate $\pi(\xi|\mathcal D,\mathcal Z)$ by $\pi(\xi|\mathcal Z)$ and obtain a Bayesian point estimate $\widehat\xi$.

2.  Approximate $\pi(\phi|\xi ,\mathcal D)$ by $\pi(\phi|\widehat \xi ,\mathcal D)$.

C1: Stage 1 {#c1-stage-1 .unnumbered}
-----------

As mentioned in Section \[sec:inference\] in the main article, the joint model $\pi(\mathcal Z, \xi)$ on the binary undirected graph is equivalent to the model introduced by [@Todeschini2016], and we will use their Markov chain Monte Carlo (MCMC) algorithm and the publicly available code SNetOC[^4] in order to approximate the posterior $\pi(\xi| \mathcal Z)$ and obtain a Bayesian point estimate $\widehat \xi$. Let $w_{*k}= W_k([0,\alpha])-\sum_i w_{ik}$ corresponding to the overall level of affiliation to community $k$ of all the nodes with no interaction (recall that in our model, the number of nodes with no interaction may be infinite). For each undirected pair $i,j$ such that $z_{ij}(T)=1$, consider latent count variables $\widetilde n_{ijk}$ distributed from a truncated multivariate Poisson distribution, see [@Todeschini2016 Equation (31)]. The MCMC sampler to produce samples asymptotically distributed according to $\pi(\xi| \mathcal Z)$ then alternates between the following steps:

1.  Update $(w_{ik})_{i=1,\ldots,V,k=1,\ldots,p}$ using an Hamiltonian Monte Carlo (HMC) update,

2.  Update $(w_{*k},a_k,b_k)_{k=1,\ldots,p},\alpha,\sigma,\tau$ using a Metropolis-Hastings step

3.  Update the latent count variables using a truncated multivariate Poisson distribution

We use the same parameter settings as in [@Todeschini2016]. We use $\epsilon=10^{-3}$ as truncation level to simulate the $w_{*k}$, and set the number of leapfrog steps $L=10$ in the HMC step. The stepsizes of both the HMC and the random walk MH are adapted during the first 50 000 iterations so as to target acceptance ratios of 0.65 and 0.23 respectively. The minimum Bayes point estimates $(\widehat w_{ik})_{i=1,\ldots,V,k=1,\ldots,p}$ are then computed using a permutation-invariant cost function, as described in [@Todeschini2016 Section 5.2]. This allows to compute point estimates of the base intensity measures of the Hawkes processes, for each $i\neq j$ $$\widehat \mu_{ij} =\sum_{k=1}^p \widehat w_{ik}\widehat w_{jk}.$$

C2: Stage 2 {#c2-stage-2 .unnumbered}
-----------

In stage 2, we use a MCMC algorithm to obtain samples approximately distributed according to $$\pi(\phi|\widehat \xi,\mathcal D)=\pi(\phi|(\widehat \mu_{ij}),\mathcal D)$$ where $\phi$ are the parameters of the Hawkes kernel. For each ordered pair $(i,j)$ such that $n_{ij}=N_{ij}(T)>0$, let $(t_{ij}^{(1)}<t_{ij}^{(2)}<\ldots<t_{ij}^{(n_{ij})})$ be the times of the observed directed interactions from $i$ to $j$. The intensity function is $$\lambda_{ij}(t) = \widehat \mu_{ij} + \sum_{\ell|t_{ji}^{(\ell)}<t} \frac{\eta}{\delta} f_\delta(t-t_{ji}^{(\ell))})
$$ where $$\frac{\eta}{\delta}  f_\delta(t-t_{ji}^{(\ell)}) =\eta \times  e^{-\delta (t-t_{ji}^{(\ell)})}.$$ We write the kernel in this way to point out that it is equal to a density function $f_\delta$, (here exponential density), scaled by the step size $\eta/\delta$. We denote the distribution function of $f_\delta$ by $F_\delta$. Then, by Proposition 7.2.III in [@Daley2008] we obtain the likelihood $$L(\mathcal D\mid \phi, (\widehat \mu_{ij}))=\prod_{(i,j)|N_{ij}(T)>0} \left [\exp(-\Lambda_{ij}(T)) \prod_{l =1}^{n_{ij}} \lambda_{ij} (t_{ij}^{(\ell)}) \right ]$$

where $$\Lambda_{ij}(t) = \int_0^t  \lambda_{ij}(u)\, du = t \widehat \mu_{ij} + \sum_{\ell | t_{ji}^{(\ell)}<t} \frac{\eta}{\delta} F_\delta(t- t_{ji}^{(\ell)} ).
$$

We derive a Gibbs sampler with Metropolis Hastings steps to estimate the parameters $\eta,\delta$ conditionally on the estimates $\widehat{\mu}_{ij}$. As mentioned in Section \[sec:model\] of the main article we follow [@Rasmussen2013] for the choice of vague exponential priors $p(\eta), p(\delta)$. For proposals we use truncated Normals with variances $\left(\sigma^2_\eta,\sigma^2_\delta \right)= (1.5,2.5) $. The posterior is given by $$\begin{aligned}
\pi(\phi \mid\mathcal D, (\widehat \mu_{ij})) \propto \exp\left(- \sum_{(i,j)|N_{ij}(T)>0}\Lambda_{ij}(T)\right)  \left [\prod_{(i,j)|N_{ij}(T)>0} \prod_{\ell=1}^{n_{ij}} \lambda_{ij}(t_{ij}^{(\ell)})  \right ] \times p(\eta)p(\delta) \\\end{aligned}$$

We use an efficient way to compute the intensity at each time point $t_{ij}^{(\ell)}$ by writing it in the form $$\lambda_{ij}(t_{ij}^{(\ell)})= \widehat \mu_{ij} + \eta S_{ij}^{(\ell)}(\delta),$$ where $$S_{ij}^{(\ell)}(\delta) =e^{-\delta t_{ij}^{(\ell)} }\sum_{k=1}^{n_{ji}} e^{\delta t_{ji}^{(k)}}  1_{t_{ji}^{(k)} < t_{ij}^{(\ell)}},$$ and then derive a recursive relationship of $S_{ij}^{(\ell)}(\delta)$ in terms of $S_{ij}^{(\ell-1)}(\delta) $. In this way, we can precompute several terms by ordering the event times and arrange them in bins defined by the event times of the opposite process.

D: Posterior consistency {#d-posterior-consistency .unnumbered}
========================

We simulate interaction data from the Hawkes-CCRM model described in Section \[sec:model\] in the main article, using parameters $p=4, \alpha=50, \sigma=0.3, \tau=1, a_k = 0.08, \phi=\left( 0.85, 3\right), T=300$. We perform the two-step inference procedure with data of increasing sample size, and check empirically that the approximate posterior $\pi(\phi|\widehat \xi, \mathcal D)$ concentrates around the true value as the sample size increases. Figure \[fig:post\_conv\] below shows the plots of the approximate marginal posterior distribution of $\delta$ and $\eta$. Experiments suggest that the posterior still concentrates around the true parameter value under this approximate inference scheme.

[![(Left) Approximate marginal posterior distribution of $\delta$ given $n$ interaction data. (Right) Approximate marginal posterior distribution of $\eta$ given $n$ interaction data. Posterior concentrates around the true value, with increasing sample size $n$.[]{data-label="fig:post_conv"}](figures/post_distr_delta "fig:"){width=".48\textwidth"}]{} [![(Left) Approximate marginal posterior distribution of $\delta$ given $n$ interaction data. (Right) Approximate marginal posterior distribution of $\eta$ given $n$ interaction data. Posterior concentrates around the true value, with increasing sample size $n$.[]{data-label="fig:post_conv"}](figures/post_distr_eta "fig:"){width=".51\textwidth"}]{}

E: Experiments {#e-experiments .unnumbered}
==============

We perform experiments in which we compare our Hawkes-CCRM model to five other competing models. The key part in all cases is the conditional intensity of the point process. which we give below. In all cases we use $\{t^{(k)}_{ij} \}_{k\geq 1}$ to refer to the set of events from $i$ to $j$, i.e. the interactions for the directed pair $(i,j)$.

Hawkes-CCRM {#hawkes-ccrm .unnumbered}
-----------

For each directed pair of nodes $(i,j), \,i \neq j$\
$N_{ij}(t) \sim \text{Hawkes Process}(\lambda_{ij}(t))$ where $ \lambda_{ij}(t)=\sum_k w_{ik} w_{jk} + \sum_{t_{ji}^{(k)}<t} \eta e^{-\delta(t-t_{ji}^{(k)})}$.

Hawkes-IRM [@Blundell2012] {#hawkes-irm .unnumbered}
--------------------------

For each directed pair of clusters $(p,q), \,p \neq q$\
$N_{pq}(t) \sim \text{Hawkes Process}(\lambda_{pq}(t))$ where $ \lambda_{pq}(t)=n_p n_q\gamma_{pq} + \sum_{t_{qp}^{(k)}<t}\eta  e^{-\delta(t-t_{qp}^{(k)})}$.\
For the details of the model see [@Blundell2012].

Poisson-IRM (as explained in [@Blundell2012]) {#poisson-irm-as-explained-in .unnumbered}
---------------------------------------------

For each directed pair of clusters $(p,q), \,p \neq q$\
$N_{pq}(t) \sim \text{Poisson}(\lambda_{pq}(t))$ where $ \lambda_{pq}(t)=n_p n_q\gamma_{pq}$.\
For the details of the model see [@Blundell2012].

CCRM [@Todeschini2016] {#ccrm .unnumbered}
----------------------

$N_{ij}(t) \sim \text{Poisson Process}(\lambda_{ij}(t))$ where $ \lambda_{ij}(t)=\sum_k w_{ik} w_{jk}$.

Hawkes {#hawkes .unnumbered}
------

$N_{ij}(t) \sim \text{Hawkes Process}\left(\lambda_{ij}(t)\right)$ where $ \lambda_{ij} (t)= \mu + \sum_{t_{ji}^{(k)}<t} \eta e^{-\delta(t-t_{ji}^{(k)})}$.

Poisson {#poisson .unnumbered}
-------

$N_{ij}(t) \sim \text{Poisson Process}(\lambda_{ij}(t))$ where $ \lambda_{ij} (t)= \mu$.

[^1]: We use the following asymptotic notations. $X_\alpha = O(Y_\alpha)$ if $\lim X_\alpha/Y_\alpha  < \infty$, $ X_\alpha = \omega(Y_\alpha)$ if $\lim Y_\alpha/X_\alpha =0$ and $X_\alpha = \Theta(Y_\alpha)$ if both $X_\alpha= O( Y_\alpha )$ and $Y_\alpha =O( X_\alpha )$.

[^2]: https://github.com/misxenia/SNetOC

[^3]: [https://snap.stanford.edu/data/]{}

[^4]: https://github.com/misxenia/SNetOC
