[1]Martina Boschi

[1,2]Faculty of Informatics, Università della Svizzera italiana, Via Giuseppe
Buffi, 13, Lugano, 6900, Switzerland

# Goodness of fit of relational event models

<EMAIL_ADDRESS>Ernst Jan Camiel Wit<EMAIL_ADDRESS>*

###### Abstract

A type of dynamic network involves temporally ordered interactions between
actors, where past network configurations may influence future ones. The
relational event model can be used to identify the underlying dynamics that
drive interactions among system components. Despite the rapid development of
this model over the past 15 years, an ongoing area of research revolves around
evaluating the goodness of fit of this model, especially when it incorporates
time-varying and random effects. Current methodologies often rely on comparing
observed and simulated events using specific statistics, but this can be
computationally intensive, and requires various assumptions.

We propose an additive mixed-effect relational event model estimated via case-
control sampling, and introduce a versatile framework for testing the goodness
of fit of such models using weighted martingale residuals. Our focus is on a
Kolmogorov-Smirnov type test designed to assess if covariates are accurately
modeled. Our approach can be easily extended to evaluate whether other
features of network dynamics have been appropriately incorporated into the
model. We assess the goodness of fit of various relational event models using
synthetic data to evaluate the test’s power and coverage. Furthermore, we
apply the method to a social study involving 57,791 emails sent by 159
employees of a Polish manufacturing company in 2010.

The method is implemented in the R package mgcv.

###### keywords:

relational event models, time-varying effects, random effects, generalized
additive models, goodness of fit, martingale residuals;

##  1 Introduction

Many real world processes can be described as dynamic networks, where edges
appear ordered in time. The vertices in these networks can be social actors,
biological species, geographical regions, hospital structures, or scientific
patents, depending on the context [1, 2, 3, 4]. And although the type of
relations can vary from social interactions, invasions of geographical regions
by non-native species, patient transfers between hospitals, or patent
citations, edges are characterized by time-specific interactions between these
vertices.

Over the past 15 years, Relational Event Models (REMs) have emerged as a
useful framework for studying the dynamics that drive these events. The
initial formulations of REMs within the domain of social analysis focused on
modelling how simple endogenous network statistics, such as reciprocity and
transitivity, influence the occurrence of these relational events [5]. The
counting process model formulation provided a principled mathematical
description of the relational event process [6]. Recent more flexible REM
formulations accommodate time-varying and random effects, with concerted
efforts to alleviate the computational cost associated with REM inference
techniques [7, 8]. Furthermore, [9] proposed an extension to REMs to
incorporate the collective participation of actors in a shared event, giving
rise to the concept of relational hyper-event models.

Despite the progress in relational event modelling, the contentious issue of
evaluating their fit to the data persists [10]. [11] proposed a simulation-
based approach that involves simulating a segment of the event sequence and
comparing the simulated events with the actual relational events. This method,
while comprehensive, is time-consuming, computationally expensive, and relies
on several assumptions and comparison metrics. Recently, [12] introduced
another simulation-based procedure that compares several non-modelled,
relevant statistics in both real and simulated data based on the fitted model.
The drawback of both approaches is that simulation-based approaches in
relational event modeling are computationally very intensive, as it
necessitates calculating endogenous statistics at each time point for all
potential dyads at risk of interacting. Alternatively, conventional Cox
regression approaches, including the assessment of Schoenfield, Deviance, and
Martingale residuals, have been adapted for REMs [2]. However, these
techniques lack formalization, and have not been extended to the latest REM
extensions mentioned above.

This paper will introduce a versatile technique for assessing the goodness of
fit (GOF) of arbitrary components of any relational event model. To achieve
this, we adapt a method from the literature on Cox regression [13, 14, 15, 16,
17]. Broadly speaking, this approach involves comparing an observed weighted
Martingale-type process over the time window of interest with its theoretical
behavior under a GOF assumption. The technique seeks to identify the
discrepancy between the observed statistic and its expected value under the
assumed model at each time point. This sequence is then accumulated over the
event sequence to produce the martingale-type process of interest [14].
Extending the existing literature, our method is also able to test complex
effects, including non-linear, time-varying, and random effects, as well as a
global assessment of a model’s adequacy. The technique can be generalized even
further to test if any particular feature or auxiliary statistic of the
system, not directly included in the model formulation, has been properly
accounted for by the model. Furthermore, the method is computationally
straightforward, avoiding the need for simulating relational events based on
the fitted model as required by simulation-based approaches.

The paper is structured as follows. In Section 2, we introduce a generalized
additive relational event model, including non-linear and random effects, as
well as a computationally efficient inference technique. Section 3 represents
the methodological core of the paper, where we describe the GOF method. For
expository purposes, we introduce the GOF test in the setting of a single
fixed effect, after which we show various extensions, including an omnibus
test and a brief discussion on the use of statistics that are not part of the
initial model formulation. In Section 4, we present a simulation study showing
the performance of the GOF test, evaluating its power and coverage. Finally,
in Section 5 we apply the GOF test to a study involving a sequence of 57,791
emails sent by 159 employees from January 2nd to September 30th, 2010, within
the context of a manufacturing company located in Poland [18].

The datasets generated and analysed during the current study are available in
the GOF GitHub repository here [19]. The repository also contains the R code
used for implementing these analyses.

##  2 Relational event modelling

A _relational event_ denotes an interaction occurring at time $t\geq 0$,
originating from a sender $s$ and directed towards a receiver $r$. Formally,
it is represented as a triplet $(s,r,t)$, where $s\in\mathcal{S}$ (the set of
possible senders) and $r\in\mathcal{C}$ (the set of possible receivers).
Relational events can be conceptualized as the manifestations of a _marked
point process_ , $\mathcal{M}=\\{(t_{k},(s_{k},r_{k}));k\geq 1\\}$. This
temporal process contains at times $t_{k}$ (_point_) when the interaction
$(s_{k},r_{k})$ (_mark_) occurs. A _counting process_ can be associated with
this marked point process, representing for each mark $(s,r)$ at any time
$t\geq 0$ the number of occurred interactions $(s,r)$ in $[0,t]$,

$N_{sr}(t)=\sum_{t_{k}\leq t}\mathbbm{1}_{\\{(s_{k},r_{k})=(s,r)\\}},$ (1)

In this paper, we make the following assumptions:

  1. 1.

$N_{sr}(0)=0$, meaning that no event is assumed to occur at time $0$.

  2. 2.

$N_{sr}$ is adapted with respect to the increasing filtration
$\mathcal{H}=\\{\mathcal{H}_{t}\\}_{t\geq 0}$, where $\mathcal{H}_{t}$ is a
sub-sigma field generated by the events that have occurred up to time $t$.

  3. 3.

No simultaneous events occur.

  4. 4.

$\forall t\in\mathbb{R}:N_{sr}(t)<\infty$ almost surely.

Under these assumptions, $N_{sr}$ is a right-continuous submartingale, and,
according to the Doob-Meyer theorem [6], it can be decomposed as follows,

$M_{sr}(t)=N_{sr}(t)-\Lambda_{sr}(t),$ (2)

where $M_{sr}(t)$ is a zero-mean Martingale, and
$\Lambda_{sr}(t)=\int_{0}^{t}\lambda_{sr}(u)\mathrm{d}u$ represents the
$\mathcal{H}_{t^{-}}$-measurable, predictable part of the counting process
$N_{sr}(t)$.

###  2.1 Stratified additive mixed effect REM formulation

REM involves modelling the _intensity function_ $\lambda_{sr}(t)$, measuring
the instantaneous hazard of the relational event $(s,r)$ at time $t$. In this
paper, we consider a rather general formulation of the hazard, involving a
stratified, additive, mixed effects structure. It considers a stratified
formulation of $\lambda_{sr}(t)$ that depends on $p$ potentially time-varying
covariates with fixed linear (FLE), or time-varying effects (TVE), $w$
covariates with non-linear effects (NLE), and $c$ random effects (RE),

$\displaystyle\lambda_{sr}(t)$
$\displaystyle=Y_{sr}(t)\lambda_{0g(s,r)}(t)\exp{\left[\bm{\beta}(t)^{T}\bm{x}_{sr}(t)+\bm{f}^{T}[\bm{v}_{sr}(t)]\mathbbm{1}_{w}+\bm{b}^{T}\bm{z}_{sr}(t)\right]}$
(3) $\displaystyle\bm{b}\sim\mathcal{N}(\bm{0},\Sigma(\bm{\theta}_{0}))$ (4)

In this model formulation, several key components are involved. First,
$Y_{sr}(t)$ is the at-risk indicator of event $(s,r)$ at time $t$. Then
$\lambda_{0g}(t)$ represents a non-negative _baseline intensity function_ ,
potentially specific to different population subgroups, indexed by
$g(s,r)\in\\{1,\ldots,G\\}$. Relational event applications offer a rich set of
possibilities for stratified formulations, which may involve the sender of the
relational event [6, 20], particular population subgroups [7], or the type of
event [21]. In this paper, we do not explicitly focus on stratification, and
for this reason we will use the notation $\lambda_{0}$ instead of
$\lambda_{0g(s,r)}$, but we will outline situations in which the distinction
among the two is actually relevant.

The _covariate processes_ $\bm{x}_{sr}\in\mathbb{R}^{p}$,
$\bm{v}_{sr}\in\mathbb{R}^{w}$, and $\bm{z}_{sr}\in\mathbb{R}^{c}$ are locally
bounded, left-continuous, adapted to $\mathcal{H}$, and thus predictable.
Covariates that are included in REMs may be _exogenous_ , consisting of
attributes external to the system of events, of one or both actors involved in
the relational event. Or they are _endogenous_ , incorporating the network
dynamics arising from the past. Based on their assumed impact on the intensity
function (FLE, TVE, NLE, or RE), the three sets of covariates $\bm{x}_{sr}$,
$\bm{v}_{sr}$, and $\bm{z}_{sr}$ are identified in their notation. The set of
$w$ variables included in the vector $\bm{v}_{sr}$ are assumed to have non-
linear effect, namely different values have different impacts on the hazard
function. The set of $w$ functions contained in $\bm{f}$ are assumed smooth
functions. We also include the assumption that $\mathbb{E}[f_{d}(v_{d})]$ is
equal to zero to avoid identifiability issues [22]. The main aim of this paper
is to devise a generic test for checking whether the covariates $\bm{x}_{sr}$,
$\bm{v}_{sr}$, and $\bm{z}_{sr}$ have been included in the model in a way that
is consistent with the data.

For inference purposes, we will express non-linear components as linear
combinations of some functional basis $\bm{\psi}$,

$f_{d}(v)=\sum_{l=1}^{q}\delta_{dl}\cdot\psi_{l}(v).$ (5)

Although the bases for the various covariates are not necessarily the same, we
suppress in the notation of the basis functions their dependence on the
covariate. Typical bases are the cubic spline basis, cyclic regression basis
and the P-spline basis.

We assume that $p_{0}\leq p$ covariates in $\bm{x}$ have a FLE: in this case
$\beta$ is a constant. For the remaining $p-p_{0}$ covariates $\bm{\beta}$ is
a time-varying function. Typically, time-varying effect is represented as a
smooth function of time, such as a thin plate spline, which can be expressed
in terms of a linear basis expansion [23]:

$\beta_{d}(t)=\sum_{l=1}^{q}\beta_{dl}\cdot\psi(||t-c_{l}||)\qquad\psi(t)=t^{2}\log{(t)}$
(6)

where $d>p_{0}$ and $c_{l}$ $(l=1,\ldots,q)$ consist of the set of control
points. Other bases can be selected as well. Applications of these two methods
can be found in [2] and [7], respectively.

Finally, $\bm{b}$ consists of the random effects, implemented as Gaussian
frailties. Random effects offer a versatile means of accounting for
heterogeneity in relational dynamics in the absence of specific covariates.
For instance, it can account for intrinsic heterogeneity in sender activity or
receiver popularity, instead of using viral factors, such as the sender’s out-
degree or the receiver’s in-degree [24]. Furthermore, random effects enhance
the flexibility of covariate effects, allowing for effects that may be, in a
certain amount, fixed across all statistical units but also specific according
to particular features within the system.

As previously mentioned, the model formulation in (3) is complex, yet highly
modular. Simpler model formulations can be implemented by appropriately
selecting from the aforementioned model components.

###  2.2 Case-control partial likelihood

Let $\mathcal{D}$ be a temporally ordered sequence composed of $n$ relational
events denoted as $e_{k}=(s_{k},r_{k},t_{k})$, where $k=1,\ldots,n$ according
to (3) and (4). For identifiability, we assume the process starts at time
$t=0$ and the last observation occurs at time $\tau=t_{n}$. Under assumptions
1-4 and conditional on the random effects, the partial likelihood (PL) is the
joint product of the multinomial conditional probabilities of observing events
in $\mathcal{D}$ at the time points in which the events actually occurred [5,
6]. Each of these probabilities is a ratio of the hazard of the dyad
experiencing the event at time $t_{k}$ ($k=1,\ldots,n$), and the sum of the
hazard rates for all the dyads at risk of interacting at that time, namely
those composing the _risk set_ $\mathcal{R}_{t_{k}}$. The computational cost
of PL is expensive, as the risk set scales with the number of possible dyads,
namely $n_{\mathcal{S}}\times n_{\mathcal{C}}$.

For this reason, _nested case-control_ (NCC) sampling [25] has been adapted to
REMs [26, 27], evaluating, at each time-point, a sample of the individuals at
risk, called the _sampled risk set_. The sampled risk set includes the dyad
involved in the event (the _case_) and $m-1$ sampled _controls_ (dyads that
could interact but did not). A new marked point process is thus driving the
relational events, $\\{\left(t_{k},((s_{k},r_{k}),SR_{k})\right);k\geq 1\\},$
where $SR_{k}\subset\mathcal{R}_{t_{k}}$ is the sampled risk set at time
$t_{k}$. The marked space of this new marked point process can be defined as
$E=\\{((s,r),SR):s\in\mathcal{S},r\in\mathcal{C},SR\in\mathcal{P}_{sr}\\}$.
Here, $\mathcal{P}$ is the power set of all possible pairs of senders and
receivers in $\mathcal{S}$ and $\mathcal{C}$, respectively. The set
$\mathcal{P}_{sr}\subset\mathcal{P}$ consists of of all elements in
$\mathcal{P}$ of size $m$ that include $(s,r)$. This new point process is
adapted to a new filtration $\\{\mathcal{F}_{t}\\}_{t\geq
0},\mathcal{F}_{t}=\mathcal{H}_{t}\cup\sigma(SR_{k};t_{k}\leq t)$, accounting
for previously observed sampled risk sets. Under the assumption of
_independent sampling_ , $\mathcal{F}_{t}$-intensity processes of $N_{sr}$ are
the same as $\mathcal{H}_{t}$-intensity processes. In practice, sampled
controls have the same failure risk as non-sampled ones [25]. As before, a
counting process is associated, numbering the observed mark and the associated
sampled risk set,

$N_{sr,SR}(t)=\sum_{t_{k}\leq
t}\mathbbm{1}\\{(s_{k},r_{k})=(s,r),SR_{k}=SR\\}.$

The original counting process can be retrieved as
$N_{sr}=\sum_{SR\in\mathcal{P}_{sr}}N_{sr,SR}$. The intensity process
$\lambda_{sr,SR}(t)=\lambda_{sr}(t)\cdot\pi_{t}(SR|(s,r))$ depends on the
probability of sampling the risk set $SR$ given the interacting dyad $(s,r)$.

The nested case-control sampled extension of the original relational event
data $\mathcal{D}$ is indicated as
$\mathcal{E}=\\{\left(t_{k},s_{k},r_{k},SR_{k}\right);k=1,\ldots,n\\}$, with
an associated NCC partial likelihood conditional on the random effects
$\bm{b}$,

$\mathcal{L}^{PS}(\bm{\beta},\bm{\delta},\bm{b}|\mathcal{E})=\prod_{k=1}^{n}\dfrac{\lambda_{s_{k}r_{k}}\left(t_{k}\right)\cdot\pi_{t_{k}}(SR_{k}|(s_{k},r_{k}))}{\sum_{(s,r)\in
SR_{k}}\lambda_{sr}\left(t_{k}\right)\cdot\pi_{t_{k}}(SR_{k}|(s,r))}$ (7)

where $SR_{k}$ represents a potentially complete or NCC-sampled risk set at
time $t_{k}$. When the sampled risk set coincides with the entire risk set
$SR_{k}=\mathcal{R}_{t_{k}}$, then $\mathcal{F}_{t}$ coincides with
$\mathcal{H}_{t}$.

This NCC partial likelihood formulation is conditioned on the random effects.
A traditional way to deal with random effects is by estimating the variance
component using restricted maximum likelihood estimation [28]. This involves
integrating out $\bm{\beta}$, and $\bm{\delta}$ from the nested case-control
partial likelihood. By replacing the estimates of $\bm{\theta}$,
$\bm{\delta}$, and $\bm{\beta}$, the random effects $\bm{b}$ can be predicted
based on their conditional expectations using the obtained estimates.

Alternatively, random effects can be included directly in the likelihood as
0-dimensional smooths [23]. As the random effects are assumed to follow a
normal distribution, their probability density a corresponds to a quadratic
penalization term involving an identity penalty matrix with dimensions equal
to the number of levels that the considered random factor may assume. We will
focus on this latter approach because it allows us to treat random effects
similarly to the other smooths involved in the model formulation.

###  2.3 Penalized maximum NCC partial likelihood

The spline formulation for the time-varying effects, the non-linear effects,
and the random effects enables to write the conditional NCC partial log-
likelihood in an additive form,

$\displaystyle\ell^{PS}(\bm{\gamma}|\mathcal{E})$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\log\left[\dfrac{\exp\left(\bm{\gamma}^{\top}\bm{h}_{s_{k}r_{k}}(t_{k})\right)\cdot\pi_{t_{k}}(SR_{k}|(s_{k},r_{k}))}{\sum_{(s,r)\in
SR_{k}}\exp\left(\bm{\gamma}^{\top}\bm{h}_{sr}(t_{k})\right)\cdot\pi_{t_{k}}(SR_{k}|(s,r))}\right],$

where the vector $\bm{\gamma}=(\bm{\beta},\bm{\delta},\bm{b})$ collects all
parameters, and the covariates and spline bases are collected in a vector
$\bm{h}_{sr}(t)=(\bm{h}^{(1)}_{sr}(t),\bm{h}^{(2)}_{sr}(t),\bm{h}^{(3)}_{sr}(t))$,
where

$\begin{array}[]{rcll}{h}^{(1a)}_{d,sr}(t)&=&x_{d,sr}(t)&\mbox{$x_{d}$ has a
FLE},d\leq p_{0}\\\
{h}^{(1b)}_{q(d-p_{0}-1)+l,sr}(t)&=&x_{d,sr}(t)\psi(|t-c_{l}|)&\mbox{$x_{d}$
has a TVE, $d\leq p_{0}$},l=1,\ldots,q\\\
{h}^{(2)}_{q(d-1)+l,sr}(t)&=&\psi_{l}(v_{d,sr}(t))&\mbox{$v_{d}$ has a
NLE},l=1,\ldots,q\\\ {h}^{(3)}_{d,sr}(t)&=&z_{d,sr}(t)&\mbox{$z_{d}$ has a
RE}\end{array}$

Both $\bm{\gamma}$ and $\bm{h}_{sr}(t)$ are elements of $\mathbb{R}^{P}$,
where $P$ is the sum of $p_{0}+(p-p_{0})\times q$, corresponding to the fixed
or time-varying effects, $w\times q$, corresponding to the non-linear effects,
and $c$, corresponding to the random effects. If some of the covariates have a
FLE, they do not have $q$ elements in the design matrix but only one. This
formulation specifically refers to non-stratified models. For stratified
models, the sampled non-event can be sampled from the same stratum $g(s,r)$ as
the event $(s,r,t)$. This cancels the nuisance parameter $\lambda_{0g(s,r)}$
from both numerator and denominator, resulting in the same expression
$\ell^{PS}$ above.

To impose smoothness on the time-varying and non-linear terms, as well as to
impose a normal prior on the random effect terms, we associate a _penalty
term_ with each of the smooth and random terms. For more details on penalty
terms, the reader is referred to [29, 23]. The _penalized maximum likelihood
estimator_ (PMLE) [30] maximizes the NCC log-likelihood, which is given as:

$\displaystyle\ell_{\lambda}^{PS}(\bm{\gamma})$ $\displaystyle=$
$\displaystyle\ell^{PS}(\bm{\gamma})-P_{\lambda}(\bm{\gamma}),$

where $P_{\lambda}(\bm{\gamma})=\sum_{l=1}^{L}\lambda_{l}\cdot
p_{l}(\bm{\gamma}_{\bm{i}_{l}})$, $L$ is the number of penalized terms,
$p_{l}$ is the $l$th penalty function, and $\lambda_{l}$ the $l$ tuning
parameter. Vector $\bm{i}_{l}\subseteq\\{1,...,P\\}$ contains the indices of
vector $\bm{\gamma}$ identifying the $l$th penalized parameter. To find the
PMLE $\hat{\bm{\gamma}}$, the penalized score is set equal to zero,
$\nabla\ell_{\lambda}^{PS}(\widehat{\bm{\gamma}})=\bm{0}$, resulting in an
important expression for the penalization term, evaluated at the solution
$\hat{\bm{\gamma}}$:

$\nabla\ell^{PS}(\hat{\bm{\gamma}})=\nabla P_{\lambda}(\hat{\bm{\gamma}}).$
(8)

The tuning parameters are chosen via generalized cross-validation [23].

###  2.4 REM inference as additive logistic regression

Throughout the rest of this paper, we will often refer to the scenario
involving case-control sampling with $m=2$, where we randomly sample only one
non-event for each event, i.e., $\pi_{t}(SR|(s,r))=\dfrac{1}{n(t)-1}$, where
$n(t)=|\mathcal{R}_{t}|$ and $SR\subset\mathcal{R}_{t}$, such that $(s,r)\in
SR$ and $|SR|=2$. In this case the NCC partial log-likelihood simplifies:

$\displaystyle\ell^{PS_{2}}(\bm{\gamma})$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\left[\bm{\gamma}^{T}\Delta\bm{h}_{k}-\log{\left(1+\exp{\left(\bm{\gamma}^{T}\Delta\bm{h}_{k}\right)}\right)}\right]$

where
$\Delta\bm{h}_{k}=\bm{h}_{s_{k}r_{k}}(t_{k})-\bm{h}_{s_{k}^{*}r_{k}^{*}}(t_{k})$
is the difference in the model matrix vector for the event and sampled non-
event. This corresponds to the likelihood of a logistic regression [7, 8] with
fixed response equal to $1$, without intercept term, and with covariates that
correspond to the difference of covariates of events and those of the sampled
non-event. Inference via the penalized NCC partial log-likelihood
$\ell_{\lambda}^{PS_{2}}$ allows for computationally efficient estimation of
model parameters while maintaining consistency.

##  3 Goodness of fit of relational event models

The goodness of fit framework that we will introduce here is based on the idea
that although the maximum likelihood is an _overall_ perfect fit to the data,
at various points in the temporal relational event process there may be
significant differences between observed statistics in the data and their
expected value according to the fitted relational model. In an ideal scenario,
where the model effectively captures all the relevant underlying dynamics,
these disparities are expected to fluctuate around zero. Although our proposal
is inspired by [14, 15, 16, 13, 17], our methodology extends beyond fixed
covariate assessments alone, enabling the evaluation of any generic process,
presumed to be left-continuous and adapted to the filtration of the event-
generating process, and to the augmented filtration incorporating sampling
information.

This section is structured as follows. In section 3.1, we introduce the
cumulative martingale residual process, which in its simplest form corresponds
to the cumulative score process. In section 3.2 we derive the asymptotic
distribution of this process. Subsequently, in section 3.3, we adapt this
process as a test statistic for assessing goodness of fit. Initially, we focus
on the basic case where the statistic is a single covariate with a fixed
linear effect before extending this to four additional statistical tests. In
section 3.4, we propose a test for covariates with time-varying or non-linear
effects. Section 3.5 deals with a GOF test for random effects. Section 3.6
considers a global test, concerning all components involved in the model
matrix. Finally, section 3.7 considers the case in which the statistic of
interest is not a covariate but instead consists of any feature of the
relational system. In these four scenarios, the distribution of the test
statistic is unknown, but can be simulated to produce empirical p-values.

###  3.1 Cumulative martingale residual process

Consider any $\mathcal{F}_{t^{-}}$-measurable statistic $\phi_{sr}$ of
interest. The aim is to check to what extent the statistic $\phi_{sr}$ is
captured by the current model formulation. We define a _martingale residual
process_ that captures the difference between the true $\phi_{sr}$ and its
expected value under (3),

$\mathcal{G}[\bm{\gamma},t]=\int_{0}^{t}\sum_{SR\in\mathcal{P}}\sum_{(s,r)\in
SR}\left[\phi_{sr}(u)-\dfrac{\Phi^{(0)}_{SR}[\bm{\gamma},u]}{S^{(0)}_{SR}[\bm{\gamma},u]}\right]\textmd{d}N_{sr,SR}(u)\quad
t\in[0,\tau]$ (9)

involving the following entities,

$\displaystyle\Phi^{(0)}_{SR}[\bm{\gamma},u]=\sum_{(s,r)\in
SR}\phi_{sr}(u)\cdot\exp{\left[\bm{\gamma}^{T}\bm{h}_{sr}(u)\right]}\cdot\pi_{u}(SR|(s,r)),$
$\displaystyle S^{(0)}_{SR}[\bm{\gamma},u]=\sum_{(s,r)\in
SR}\exp{\left[\bm{\gamma}^{T}\bm{h}_{sr}(u)\right]}\cdot\pi_{u}(SR|(s,r)).$

Given an event sequence $\mathcal{E}$, including the NCC sampled risk sets, we
define a weighted time-normalized observed martingale residual process for
$u\in[0,1]$ as

$\displaystyle{G}[{\bm{\gamma}},u|\mathcal{E}]$ $\displaystyle=$
$\displaystyle\sum_{k\leq\lfloor
nu\rfloor}w_{s_{k}r_{k}}(t_{k})\cdot\left[\phi_{s_{k}r_{k}}(t_{k})-\dfrac{\Phi^{(0)}_{SR_{k}}[{\bm{\gamma}},t_{k}]}{S^{(0)}_{SR_{k}}[{\bm{\gamma}},t_{k}]}\right]$
(10) $\displaystyle=$ $\displaystyle\sum_{k\leq\lfloor
nu\rfloor}w_{s_{k}r_{k}}(t_{k})\cdot\phi_{s_{k}r_{k}}(t_{k})\left[1-\nabla{\Lambda}_{s_{k}r_{k}}(t_{k})\right]$
(11)

where
$\nabla{\Lambda}_{s_{k}r_{k}}(t_{k})=\exp(\gamma^{\top}\bm{h}_{s_{k}r_{k}}(t_{k}))\cdot\pi_{t_{k}}(SR_{k}|(s_{k},r_{k}))/\sum_{(s,r)\in
SR_{k}}\exp{\left[\bm{\gamma}^{T}\bm{h}_{sr}(t_{k})\right]}\cdot\pi_{t_{k}}(SR_{k}|(s,r))$
is the $k$th increment in the cumulative intensity process due to observation
$k$. This process $G$ transforms $\mathcal{G}$ to $n$ equally spaced time-
points in $[0,1]$. The weights $w_{s_{k}r_{k}}(t_{k})\in[0,1]$ can be selected
to give the process more efficient small sample properties, by down-weighing
high-variance increments.

When the process $\phi_{sr}$ corresponds to one element of the design matrix,
namely $\phi_{sr}=h_{d,sr}$ and $w_{sr}=1$, then (10) evaluated at $u=1$
becomes

$\displaystyle{G}[{\bm{\gamma}},1|\mathcal{E}]$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{n}\left[h_{d,s_{k}r_{k}}(t_{k})-\dfrac{\sum_{(s,r)\in
SR_{k}}{h}_{d,sr}\exp{\left[{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}}{\sum_{(s,r)\in
SR_{k}}\exp{\left[{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}}\right]$
$\displaystyle=$ $\displaystyle\nabla_{d}\ell^{PS}(\bm{\gamma}),$

and is identical to the $d$th score component of the unpenalized NCC partial
log-likelihood. By substituting the maximum penalized likelihood estimate
$\widehat{\bm{\gamma}}$ and using (8), we get

$\bm{G}[\widehat{\bm{\gamma}},1|\mathcal{E}]=\nabla\ell^{PS}(\widehat{\bm{\gamma}})=\nabla
P_{\lambda}(\widehat{\bm{\gamma}}).$ (12)

We can identify several specific cases. If $\phi_{sr}$ is a covariate
$h_{d_{1},sr}$ with a FLE, its corresponding penalty term is equal to zero.
Therefore, the $d$th component of process (3.1) returns to $0$ at the end of
the observational time. This is not the case for the components of the model
matrix that refer to covariates with NLE, TVE, or RE, where the penalty term
is indeed different from zero. Furthermore, in such circumstances, each
covariate $\bm{\phi}_{sr}=\bm{h}_{\bm{i},sr}$ involves more than one component
in the design matrix, as indicated by the indices $\bm{i}$. The associated
process $\bm{G}[\hat{\bm{\gamma}},\cdot|\mathcal{E}]$ is consequently
multidimensional.

When dealing with NCC with $m=2$, the martingale residual process associated
with the covariates $\phi_{sr}=\bm{h}_{sr}$ simplifies even further,

$\bm{G}[{\bm{\gamma}},u]=\sum_{k\leq\lfloor
nu\rfloor}w_{s_{k}r_{k}}\cdot\left[1-\text{logistic}\left({\bm{\gamma}}^{T}\Delta\bm{h}_{k}\right)\right]\cdot\Delta\bm{h}_{k}.$

###  3.2 Asymptotic distribution of martingale residual process

As we want to use $G$ for testing purposes, we want to identify its
distribution under the assumption that the model fits the data. We consider
both the simpler case in which $\phi_{sr}$ is a FLE covariate, but also when
$\bm{\phi}$ is multidimensional, corresponding to TVE, NLE, and RE covariates.
Since the first case can be seen as a particular case of the second, we
introduce here the results by referring to a multivariate process
$\bm{G}[{\bm{\gamma}},\cdot|\mathcal{E}]$.

Under the true parameter $\bm{\gamma}_{0}$ and mild regularity conditions,
$n^{-\frac{1}{2}}\times\bm{G}[{\bm{\gamma}}_{0},\cdot]$ converges in
distribution to a Gaussian process with continuous paths with zero mean and
covariance function
$\mathbb{C}[\bm{G}(\bm{\gamma}_{0},t),\bm{G}(\bm{\gamma}_{0},u)]=\min{(t,u)}\cdot\bm{J}_{\bm{G}[\bm{\gamma}_{0}]},\
t,u\in[0,1]$ [31], where
$\bm{J}_{\bm{G}[\bm{\gamma}_{0}]}=\mathbb{V}(\bm{\phi}_{sr}(t)\left[1-\nabla{\Lambda}_{sr}(t|\bm{\gamma}_{0})\right])$
is defined as the variance of an individual martingale residual contribution.
Whenever $\bm{J}_{\bm{G}[\bm{\gamma}_{0}]}$ is non-singular, we define
$\bm{W}[\bm{\gamma}_{0},u]=\bm{J}_{\bm{G}[\bm{\gamma}_{0}]}^{-\frac{1}{2}}\times
n^{-\frac{1}{2}}\times\bm{G}[{\bm{\gamma}}_{0},u]$ for $u\in[0,1]$. Then,
$\bm{W}[{\gamma}_{0},\cdot]$ converges in distribution to a standard Brownian
motion $\bm{Z}(\cdot)$.

Given a consistent estimate $\widehat{\bm{\gamma}}$ and a consistent, and non-
singular covariance matrix estimate $\widehat{\bm{J}}$ for
$\bm{J}_{\bm{G}[\bm{\gamma}_{0}]}$, we have

$\widehat{\bm{W}}[\hat{\bm{\gamma}},u]=\widehat{\bm{J}}^{-\frac{1}{2}}\times
n^{-\frac{1}{2}}\times\bm{G}[\hat{\bm{\gamma}},u]\qquad u\in[0,1]$ (13)

converges to a standard Brownian bridge $\bm{Z}^{0}(\cdot)$, where
$\bm{Z}^{0}(u)=\bm{Z}(u)-u\bm{Z}(1)$.

Under some regularity conditions the penalized maximum NCC partial likelihood
estimate $\widehat{\gamma}$ is consistent, and some transformation of the
observed information matrix $\mathcal{I}(\widehat{\gamma})/n$ is a consistent
invertible estimate of the variance of process $G$ increments. This will allow
us to derive an effective test-statistic for the GOF of an FLE in section 3.3.
However, there are a number of improvements that are necessary to make
$\widehat{\bm{W}}$ a viable GOF test-statistic of TVE, NLE, and RE for finite
$n$, discussed in the subsequent sections.

###  3.3 GOF test for single covariate with fixed linear effect

In this section, we derive a GOF statistic for the case $\phi_{sr}=x_{d,sr}$,
namely $\phi_{sr}$ is a covariate that is part of the model formulation with a
FLE. We recall that in this scenario ${G}[\hat{\bm{\gamma}},u|\mathcal{E}]$
coincides with the corresponding element of the cumulative score vector, i.e.
$\nabla_{d}{\ell^{PS}(\hat{\bm{\gamma}})}$. The variance of individual
contribution to the $d$th component of the score process is given as the $d$th
diagonal element of the Fisher information matrix
$\mathcal{I}[{\bm{\gamma}}]_{d,d}$ divided by $n$, which can be estimated via
the observed information matrix,
$\widehat{J}=\frac{\mathcal{I}[\hat{\bm{\gamma}}]_{d,d}}{n}$.

Therefore,
$\widehat{W}[\hat{\bm{\gamma}},\cdot]={G}[\hat{\bm{\gamma}},\cdot|\mathcal{E}]\times\sqrt{\mathcal{I}[\hat{\bm{\gamma}}]_{d,d}}$
under the null hypothesis of a correctly specified effect converges to a
standard unidimensional Brownian bridge $Z^{0}(\cdot)$. Deviations from the
null can be tested by evaluating, for example, the largest value of
$\widehat{W}[\hat{\bm{\gamma}},\cdot]$,

$T_{x}=\sup_{u\in[0,1]}\lvert\widehat{W}[\hat{\bm{\gamma}},u]\rvert$ (14)

The distribution of the supremum of a Brownian Bridge has a Kolmogorov
distribution. Given the distribution of the univariate test statistic $T_{x}$
under the null hypothesis of interest, i.e., the adequacy of the chosen model,
we can compute exact value for the p-value when observing $T_{x}=t_{x}$,

$\mbox{p-value}(t_{x})=2\sum_{k=1}^{\infty}(-1)^{k-1}e^{-2k^{2}t_{x}^{2}}.$

###  3.4 GOF test for covariates with non-linear/time-varying effect

When testing for the goodness of fit of a covariate with a non-linear or time-
varying effect, we consider the statistics we introduced in section 2.3. Let
$\bm{i}_{d}$ correspond to the indices in the model matrix that represent the
$d$th covariate of interest. In particular, for time-varying covariates we
included the $q$-dimensional vector
$\bm{h}_{\bm{i}_{d},sr}(t)=\bm{\psi}(t)x_{d,sr}(t)$ in the model matrix,
whereas for non-linear effects this corresponded to
$\bm{h}_{\bm{i}_{d},sr}(t)=\bm{\psi}(v_{d,sr}(t))$.

We consider the multivariate martingale residual process
$\bm{G}[\hat{\bm{\gamma}},\cdot]$ for $\bm{\phi}_{sr}=\bm{h}_{\bm{i},sr}$. For
finite $n$, the process $\bm{G}[\hat{\bm{\gamma}},\cdot]$ arrives at
$\nabla_{\bm{i}}\ell^{PS}(\hat{\bm{\gamma}})$. To improve the Brownian bridge
approximation, we define a centered process by subtracting, for each term, a
fraction of the penalty term,

$\bm{G}_{k}[\hat{\bm{\gamma}},t_{k}]=\left[\bm{h}_{\bm{i},s_{k}r_{k}}(t_{k})-\dfrac{\sum_{(s,r)\in
SR_{k}}\bm{h}_{\bm{i},sr}\exp{\left[\hat{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}}{\sum_{(s,r)\in
SR_{k}}\exp{\left[\hat{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}}\right]-\dfrac{\nabla_{\bm{i}}\ell^{PS}(\hat{\bm{\gamma}})}{n}$
(15)

Since $\bm{J}$ is defined as the variance of an individual contribution
$\bm{G}_{k}[{\bm{\gamma}}_{0}]$, a candidate for its estimation is the
empirical variance-covariance matrix at the end of the observational period,
namely
$\widehat{\bm{J}}=n^{-1}\times\sum_{k=1}^{n}{{\bm{G}}}_{k}[\hat{\bm{\gamma}},t_{k}]{{\bm{G}}}_{k}[\hat{\bm{\gamma}},t_{k}]^{T}$
[31]. These two elements now allow us to define a normalized process for
$u\in[0,1]$ if the matrix $\widehat{\bm{J}}$ is invertible,

$\widehat{\bm{W}}[\hat{\bm{\gamma}},u]=\widehat{\bm{J}}^{-\frac{1}{2}}\times
n^{-\frac{1}{2}}\times\sum_{k\leq\lfloor
nu\rfloor}\bm{G}_{k}[\hat{\bm{\gamma}},t_{k}|\mathcal{E}].$ (16)

This process starts and ends at the origin for each value of $n$, and
converges to a $q$-dimensional Brownian bridge
$\bm{Z}^{0}=\begin{bmatrix}Z_{1}^{0}(\cdot)&\ldots&Z_{q}^{0}(\cdot)\end{bmatrix}$.
Following [32], we define the following test statistic,

$T_{\psi}=\sup_{u\in[0,1]}\lVert\widehat{\bm{W}}[\hat{\bm{\gamma}},u]\rVert^{2}$
(17)

There is no closed form for the theoretical distribution of the supremum of
the squared norm of a multivariate Brownian bridge,
$\sup_{u\in[0,1]}\lVert{{\bm{Z}}}^{0}(u)\rVert^{2}$, but it can easily be
simulated and the empirical p-value can be computed as the fraction of
statistics that are larger or equal to the observed value in (17).

###  3.5 GOF test for random effects

As we interpret random effects as $0$-dimensional splines, evaluating the
goodness of fit follows a similar approach to that outlined in section 3.4.
However, some special considerations are necessary. Typically, random effects
correspond to a set of entries $\bm{i}$ of the design matrix. By identifying
the statistic of interest $\bm{h}_{\bm{i}_{d},sr}=\bm{z}_{d}$, where
$\bm{z}_{d}$ is the $d$th random covariate vector, we can calculate the same
contributions $\bm{G}_{k}[\hat{\bm{\gamma}},t_{k}]$ as defined in (15). These
vectors can be combined, as in (16), to obtain a $|\bm{i}_{d}|$ dimensional
process $\widehat{\bm{W}}[\hat{\bm{\gamma}},\cdot]$. Again, similar to (17),
we can define the test-statistic $T_{z}$, which under the null would be the
supremum of the squared norm or a standard $|\bm{i}|$ dimensional Brownian
bridge.

Due to the robustness of random effect estimation, the GOF test tends to have
low power for detecting minor distributional deviations. Still, as we will
show in section 4, certain forms of misspecification, especially those
involving temporal misspecification of $\bm{z}_{sr}(t)$, have adequate power
and can be actually detected in moderate data settings.

###  3.6 Global GOF test of the relational event model

The versatility of the cumulative martingale residual approach allows for its
extension to a comprehensive _omnibus test_. This test aims to confirm whether
the model specification is appropriate _overall_. As such, one would perform
this test first after fitting a model. In case it gets rejected, one would
proceed to the individual tests discussed before in order to diagnose exactly
what part of the model has been misspecified.

We consider a collection of $L$ covariates with a combination of possible
effects (FLE, TVE, NLE or RE), each associated with a set of $\bm{i}_{l}$
columns in the model matrix ($l=1,\ldots,L$). For each covariate $l$ we can
define the $|\bm{i}_{l}|$-dimensional normalized process
$\widehat{\bm{W}}[\hat{\bm{\gamma}},\cdot]$, the associated test statistics
$T_{l}$ and p-values $P_{l}$. There is an extensive literature on how to
combine p-values. Given that the p-values are correlated, we consider the
Cauchy combination test [33], which is a powerful test under arbitrary
dependency structures that has the additional advantage of having an analytic
p-value calculation. The global omnibus test-statistic is defined as

$T_{o}=\frac{1}{L}\sum_{l=1}^{L}\tan\left(\pi(0.5-P_{l})\right).$ (18)

Given the observed test-statistic $T_{o}=t_{o}$, the p-value for this global
test can be calculated explicitly as,

$\mbox{p-value}(t_{o})=1/2-\arctan(t_{o})/\pi.$

###  3.7 GOF test for auxiliary statistics

A model attains additional epistemological value when it transcends its
immediate formulation by predicting phenomena not explicitly included within
its original parametrisation. This emergent predictive power suggests a deeper
correspondence between the model and the complexities of the real world
phenomenon, indicating that the model captures essential truths about the
system it represents. For this reason, in addition to considering testing for
the GOF of model components, in this section we consider testing for arbitrary
auxiliary statistics $\phi$. For instance, we can consider using auxiliary
statistics as discussed by [12]. These are relational summary statistics of
the data, such as higher-order interactions, that provides a relevant
complementary perspective and may vary depending on the context.
Alternatively, we can consider a possible stratification variable of the
actors or events involved in the system. We would then be interested in
determining whether this stratification has been adequately accounted for by
the model definition that does not explicitly include this stratification.

Given a univariate auxiliary statistic $\phi_{sr}$, we define the process
$G[\hat{\bm{\gamma}},\cdot]$ according to (10) and an associated test-
statistic,

$T_{\phi}=\sup_{u\in[0,1]}|G[\hat{\bm{\gamma}},u]|$

The main challenge is to find the distribution of $T_{\phi}$ under the null
hypothesis. In fact, when $\phi$ is not a model component, the distributional
properties outlined in section 3.2, related to the score process, are lost. If
the model is correctly specified, [17] suggests approximating the distribution
of $G[\hat{\bm{\gamma}},u|\mathcal{E}]$ by simulating multiple replicates of
the process $G^{*}[\hat{\bm{\gamma}},u|\mathcal{E}]$, which only involves
sampling i.i.d. standard normal distributed $N_{k},k=1,\ldots,n$,

$\displaystyle G^{*}[\hat{\bm{\gamma}},u|\mathcal{E}]=\sum_{k\leq\lfloor
nu\rfloor}\left[\phi_{s_{k}r_{k}}(t_{k})-\dfrac{\Phi^{(0)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]}{S^{(0)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]}\right]\cdot
N_{k}$ (19) $\displaystyle-\sum_{k\leq\lfloor
nu\rfloor}\dfrac{\Phi^{(0)*}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]^{T}}{S^{(0)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]}\times\mathcal{I}[\hat{\bm{\gamma}}]^{-1}\cdot\sum_{k=1}^{n}\left[\bm{h}_{s_{k}r_{k}}(t_{k})-\mathbb{E}^{U}_{SR_{k}}\left[\hat{\bm{\gamma}},t_{k}\right]\right]\cdot
N_{k}$

where,

$\displaystyle\Phi^{(0)*}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]=\sum_{(s,r)\in
SR_{k}}\left(\bm{h}_{sr}(t_{k})-\mathbb{E}^{U}_{SR_{k}}\left[\hat{\bm{\gamma}},t_{k}\right]\right)\cdot\phi_{sr}(t_{k})\cdot\exp{\left[\hat{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}\cdot\pi_{t_{k}}(SR_{k}|(s,r))$
(20) $\displaystyle S^{(1)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]=\sum_{(s,r)\in
SR_{k}}\bm{h}_{sr}(t_{k})\cdot\exp{\left[\hat{\bm{\gamma}}^{T}\bm{h}_{sr}(t_{k})\right]}\cdot\pi_{t_{k}}(SR_{k}|(s,r))$
$\displaystyle\mathbb{E}^{U}_{SR_{k}}\left[\hat{\bm{\gamma}},t_{k}\right]=\frac{S^{(1)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]}{S^{(0)}_{SR_{k}}[\hat{\bm{\gamma}},t_{k}]}.$

$\mathbb{E}^{U}_{SR}\left[\hat{\bm{\gamma}},u\right]$ corresponds to the
expectation of the process $\bm{h}_{sr}$ under the estimated model. The
expression
$\sum_{k=1}^{n}\left[\bm{h}_{s_{k}r_{k}}(t_{k})-\mathbb{E}^{U}_{SR_{k}}\left[\hat{\bm{\gamma}},t_{k}\right]\right]$
is another way of writing the score (12). This approach has been shown to
converge to the distribution of $G[\hat{\bm{\gamma}}]$.

Given $b=1,\ldots,B$ samples $G_{b}^{*}[\hat{\bm{\gamma}},\cdot]$ of the
simulated process along with their associated test-statistic $t^{*}_{\phi,b}$,
we can obtain an empirical p-value for $T_{\phi}=t_{\phi}$ by adapting the
statistical testing procedure from [17] as

$\widehat{\mbox{p-value}}(t_{\phi})=\frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}\\{t_{\phi,b}^{*}\geq
t_{\phi}\\}.$

Alternatively, under the assumptions stated in [14], it is possible to prove
that $n^{-\frac{1}{2}}\times G[\hat{\bm{\gamma}},\cdot]$, even in its general
formulation, weakly converges in distribution to a zero-mean Gaussian process,
which leads to the possibility of a formal $\chi^{2}$ statistical test by
using $\widehat{W}[\hat{\bm{\gamma}},1]$.

##  4 Simulation studies

This section provides an empirical evaluation of the goodness-of-fit tests
presented in section 3. We will examine both coverage and power, two key
factors essential for determining the effectiveness of a statistical test. To
achieve this, we will investigate several scenarios, explaining the test’s
adaptability and issues in each.

Relational event occurrences are simulated according to a marked counting
process described in (3). Technical details can be found in the Supplementary
Materials. We consider a data generation process (DGP), which may involve
linear, time-varying, non-linear, and random effects of various factors. Each
simulation setting is repeated $n_{\textmd{sim}}=500$ times. In order to study
coverage (section 4.1) and power (section 4.2) of the testing procedure, we
examine situations where the fitted model is correctly specified and where it
is not, respectively. In section 4.3, we study the omnibus test to evaluate
the relational event model in its entirety.

###  4.1 Coverage of goodness-of-fit test

When the model is correctly specified, we expect that p-values of the GOF test
are uniformly distributed between $0$ and $1$. In particular, the proportion
of rejections obtained by a correct statistical testing procedure is
approximately equal to the significance level, which, for simplicity, is set
equal to $\alpha=0.05$ in the simulations.

####  4.1.1 Coverage of NLE GOF test with increasing sample sizes

In the first simulation study, we investigate how the coverage of the GOF test
for non-linear effects evolves with increasing sample size,
$n=1000,5000,10000,$ and $50000$. Relational event models capture dynamics
like reciprocity, common in real-world interactions such as email exchanges or
trading favours. In our simulation study, we focus on reciprocity as a driving
factor. The DGP incorporates reciprocity in an intrinsic non-linear way,
namely as a non-linear function of the time since the last reciprocal event,
which allows for a higher hazard closer to this event, which reduces over
time. We estimate the non-linear effect using a $10$-dimensional thin plate
regression spline, which smoothly varies the reciprocity effect over time,
reflecting the flexibility that characterizes the DGP.

|
---|---
(a) | (b)
Figure 1: Coverage of the statistical procedure. _Left_. Empirical p-value
distribution under the null hypothesis of correctly incorporating a non-linear
reciprocity effect for different sample sizes ($n=1000,5000,10000,50000$) with
proposed GOF test in (17). Reciprocity has been fitted by employing a
$9$-dimensional thin plate regression spline. The test achieves almost perfect
coverage over a wide range of sample sizes. _Right_. Empirical p-value
distribution under the null hypothesis of correctly incorporating a random
sender effect for different numbers of actors
($n_{\mathcal{S}}=n_{\mathcal{C}}=10,50,100,150$) with the proposed GOF test
in section 3.5. The test tends to be somewhat conservative for a low number of
actors, whereas they appear more uniform when the number of actors is larger.

For each simulated dataset, we obtain a model fit and evaluate the
multivariate observed martingale residual process
$\bm{G}[\hat{\bm{\gamma}},\cdot]$, in (15) and calculate test statistics
following (17). Since the model includes only reciprocity, the test involves
the entire score vector. The empirical $p$-values are found by simulating
$9$-dimensional Brownian bridges, reproducing the theoretical behaviour of the
observed score process. Figure 1a shows the distribution of $p$-values across
$n_{\textmd{sim}}=500$ simulations for different sample sizes, relative to the
expected black line representing the $U(0,1)$ distribution. We can conclude
that the GOF test coverage performance is robust across different sample
sizes. This is confirmed proportion of rejections that is approaching the
expected value of $0.05$.

####  4.1.2 Coverage of RE GOF test with increasing number of actors

In this simulation study, we investigate how the coverage of the GOF test for
random effects proposed in section refsubsubsec:random-effects evolves with an
increasing number of actors, $p=10,50,100$, and $150$. In the simulation,
events are exclusively driven by the sender’s intrinsic properties. In other
terms, the event generation is driven by an intensity function that only
depends on a sender-specific parameter, that has been sampled from a Gaussian
distribution. For each simulated data, we apply a model incorporating sender-
related random effects, fitted as $0$-dimensional splines. This model is
correctly specified, so we expect the empirical $p$-values to be distributed
uniformly.

As Figure 1b shows, the goodness of fit test for random effects seems to be
slightly too conservative when there is only a small number of actors in the
system. As soon as $p=n_{\mathcal{S}}=n_{\mathcal{C}}$ increases, the GOF test
achieves nominal coverage, as indicated by the fact that the distribution
converges to a $\mathcal{U}(0,1)$ distribution.

###  4.2 Power of goodness-of-fit test

The power of a statistical test is its ability to correctly reject a null
hypothesis when it is false. The power of the test can be evaluated by
applying it to a misspecified model.

####  4.2.1 Detecting FLE misspecification with increasing sample size

In this section, we evaluate how the GOF test can detect a misspecified linear
effect as an increasing function of the sample size, $n=1000,5000,10000,$ and
$50000$. As in the previous section, we simulate a process with a non-linear
function of reciprocity, defined as the inter-arrival time since the last
reciprocal event occurred. In this section, we fit a misspecified model by
assuming a linear effect for reciprocity.

|
---|---
(a) | (b)
Figure 2: Power of the statistical procedure.. _Left_. Incorporating
reciprocity linearly. For each size of the experiment
($n=1000,5000,10000,50000$), the picture reports the empirical distribution of
empirical $p$-values. Since we test a covariate with a fixed effect, we can
rely on (14). The power increases as the size increases. _Right_. Events are
driven by sender intrinsic properties. In this case, for each number of actors
($n_{\mathcal{S}}=n_{\mathcal{C}}=10,50,100,150$), the model inadequately
incorporates a different impact of time of the event according to the sender
involved. Test follows (17). Although the misspecification is clearly
detected, the power tends to reduce as the number of entities increases.

Figure 2a shows a clear increase in power as the number of events increases.
This allows us to conclude that the test becomes more reliable and robust with
larger sample sizes. There is another issue that exacerbates this point: as
the event sequence increases, the time-sequence lengthens, which will affect
for a number of events the time since the last reciprocal event. This means
that the linear effect for reciprocity simply becomes a worse approximation
with increasing $n$. This further necessitates using a smooth function capable
of accommodating the actual non-linear effect.

####  4.2.2 Detecting RE misspecification with increasing actor numbers

In this section, we test the ability of the goodness-of-fit test for random
effects to detect misspecification for increasing network sizes,
$p=10,50,100$, and $150$. Again, as in the previous section, the simulation
consists of a random intercept model for each of the senders of the events.
Instead, the fitted model considers a random slope model for each of the
senders as a function of time, incorrectly assuming that the impact of the
interactions is affected by the time of the interaction.

Figure 2b shows that the misspecification in this case is clearly detected.
Perhaps somewhat surprisingly, as the number of actors increases there is a
reduction in the power of the test. What is happening seems to be that as the
number of relational events $n$ is kept fixed in the simulations, the increase
in the number of actors reduces the overall accuracy of the model fit and
therefore the power of the test. Therefore, it may be that the power reduction
is due to the increase in the effective number of parameters in the model,
rather than to the misspecification itself.

The conducted simulation studies show that our technique has high power and
flexibility, identifying several forms of misspecification. The result is
confirmed in the following section, where all these components are considered
in an omnibus test.

###  4.3 Global testing of a relational event model

We conclude with an example of the application of the global goodness-of-fit
test. The omnibus test introduced in section 3.6 aims at confirming whether or
not the model can be deemed to have an appropriate fit to the data _overall_.
In this simulation study, we rely on $4$ DGPs. The DGPs progressively increase
in complexity, by adding one of the following additional covariates:

  1. 1.

Time-based covariate for reciprocity; in this case, the reciprocity measured
as defined above is assumed to have a linear impact on the intensity function;

  2. 2.

Exogenous covariate, sampled from an exponential distribution, with a linear
effect;

  3. 3.

Exogenous covariate, sampled from a Gaussian distribution, with a time-varying
effect;

  4. 4.

Random intercept sender effect.

So, the first DGP incorporates only reciprocity, while the last one includes
all the four covariates mentioned above.

To evaluate coverage, we fit the correctly specified model for each of these
four datasets and, as before, we look to the empirical distribution of the
$p$-values. In this scenario, we rely on the statistical test-statistic
$T_{o}$. We note that the first two covariates are uni-dimensional, the third
one is a $9$-dimensional thin plate regression spline, and the fourth
covariate is $50$-dimensional random factor.

|
---|---
(a) | (b)
Figure 3: Omnibus statistical procedure.. _Left._ For each model including a
different number of groups of elements ($L=1,2,3,4$), we show the empirical
distribution of the simulated p-values related to the test of the overall
adequacy of the relational event model according to the test in (18). _Right._
In the scenario in which $2$ elements are incorporated in the model
formulation, four model are tested: the correctly specified (CS) and those
with one or both elements misspecified (MS).

Figure 3a confirms the reliability of the omnibus test, even if we outline
that, especially when including the random effect, the test tends to be
slightly conservative. We think that this may be because the test makes no
assumptions about the dependence between the statistics. Since it is quite
common to deal with correlated statistics in REMs, this lack of assumption is
relevant even if it comes at the price of a slightly conservativeness of the
test. Figure 3b shows four possibilities for a model with $L=2$ covariates:
both covariates are correctly specified, one is incorrectly specified or both
are incorrectly specified. The completely misspecified model is most easily
identified, whereas the models with one misspecified element have a lower
power. The Supplementary materials report all other graphical results.

The purpose of the global test is that once a model has been identified as
misspecified, one can go further and investigate the source of the problem.
For example, in the scenario in which $2$ elements are incorporated in the
model formulation, it might be that only one of them or both of them are not
modeled properly. In the Supplementary Materials, the reader can find for each
of the fitted models both the global and the group-specific tests of adequacy.

##  5 Email communication in a manufacturing company

This section examines the email exchanges among employees of a manufacturing
company located in Poland [18]. This dataset has previously been analyzed to
investigate the temporal component of social dynamics [21]. The original
dataset, detailed in [34], comprises over 82,000 emails sent between January
2nd, 2010, and September 30th, 2010. After a preliminary analysis aimed at
removing duplicates and self-sent emails, the final dataset consists of
$n_{\mathcal{E}}=57791$ relational events involving
$n_{\mathcal{S}}=n_{\mathcal{C}}=159$ employees.

###  5.1 Incorporating temporal structure in relational dynamics

[21] focused their analysis on the temporal aspects of reciprocity,
repetition, and cyclic closure. Here, we additionally include transitive
closure. Reciprocity refers to the mutual exchange of goods, services, or
behaviours between actors, relying on a give and take principle that
encourages collaboration among the parties involved. In the context of email
communication, this typically manifests as a tendency to respond courteously
to emails. Repetition reflects the tendency to maintain established patterns
in communication. In email exchanges, this might involve following up with
previously contacted counterparts. Cyclic closure and _transitive closure_ are
examples of triadic closure, which involves an event between two actors that
have previously been in touch with a third actor. In cyclic closure the event
$(s,r)$ follows the events $(r,k)$ and $(k,s)$, whereas transitive closure it
follows the events $(s,k)$ and $(k,r)$. Triadic closure is known for its role
in path abbreviation to enable contact between previously unconnected
individuals.

Once the endogenous covariates of interest are defined, the next concern is
their quantification. Table 1 outlines two potential approaches for
quantifying each endogenous covariate, highlighting different choices for
assessing the presence and magnitude of these phenomena. The initial approach
involves a simple identity function, determining whether the event of interest
occurred before the current time. The second method incorporates a temporal
dimension by analyzing the time interval between the occurrences of current
time and the time of the relevant event in the past. In order to deal with the
possibility that the event may not have occurred at all, we generalize the
expression for the time-interval as follows:

$\Delta\tau_{sr}^{*}=\exp{[-b\times(t-t^{*}(s,r))]}$ (21)

where $t^{*}(s,r)$ represents the time of the event of interest. For example,
for reciprocity this could be the last event $(r,s)$ prior to $t$, i.e.,
$t^{*}(s,r)=\max_{s_{k}=r,r_{k}=s,t_{k}<t}t_{k}$. Practically, it can be
computed even when the event of interest has not occurred, yielding a value of
$0$. The tuning parameter $b$ allows for adjustment of the time difference
values to improve uniformity of the model matrix elements in the interval
$[0,1]$. Although the time-scale may not be immediately interpretable, the
effect can be easily transformed to the original time-scale.

Furthermore, several definitions and conceptualizations of these endogenous
covariates can be proposed and must then be chosen and evaluated: first,
considering which events of interest to include. For instance, when multiple
reciprocal events exist in the past, we might choose the first one, the last
one, or the average time of the occurred ones; in this work, we refer to the
most recent one. Additionally, we may include a lifespan in the decay function
of the time window. Finally, we need to decide about the closure of the
mentioned dynamics: for example, in our conceptualization, responding to a
reciprocal event makes that event no longer reciprocal. This is a choice and
can be potentially be evaluated by means of a goodness-of-fit test. Some
potential and assumed choices are described in Table 1.

Relational Dynamics Interpretation Representation Identity Covariate Time-
Related Covariate Reciprocity If the contribution of reciprocity is positively
estimated, then there is a tendency to send emails to counterparts who have
previously initiated communication. $s$$r$$s$$r$$t$$t_{rs}^{l}$
$t^{l}_{sr}<t^{l}_{rs}<t$ $\displaystyle\textmd{rec}^{\textmd{id}}_{sr}(t)=$
$\displaystyle 1_{\left[t>t^{l}_{rs}\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{rs}<t$ $\displaystyle\textmd{rec}^{\textmd{time}}_{sr}(t)=$
$\displaystyle\exp{\left[-\left(t-t^{l}_{rs}\right)\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{rs}<t$ Repetition If the contribution of repetition is
positively estimated, then there is a tendency to send emails to counterparts
already contacted before. $s$$r$$s$$r$$t$$t_{sr}^{l}$ $t^{l}_{sr}<t$
$\displaystyle\textmd{rep}^{\textmd{id}}_{sr}(t)=$ $\displaystyle
1_{\left[t>t^{l}_{sr}\right]}$ $\displaystyle t^{l}_{sr}<t$
$\displaystyle\textmd{rep}^{\textmd{time}}_{sr}(t)=$
$\displaystyle\exp{\left[-\left(t-t^{l}_{sr}\right)\right]}$ $\displaystyle
t^{l}_{sr}<t$ Cyclic Closure If the contribution of cyclic closure is
positively estimated, it indicates a tendency for a sender to correspond with
counterparts who have previously initiated communication with another
employee, and subsequently, that employee has communicated with the sender.
$s$$r$$k$$t$$t_{rk}^{l}$$t_{ks}^{l}$ $t^{l}_{sr}<t^{l}_{rk}<t^{l}_{ks}<t$
$\displaystyle\textmd{cyc}^{\textmd{id}}_{sr}(t)=$ $\displaystyle
1_{\left[t>t^{l}_{ks}\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{rk}<t^{l}_{ks}<t$ $\displaystyle k\text{: common
counterpart}$ $\displaystyle\textmd{cyc}^{\textmd{time}}_{sr}(t)=$
$\displaystyle\exp{\left[-\left(t-t^{l}_{ks}\right)\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{rk}<t^{l}_{ks}<t$ $\displaystyle k\text{: common
counterpart}$ Transitive Closure If the contribution of transitive closure is
positively estimated, it indicates a tendency for a sender to directly
correspond with counterparts, previously indirectly contacted by means of an
intermediator. $s$$r$$k$$t$$t_{sk}^{l}$$t_{kr}^{l}$
$t^{l}_{sr}<t^{l}_{sk}<t^{l}_{kr}<t$
$\displaystyle\textmd{trs}^{\textmd{id}}_{sr}(t)=$ $\displaystyle
1_{\left[t>t^{l}_{kr}\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{sk}<t^{l}_{kr}<t$ $\displaystyle k\text{: common
counterpart}$ $\displaystyle\textmd{trs}^{\textmd{time}}_{sr}(t)=$
$\displaystyle\exp{\left[-\left(t-t^{l}_{kr}\right)\right]}$ $\displaystyle
t^{l}_{sr}<t^{l}_{sk}<t^{l}_{kr}<t$ $\displaystyle k\text{: common
counterpart}$

Table 1: Table of Relational Dynamics

###  5.2 GOF testing of temporal dynamics of email communication

The goal of this section is to determine whether and how our goodness-of-fit
test allows us to understand if some of these choices concerning the effect
structure of the covariates may be more or less adequate for fitting the REM
to the data at hand. We adopt the notation as outlined in Table 1. Our first
objective is to assess the optimal relational event model that incorporates
elements of reciprocity, repetition, cyclic, and transitive closure.

Suppose, for a moment, that the precise timestamp of the emails is unknown,
and only their sequential order is available. Such a scenario is often
encountered in practical situations [5]. In this case, it remains feasible to
assess the endogenous dynamics for all four covariates using an identity
function. The fitting procedure leads us to positive estimates for
reciprocity, repetition and transitive closure, while a negative one for
cyclic closure. The global test (18) is rejected, which leads us to evaluate
whether each of the four covariates have been properly modeled, relying on the
computation of the one-dimensional statistics in (14). For each of the effects
the p-values are all smaller than $\alpha=0.05$, suggesting that all elements
in the model are incorrectly included.

Within this context, we have however access to the timestamps of the events.
Consequently, we may choose to evaluate more sophisticated versions of the
covariates, specifically those integrating temporal information. Potential
model formulations become more diverse, even under the assumption that each
relational dynamic is exclusively modelled either as an identity or as a
linear function of time. Model selection can be conducted by assessing the
AIC, as outlined by [35]. The best model in terms of AIC is given by
$\lambda_{sr}(t)=Y_{sr}(t)\times\lambda_{0}(t)\times\exp{\left[f_{sr}(t)\right]}$,
where

$f_{sr}(t)=\beta_{\textmd{rec}}\textmd{rec}^{\textmd{time}}_{sr}(t)+\beta_{\textmd{rep}}\textmd{rep}^{\textmd{id}}_{sr}(t)+\beta_{\textmd{trs}}\textmd{trs}^{\textmd{time}}_{sr}(t)+\beta_{\textmd{cyc}}\textmd{cyc}^{\textmd{id}}_{sr}(t)$
(22)

However, the global and individual GOF tests still reject the overall model,
as well as each term individually at the $\alpha=0.05$ level. We thus need to
look for an alternative and richer model formulations with the information at
our disposal.

Following [24], it might be relevant to include in the model formulation
random effects accounting for the sender activity and the receiver popularity.
In this way, we can effectively model the individual heterogeneity to send (or
receive) emails. Furthermore, we consider a non-linear effect for all four
temporal versions of the endogenous covariates, i.e.,

$f_{sr}(t)=f_{1}(\textmd{rec}^{\textmd{time}}_{sr}(t))+f_{2}(\textmd{rep}^{\textmd{time}}_{sr}(t))+f_{3}(\textmd{trs}^{\textmd{time}}_{sr}(t))+f_{4}(\textmd{cyc}^{\textmd{time}}_{sr}(t))+b^{\textmd{act}}_{s}+b^{\textmd{pop}}_{r}$
(23)

Each of the endogenous variables are fitted as thin plate regression splines.
The random effects are tested as non-properly modeled, but their inclusion
leads to an improvement of the four covariates that are globally tested as
adequate (p-value = $0.07$). This is consistent with the discussion in [24].
Figure 4a reports the GOF curves $\hat{\bm{W}}[\hat{\bm{\gamma}},u]$ for the
individual model components over the normalized time $u\in[0,1]$. The
individual test for repetition has a rather low p-value ($0.019$), suggesting
that the fit is not optimal. We have not optimized the temporal scale
parameter in (21), and simply selected $b=1$. By adapting this the model could
possibly be further improved.

|
---|---
(a) | (b)
Figure 4: _Left_. GOF process $\hat{\bm{W}}[\hat{\bm{\gamma}},\cdot]$
evaluated for each of the four endogenous covariates included in model (23)
with non-linear effects. Although the global test does not find any major
evidence against the adequacy of the model, it might be that repetition is
still not fully accurately incorporated in the model. It might be possible by
changing the time-scale parameter $b$ in (21). _Right_. Non-linear impact of
reciprocity on the intensity of email occurrence. Reciprocity is represented
here by the corresponding inter-arrival times of emails. Emails have a
tendency to be answered quickly, but this is followed by a steep decline.
There is another peak after 16 hours, probably corresponding to the morning
after receiving the email the day before. To enhance readability, the
background is shaded according to working hour blocks, aligning with typical
working and resting hours.

###  5.3 Final model interpretation

Model interpretation is not the core focus of this paper. Nevertheless, the
leading idea of this work is that one should exclusively interpret those model
formulations that are assessed as having an adequate model fit. In this case,
we focus on the interpretation of the model in (23), even if an improvement is
still required in terms of the modeling of random effects. Specifically, we
focus here on reciprocity. The other non-linear endogenous effects can be
found in the Supplementary Materials.

Figure 4b depicts the estimated function for reciprocity plotted against the
time from the last observed reciprocal event. We highlight time from last
reciprocal event in a number of intervals. It can be seen that once an email
is sent, the likelihood of receiving an answer decreases quickly as time
passes during the remainder of the same workday (around $8$ hours). After this
minimum, the likelihood of receiving an answer has another peak at $16$ hours
after the email, likely because people postponed the response to the following
day. A similar pattern can be observed after two workdays following the email,
although the increase after the minimum is smaller in absolute value. Around
$64$ hours after the email, the chance of receiving an answer continues to
decrease without any further increase.

The model includes random sender and receiver effects. The estimated standard
deviations are $\hat{\sigma}_{\textmd{act}}=0.97$ and
$\hat{\sigma}_{\textmd{pop}}=1.54$, suggesting that particularly the receiver
has a strong influence on whether a particular event happens. It would be
interesting to cross-list this sender and receiver heterogeneity with some
kind of exogenous information concerning the actors involved. Unfortunately,
information reported in [34] does not provide obvious clues about the
interpretation of the random effects. Furthermore, the rejection of the null
of GOF tests for these effects does not allow us to interpret them.

While acknowledging that our final model is not exhaustive, we have presented
a computationally efficient method for assessing the adequacy of several
potential relational event models that may be fitted to these relational data.
There remains the possibility of achieving a more adequate model formulation
by appropriately incorporating additional causal determinants.

##  6 Discussion

As has been clear from the empirical analysis there is an intricate
relationship between model fit and model adequacy. In the context of assessing
the goodness of fit for relational event models, some literature still relies
on likelihood-based information criteria, such as AIC or BIC. These methods,
while useful for selecting the most predictive or best fitting model from a
set of potential candidates, primarily focus on explained variation rather
than the adequacy of fit [36, 24].

Model selection is a crucial step in relational event modeling, where the
potential dynamics influencing event occurrences can be numerous and complex.
For instance, [20] highlights a significant number of statistical measures
that can be incorporated into a REM. Likelihood-based methods can be employed
to compare these various options effectively. Additionally, model selection is
essential for evaluating the importance of including random effects, which aim
to explain the intrinsic characteristics of the actors or event types. This
aspect is particularly emphasized in [24]. Nevertheless, selecting the best
model based on some information criteria or feature selection method does not
guarantee that the model is adequate. We provided an example of this
phenomenon in our application section. Model adequacy is crucial because a
model that is not adequate cannot be interpreted.

In practice, any goodness of fit procedure should commence from the global
omnibus test. If this test is rejected, the next step is to inspect the
individual model components. This inspection can identify where the problem in
the model formulation occurred. The temporal structure of the martingale
residual process gives the analyst a lead to identify where the problem
occurs. In this case, they can examine the specific time frame where the
problem arises and speculate on potential causes and improvements.

##  7 Conclusions

This paper has introduced a comprehensive method for assessing the goodness of
fit of relational event models, including their individual components. We
propose a procedure that can be easily applied after fitting the relational
event model, without requiring additional assumptions about what needs to be
tested. This approach focuses on evaluating the actual components of the model
formulation. Furthermore, the method can be extended to test for the GOF of
arbitrary auxiliary statistics. Therefore, we believe our method can be
directly compared to that described by [12]. Our approach, however, offers the
possibility of not needing to simulate the entire DGP according to the fitted
relational event model, as it is based on computationally more efficient
approaches.

Indeed, the advantage of our method, compared to existing GOF techniques, lies
in its low computational overhead, particularly when testing a component
included in the model formulation. Generally speaking, the primary
computational cost involves simulating the theoretical asymptotic behaviour of
the process. This is significantly less demanding than the computational cost
required to simulate events such as in [11] or [12], where the simulation of
events requires the computation and the updating of endogenous covariates for
all potential events (proportional to the square of the number of actors) at
each time point of interest.

When one does not reject the GOF of a component of a model, we will usually
pretend that the corresponding element is properly specified. However, there
might still be misspecifications. It could be that the misspecification is
smaller than can be detected by the available data. Or, it could be that
omitted confounders might affect the structure and adequacy of the included
model component. This is an important consideration for both testing and
interpreting the model. For example, in our application study, we interpret
the reciprocity effect by referring to the impact of the work-time may have on
event patterns. We did not explicitly include the time-of-day in the model
formulation (as it is part of the baseline hazard), resulting in the fact that
it now "appears” in the other effects that are included in the model.

##  8 List of Abbreviations

  1. 1.

REM: relational event model;

  2. 2.

GOF: goodness of fit;

  3. 3.

FLE: fixed linear effect;

  4. 4.

NLE: non-linear effect;

  5. 5.

TVE: time-varying effect;

  6. 6.

RE: random effect;

  7. 7.

PL: partial likelihood;

  8. 8.

NCC: nested case control;

  9. 9.

PMLE: penalized maximum likelihood estimator;

  10. 10.

DGP: data generation process;

  11. 11.

CS: correctly specified;

  12. 12.

MS: misspecified;

##  9 Funding

This work was supported by funding from the Swiss National Science Foundation
(grant 192549).

##  10 Authors’ contributions

MB and EW designed the study and developed the methodology. MB performed the
computational analyses in the manuscript. MB wrote the manuscript with
contributions from EW. Both authors read and approved the manuscript.

##  11 Competing Interests

The authors declare that they have no competing interests.

##  12 Acknowledgments

Not applicable.

## References

  *   * Uzaheta et al. [2023] Uzaheta, A., Amati, V., Stadtfeld, C.: Random effects in dynamic network actor models. Network Science, 1–18 (2023)
  * Juozaitienė et al. [2023] Juozaitienė, R., Seebens, H., Latombe, G., Essl, F., Wit, E.C.: Analysing ecological dynamics with relational event models: the case of biological invasions (2023)
  * Vu et al. [2017] Vu, D., Lomi, A., Mascia, D., Pallotti, F.: Relational event models for longitudinal network data with an application to interhospital patient transfers 36(14), 2265–2287 (2017) https://doi.org/10.1002/sim.7247 . Accessed 2022-09-22
  * Filippi-Mazzola et al. [2023] Filippi-Mazzola, E., Bianchi, F., Wit, E.C.: Drivers of the decrease of patent similarities from 1976 to 2021. PLOS ONE 18(3), 1–13 (2023) https://doi.org/10.1371/journal.pone.0283247
  * Butts [2008] Butts, C.T.: 4\. a relational event framework for social action 38(1), 155–200 (2008) https://doi.org/10.1111/j.1467-9531.2008.00203.x . Accessed 2022-09-22
  * Perry and Wolfe [2013] Perry, P.O., Wolfe, P.J.: Point process modelling for directed interaction networks 75(5), 821–849 (2013) https://doi.org/10.1111/rssb.12013 . Accessed 2022-09-22
  * Boschi et al. [2023] Boschi, M., Juozaitienė, R., Wit, E.-J.C.: Smooth alien species invasion model with random and time-varying effects. arXiv preprint arXiv:2304.00654 (2023)
  * Filippi-Mazzola and Wit [2024] Filippi-Mazzola, E., Wit, E.C.: A stochastic gradient relational event additive model for modelling us patent citations from 1976 to 2022. Journal of the Royal Statistical Society Series C: Applied Statistics, 023 (2024)
  * Lerner and Lomi [2022] Lerner, J., Lomi, A.: A dynamic model for the mutual constitution of individuals and events. Journal of Complex Networks 10(2), 004 (2022)
  * Bianchi et al. [2023] Bianchi, F., Filippi-Mazzola, E., Lomi, A., Wit, E.C.: Relational Event Modeling (2023)
  * Brandenberger [2019] Brandenberger, L.: Predicting network events to assess goodness of fit of relational event models 27(4), 556–571 (2019) https://doi.org/10.1017/pan.2019.10 . Accessed 2022-09-23
  * Amati et al. [2024] Amati, V., Lomi, A., Snijders, T.A.: A goodness of fit framework for relational event models. Journal of the Royal Statistical Society Series A: Statistics in Society, 016 (2024)
  * Lin et al. [1993] Lin, D.Y., Wei, L.-J., Ying, Z.: Checking the cox model with cumulative sums of martingale-based residuals. Biometrika 80(3), 557–572 (1993)
  * Marzec and Marzec [1997a] Marzec, L., Marzec, P.: Generalized martingale-residual processes for goodness-of-fit inference in cox’s type regression models. The Annals of Statistics 25(2), 683–714 (1997)
  * Marzec and Marzec [1997b] Marzec, L., Marzec, P.: On fitting cox’s regression model with time-dependent coefficients. Biometrika 84(4), 901–908 (1997)
  * Marzec and Marzec [1998] Marzec, L., Marzec, P.: Testing based on sampled data for proportional hazards model. Statistics & probability letters 37(3), 303–313 (1998)
  * Borgan and Zhang [2015] Borgan, Ø., Zhang, Y.: Using cumulative sums of martingale residuals for model checking in nested case-control studies. Biometrics 71(3), 696–703 (2015)
  * Michalski et al. [2011] Michalski, R., Palus, S., Kazienko, P.: Matching organizational structure and social network extracted from email communication. In: Business Information Systems: 14th International Conference, BIS 2011, Poznań, Poland, June 15-17, 2011. Proceedings 14, pp. 197–206 (2011). Springer
  * [19] GOF GitHub Repository. https://github.com/martinaboschi/GOF.git
  * Bianchi and Lomi [2024] Bianchi, F., Lomi, A.: From ties to events in the analysis of interorganizational exchange relations (2024)
  * Juozaitienė and Wit [2022] Juozaitienė, R., Wit, E.C.: Non-parametric estimation of reciprocity and triadic effects in relational event networks. Social Networks 68, 296–305 (2022)
  * Hastie and Tibshirani [1990] Hastie, T.J., Tibshirani, R.J.: Generalized Additive Models. Chapman & Hall/CRC Monographs on Statistics & Applied Probability, (1990)
  * Wood [2017] Wood, S.N.: Generalized additive models: An introduction with r (2017)
  * Juozaitienė and Wit [2024] Juozaitienė, R., Wit, E.C.: Nodal heterogeneity can induce ghost triadic effects in relational event models. Psychometrika 89(1), 151–171 (2024)
  * Borgan et al. [1995] Borgan, O., Goldstein, L., Langholz, B.: Methods for the analysis of sampled cohort data in the cox proportional hazards model. The Annals of Statistics, 1749–1778 (1995)
  * Vu et al. [2015] Vu, D., Pattison, P., Robins, G.: Relational event models for social learning in moocs. Social Networks 43, 121–135 (2015) https://doi.org/10.1016/j.socnet.2015.05.001
  * Lerner and Lomi [2020] Lerner, J., Lomi, A.: Reliability of relational event model estimates under sampling: How to fit a relational event model to 360 million dyadic events. Network science 8(1), 97–135 (2020)
  * Pinheiro and Bates [2006] Pinheiro, J., Bates, D.: Mixed-effects Models in S and S-PLUS, (2006)
  * Wood [2003] Wood, S.N.: Thin-plate regression splines. Journal of the Royal Statistical Society (B) 65(1), 95–114 (2003)
  * Oelker [2015] Oelker, M.-R.: Penalized regression for discrete structures. PhD thesis, lmu (2015)
  * Zeileis and Hornik [2007] Zeileis, A., Hornik, K.: Generalized m-fluctuation tests for parameter instability. Statistica Neerlandica 61(4), 488–508 (2007)
  * Hjort and Koning [2002] Hjort, N.L., Koning, A.: Tests for constancy of model parameters over time. Journal of Nonparametric Statistics 14(1-2), 113–132 (2002)
  * Liu and Xie [2020] Liu, Y., Xie, J.: Cauchy combination test: A powerful test with analytic p-value calculation under arbitrary dependency structures. Journal of the American Statistical Association 115(529), 393–402 (2020)
  * Michalski et al. [2014] Michalski, R., Kajdanowicz, T., Bródka, P., Kazienko, P.: Seed selection for spread of influence in social networks: Temporal vs. static approach. New Generation Computing 32, 213–235 (2014)
  * Wit et al. [2012] Wit, E., Heuvel, E.v.d., Romeijn, J.-W.: ‘all models are wrong…’: an introduction to model uncertainty. Statistica Neerlandica 66(3), 217–236 (2012)
  * Schecter et al. [2018] Schecter, A., Pilny, A., Leung, A., Poole, M.S., Contractor, N.: Step by step: Capturing the dynamics of work team process through relational event sequences. Journal of Organizational Behavior 39(9), 1163–1181 (2018)