A Bayesian actor-oriented multilevel relational event model with hypothesis testing procedures

Relational event network data are becoming increasingly available. Consequently, statistical models for such data have also surfaced. These models mainly focus on the analysis of single networks, while in many applications, multiple independent event sequences are observed, which are likely to display similar social interaction dynamics. Furthermore, statistical methods for testing hypotheses about social interaction behavior are underdeveloped. Therefore, the contribution of the current paper is twofold. First, we present a multilevel extension of the dynamic actor-oriented model, which allows researchers to model sender and receiver processes separately. The multilevel formulation enables principled probabilistic borrowing of information across networks to accurately estimate drivers of social dynamics. Second, a flexible methodology is proposed to test hypotheses about common and heterogeneous social interaction drivers across relational event sequences. Social interaction data between children and teachers in classrooms are used to showcase the methodology.


Introduction
Current technological advancements, combined with constant development of new communication applications, has originated huge amounts of data on social interactions (Eagle & Pentland 2003).This has enabled researchers to build and test increasingly rich and complex models of human interaction using time-stamped data.Social network analysis is among the fields that have benefited the most from having access to rich temporal interaction data.In this paper, we specifically focus on the Relational Event network model (Butts & Marcum 2017;Stadtfeld & Block 2017;Mulder & Leenders 2019) A relational event network is comprised of multiple interactions among a finite set of actors.Observation units are called "relational events", which are defined as discrete instances of interactions among social entities along a timescale (Butts & Marcum 2017).In a directed network, events contain a clear indication of sender and receiver.Thus, every observation unit displays information on which actor was the sender, which actor was the receiver, and at which point in time the interaction occurred.The stream of events in a given network is often referred to as a relational event sequence.
One of the traditional approaches to social network analysis is based on Markovian random graphs theory and involves the aggregation of events into a graph (Van Der Hofstad 2009).Transitions on the graph structure are then modeled via sufficient statistics, which are functions of the observed network (Frank & Strauss 1986;Hanneke, Fu, & Xing 2010;Lusher, Koskinen, & Robins 2013).A different approach was taken by Butts (2008), who employed survival analysis concepts and introduced the relational event model which profits from the temporal structure of relational events.In this framework, the main goal is to model the rates of communication between sender and receivers via a log-linear function, without requiring the data to be aggregated.This is done by employing sufficient statistics that capture important social patterns and actor-specific covariates.
Since then, the relational event model has gained increasing popularity and received multiple extensions.For instance, Vu et al. (2011) introduced a model with time-varying parameters, based on the additive Aalen model (Aalen 1989), and developed an algorithm for online inference in social networks.Perry & Wolfe (2013) used a partial likelihood approach to modeling who the receiver will be for a given sender.Vu et al. (2015) implemented case-control sampling to decrease the number of computations in estimating relational event models.Later, the control-case sampling approach was further explored by Lerner & Lomi (2020) to estimate relational event models in large networks.Stadtfeld, Hollway, & Block (2017a) and Stadtfeld & Block (2017) introduced a two-step model, where the first step consists of modeling the activity rate of the sender and the second one features the choice of the receiver conditional on the sender.Mulder & Leenders (2019) proposed a way to analyze the temporal evolution of effects in the social network, by estimating effects in different subsets of the data determined by overlapping intervals, which they called moving windows, this approach was further developed by Meijerink et al. (2022).
These contributions have mainly focused on the statistical analysis of a single relational event sequence (with a notable exception of DuBois et al. (2013)) while in practice we also often observe multiple independent relational event sequences which show similar but not identical social interaction behavior.For example, Blonder & Dornhaus (2011) analyses information flow in multiple relational event networks of ant colonies, DuBois et al. (2013) studies several relational event sequences of high-school students interactions and Kiti et al. (2016) examines social contact in numerous networks of Kenyan households.The gold standard for analyzing such clustered (or hierarchical) data is using a multilevel or mixed effects approach.In a relational event modeling approach, this implies that a probability distribution is specified for the network effects (or drivers) of social interaction behavior (such as inertia, reciprocity, transitivity, or group effects) across clusters.A multilevel approach will result in 'pooled' estimates of the clusterspecific network effects where information about social interaction behavior in other clusters is borrowed to improve the estimation (Gelman & Hill 2006).We achieve this by modeling the relational event sequences in a Bayesian hierarchical framework via Markov chain Monte Carlo (MCMC) methods.Therefore, our first contribution is a multilevel extension for the dynamic actor-oriented model (Stadtfeld et al. 2017a) for independent relational event sequences.In this framework, a relational event between a sender and a receiver in a network is modeled by separately modeling (i) when an actor decides to initiate an event as a sender given the past history, and (ii) modeling the receiver that is chosen given the sender and the past event history between the actors.Thereby, this approach differs from the dyadic approach of Butts (2008), and the hierarchical extension by DuBois et al. (2013), where the time, sender, and receiver are jointly modeled using the same set of parameters (Stadtfeld et al. 2017b).The actor-oriented approach is therefore more flexible by separately modeling the behavior of actors to become a sender and the behavior of senders to choose a receiver in contrast to the dyadic approach.A second important difference is that we explicitly model the full unstructured covariance matrix of the network effects across sequences.Thereby the dependency between interaction behavior across sequences, such as inertia and reciprocity (which is generally nonzero), is included in the analysis resulting in improved estimation using the concept of pooling.A final important difference we note here between a dyadic relational event model and an actor-oriented approach is the computational burden, which is considerably larger in the case of a dyadic approach.
This can be explained by the size of the risk set (i.e., the possible events that can be observed at every step), which is equal to N (N − 1) in the case of a network consisting of N actors where N (N − 1) dyads are at risk for each event, while the risk set of the sender model in an actor-oriented model consists of N actors, and the risk set of the receiver model consists of N − 1 actors (when assuming that the sender cannot be equal to the receiver).In the case of moderate to large networks, this difference can have a huge impact on the computational costs when using a dyadic approach.
Our second contribution is an extensive set of statistical tests using the Bayes factor (Jeffreys 1961;Kass & Raftery 1995) for evaluating hypotheses about social interaction behavior in the case of independent relational event sequences.The first test can be used to assess whether specific interaction behavior, such as inertia (which quantifies the degree of habitual behavior between actors to keep sending messages of one actor towards another actor), is equal or differs across sequences (i.e., whether it should be modeled as a fixed effect or as a random effect).In the case of different interaction behavior across sequences, a second test is proposed to assess the degree of heterogeneity (or variability) of network effects on social interaction behavior, e.g., to assess whether inertia behavior varies more across clusters than reciprocity, or whether the impact of a group effect (e.g., whether the sender is the teacher in a class room or not) on social interaction dynamics or more or less heterogeneous than another group effect (e.g., whether the sender is male or female).Third, a test is proposed for evaluating hypotheses with equality and order constraints on the average relative importance across clusters (i.e., on the fixed effects) of drivers on social interaction dynamics.
The methodology is flexible to also allow for testing the absolute values of network effects, which is useful when the exact sign is unknown before observing the data.These tests build on and further extends recent developments on Bayes factor hypothesis testing (Mulder et al. 2021).
Our methods have been implemented in R (R Core Team 2017), interfacing with Stan through the rstan package (Stan Development Team 2018).We provide the full code to facilitate the fitting of other relational event models (see link suppressed for peer review).
The remainder of this paper is divided as follows: in Section 2 we describe the relational events framework and present the hierarchical actor-oriented relational event model.This Section also introduces the Bayesian specification of the model by presenting prior distributions and briefly discussing a reparametrization that allows more efficient sampling; we detail the hypothesis testing methods in Section 3 and conduct the empirical application in Section 4. The paper ends with a discussion in Section 5, where we go through some limitations of the model and possible routes for future research.

The relational event framework
The relational event model (REM) is used to model the rate of interactions in dynamic social networks (Butts 2008).
In this framework, events happen among actors at particular points in time, being represented by tuples in the form et = (s, r, t), where s is the sender, r is the receiver, and t is the time point of the interaction.Thus et is called a relational event.Therefore, we have K clusters (in this paper we will use cluster, groups, relational event sequence and network interchangeably), each of which with N1, . . ., NK social actors allowed to be senders or receivers at any given time.At time t, in cluster k, sender s ∈ {1, 2, . . ., N k } interacts with receiver r ∈ {1, 2, . . ., N k | r = s}, forming a dyad (s, r) ∈ R k (t), where R k (t) = (s, r) : s, r ∈ {1, 2, . . ., N k }, s = r is called the risk set, which comprises all possible dyads (s, r) at a particular point in time.
The main idea consists of assuming that, for each cluster k, for k = 1, . . ., K, we are able to observe an ordered sequence of M k dyadic events among N k individuals on the time window [0, τ k ) ∈ IR + .Thus a relational event sequence for cluster k is formally defined as where tm is the time at which the m th event occurred, in cluster k.Then, following Butts (2008), who borrows concepts from survival analysis, the relational event history E k is modeled as a stochastic process, where the rate of events from sender s to receiver r , with (s, r) ∈ R k (t), is given by the intensity λsr t|E k .This intensity function has the form of a Cox proportional hazards model (Cox 1972).Moreover, the intensity is assumed constant between subsequent events and the waiting times conditionally exponentially distributed.This amounts to the well know piecewise constant exponential model for survival data (Friedman 1982).Thus, the survival function is given by

The actor-oriented multilevel relational event model
In this paper, we focus on an alternative relational event model that was proposed by Stadtfeld et al. (2017a) and Stadtfeld & Block (2017).Their approach models the receiver given the sender, similar to Perry & Wolfe (2013).Conceptually, it builds on the same tradition as the stochastic actor-oriented model (Snijders 1996), where the evolution of the network is assumed to be a product of actors' individual behaviors as they constantly seek to maximize their own utilities.The basic framework consists of a two-step approach based on log-linear predictors.
First, the waiting time until an actor becomes active is modeled.After that, a multinomial choice model (McFadden 1973) is employed to determine the choice of the receiver by the active sender actor.These two steps are assumed to be conditionally independent given the available information about the past up to that event.
The Bayesian hierarchical approach we will develop for multiple relational event sequences (denoted as "clusters") will build on this model.At each point in time, in cluster k, sender s ∈ {1, 2, . . ., N k } starts an interaction with intensity λs t|E k .This intensity is directly proportional to the probability of a given actor s to be the next sender.

This probability is given by
where M k is the number of events in cluster k.Next, this sender chooses the receiver r ∈ {1, 2, . . ., N k |r = s}, forming the dyad (s, r) ∈ R k (t), with intensity λ r|s t|s, E k .The receiver intensity represents the rate at which actor s chooses actor r to form a dyad, which is proportional to the probability of observing dyad (s, r) ∈ R k (t) as the next one in the sequence.This probability is given by P rm = r|sm = s, E k = λ r|s t|s, E k / h∈k λ h|s t|s, E k , ∀ h ∈ {1, 2, . . ., N k | h = s} and m ∈ {1, 2, . . ., M k }.Then, these intensities for cluster k, for the sender and receiver steps of this model, will be given by the following log-linear functions: where φ and ψ are a vectors of fixed-effect parameters, γ k and β k are a vectors of random-effect parameters for cluster k, with k = 1, . . ., K. The vectors of statistics zs(t) and xs(t) are associated with actor s ∈ {1, 2, . . ., N k }, whereas zsr(t) and xsr(t) are vectors of statistics associated with dyad (s, r) ∈ R k (t).
Assuming this log-linear form will allow us to conduct inference at the actor level, unveiling effects that make actors more (or less) prone to start interactions or more (or less) likely to be chosen as the next receiver.Also, this idea helps us limit the size of the set of possible dyads that need to be analyzed at each point in time: if actor s is the one starting an interaction at time t, then all dyads where s is not the sender become impossible to happen.
Thus, for the actor-oriented model, the risk set, R k (t), will have size N k − 1 (given the sender), whereas for a dyadic model, such as in DuBois et al. (2013), the risk set has size (N k − 1)N k , which can easily become massive for medium to large networks.
The likelihood for the actor-oriented model is given by where θ is the vector containing all parameters and Z and X are matrices with fixed and random effects covariates, respectively.The time of the last observed event in cluster k is denoted by tM k and τ k the end of the observation period.
In most empirical applications, it is assumed that tM k = τ k .This way the last part of the likelihood is equal to 1, since due to the piecewise constant exponential assumption we have Ss m tm − tm−1|E k = exp − (tm − tm−1)λs tm|E k .In this setting, at time tm, P s = sm|E k = λs h∈k λ h , ∀ s ∈ k is the probability of sender s being active.The probability of actor r being the receiver given that s is the sender is given by The likelihood in equation ( 2) is a product of the two pieces.The first one, representing the sender model, is a piece-wise constant exponential likelihood and the second one, representing the receiver model, is a multinomial likelihood In the literature, it has been shown that both of these models are special cases of the Poisson model.Holford (1980) and Laird & Olivier (1981) are examples of cases where the Poisson representation of the piece-wise exponential model is discussed.Baker (1994) extensively discusses the multinomial-Poisson transformation and how their likelihoods yield identical estimates.For proofs, see Appendix A. The Poisson regression model has been extensively studied, and it is well understood in the statistical literature (Frome 1983;Consul & Famoye 1992;Hayat & Higgins 2014).
The second level of the hierarchical actor-oriented relational event model specifies the multivariate distributions of network effects across the K sequences.We follow the standard approach in multilevel modeling by assuming multivariate normal distributions for β k under the sender model and for γ k under the receiver model, i.e., for k = 1, . . ., K. The mean vectors quantify the overall, global effect of the network effects and the unstructured covariance matrices quantify the variability of the effects across clusters and the dependency structure between the effects across clusters.Thus, when estimating these distributions when fitting the multilevel model using independent relational event sequences, the estimated means and (co)variances of the second level are used to improve the estimates of the specific parameters in the separate sequences (especially in the case of short sequences).

Prior specification
In a Bayesian approach, prior distributions have to be specified which reflect our uncertainty about the model parameters before observing the data.Throughout this paper we shall work with vague, noninformative priors (which are completely dominated by the data).
It has been shown that a half-Cauchy prior for the random effects standard deviation results in desirable estimates in the case of multilevel data with few clusters (Gelman & Hill 2006).Furthermore, by setting very large prior scale parameters, we can obtain approximately flat priors for the standard deviations.Further note that a half-Cauchy prior for the standard deviation corresponds to a F distribution on the variance (Mulder & Pericchi 2018), which is a common distribution for modeling variance components.
The LKJ prior is defined as LKJCorr(Ω|η) ∝ det(Ω) η−1 , with η ∈ IR + .This distribution allows us to sample uniformly from the space of positive definite correlation matrices and has a behavior similar to the beta distribution (Wang, Wu, & Chu 2018;Lewandowski, Kurowicka, & Joe 2009).For example, when η = 1 it has a uniform behavior, when η < 1 it favors stronger correlation, whereas when η > 1 it favors weaker correlation.For the random effects covariance matrix under the sender Σγ, the same prior is specified.
Finally, multivariate normal priors are specified for the fixed effects and random effects, i.e., For the prior covariance matrices, diagonal matrices are specified with very large variances.This builds upon existing knowledge of weakly informative prior specification for mixed effects generalized linear models, which have become the standard approach under this class of mixed effects models (Gelman et al. 2013;Gelman & Hill 2006).

Reparameterization to improve Bayesian computation
Due to the hierarchical structure of the data, the random-effects parameters β k and γ k are highly correlated with the population parameters µ, ζ and Σ.This introduces severe computational inefficiencies of the sampling process.When the data are sparse, which is a characteristic of most social network data, the geometry of the posterior distribution makes it very difficult to sample from the highest posterior density areas.Betancourt & Girolami (2015) called these issues pathologies of the hierarchical model.Therefore, to ease the burden on the sampler, we take advantage of the multivariate normal structure of the random effects and apply a non-centered linear transformation to those parameters.
Lemma 1.Let β ∼ N (µ, Σ), where β ∈ IR p .Then, with µ ∈ IR p and A being a p × p matrix, such that AA = Σ, one can write β = AZ + µ, where Z ∈ IR p and Z ∼ N (0, I), where I is a p × p identity matrix.
This transformation can be applied to both β and γ.A natural candidate for matrix A is the Cholesky factorization of the covariance matrix Σ.For details see appendix B.
This reparameterization is more efficient for two reasons.First, it reduces the dependency between the random effects parameters and the population parameters by sampling from independent standard normal distributions.This simplifies the geometry of the posterior and avoids inverting Σ at every evaluation of the multivariate normal density (Carpenter et al. 2017).Therefore, we can safely and efficiently transform the random-effects parameters without causing any change in the prior specification of population parameters.
The model has been implemented in Stan, a probabilistic programming language that employs Hamiltonial Monte Carlo (HMC) algorithms to sample from posterior distributions.The advantage of HMC methods is that they avoid the random walk behavior and the sensitivity to posterior correlations that plague many Bayesian applications (Hoffman & Gelman 2014), including hierarchical models.This particularities allow HMC to generally converge to high dimensional distributions much faster than Metropolis-Hastings or Gibbs sampler methods (Hoffman & Gelman 2014;Betancourt & Girolami 2015).
3 Bayesian hypothesis testing under the actor-oriented multilevel relational event model In multilevel analysis, researchers are usually interested in testing which theories receive the most support from the observed data, so that inferences about the population can be conducted.This kind of analysis is carried out through a process called hypothesis testing.From a Bayesian perspective those tests are usually performed by the computation of Bayes factors (BF) (Jeffreys 1961).The BF is given by the ratio of marginal likelihoods under the parameter space of competing hypotheses, which quantifies the probability of observing the data under one hypothesis relative to another hypotheses, and thereby providing a quantification of the relative evidence in the data between the hypotheses.Let E be the observed data, θ a vector of parameters in the space Θ, and H0 ∈ Θ0 be a hypothesis that will be tested against H1 ∈ Θ1, then BF01 is expressed as where m(E|Hi) = θ i ∈Θ i p(E|θi)p(θi)dθi, for i = 0, 1, with p(E|θi) being the likelihood and p(θi) the prior.Also, Θ0 ∩ Θ1 = ∅, with both Θ0 and Θ1 being subsets of Θ. Kass & Raftery (1995) provide a rule-of-thumb for Bayes Factors interpretation.In their setting, the evidence provided by the BF01 in favor H0 can be seen as "insufficient" if 1 < BF01 < 3, "positive" if BF01 > 3, "strong" if BF01 > 20, and "very strong" if BF01 > 150.These are rough guidelines to aid the interpretation of Bayes factors and should not be used as strict cut-off values.
Furthermore, when prior probabilities of the hypotheses have been formulated before observing the data, the prior odds can be updated with the Bayes factor to obtain the posterior odds of the hypotheses, i.e., P (H0|E) P (H1|E) = BF01 × P (H0) P (H1) .
For example, when both hypotheses are equally likely a priori, i.e., P (H0) = P (H1) = 1 2 , the posterior probabilities of the hypotheses are given by P (H0|E) = BF 01 BF 01 +1 and P (H1|E) = 1 BF 01 +1 .These posterior probabilities quantify the probabilities that each of these hypotheses are true after observing the data (when assuming that one of these hypotheses is true).Because Bayes factors are consistent under mild conditions, the evidence for the true hypothesis goes to infinity and the posterior probability of the true hypothesis goes to 1 as the sample goes to infinity.
Below we introduce Bayes factors for three types of hypothesis testing problems for the multilevel relational event modeling literature.Test I is useful for testing the homogeneity of network effects across independent sequences.
Test II is useful for testing the degree of heterogeneity of network effects across sequences using order constraints on random effects variances.Test III is useful to test equality and order constraints on fixed effects.Prior specification and numerical computation, which are important aspects when computing Bayes factors and posterior probabilities of the hypotheses, are discussed for each test separately.

Test I: Testing for homogeneous social interaction behavior across sequences
When building multilevel relational event models, a central question is whether a driver of specific social interaction behavior (as quantified via the network effects) is constant over the sequences or whether it varies across the sequences.
Testing this is important from a statistical point of view (i.e., to keep the model as parsimonious as possible, and thus to avoid an enormous overparameterization using different effects across all K sequences) but also from a substantive point of view (i.e., to understand which drivers of social interaction behavior in relational event data are constant across sequences and which (highly) differ).Thereby, testing this contributes to a better understanding of the heterogeneity of social interaction behavior.
The hypothesis test for the p-th random network effect in the receiver model can be formulated as H0 : σ β,p = 0 versus H0 : σ β,p > 0, or equivalently as, H0 : βp,1 = . . .= βp,K versus H1 : not H0.The second formulation of the hypothesis is simpler as we avoid the need of testing whether the variance is equal to the boundary of 0, but instead we only need to test whether the random effects are equal across the K clusters.Using the fact that the distribution of the random effects effectively serves as a prior for the random effects on the first level, we can simply compute a Savage-Dickey density ratio using the estimates of the random effects distribution from the data, i.e., . ., K, are the differences between the random effects.
To compute the posterior and prior densities at the null value, we use the fact that our actor-oriented multilevel relational event model can be written as two independent mixed effects Poisson regression models having specific forms (Appendix A).Thus the two independent components of the model belong to the family of generalized linear mixed effects models, which is well understood in the statistical literature (e.g., Gelman & Hill 2006;Bürkner 2017).
Moreover, the model results in unimodal posteriors when using weakly informative priors, as done in this paper.
We illustrate the accuracy of this Gaussian approximation later in this paper.The Gaussian approximation of the posterior is then used to compute the numerator (then having an analytic expression).Moreover, the prior distribution, which is needed to compute the denominator, follows a multivariate Gaussian distribution by definition because it is the second level of our multilevel model, β k ∼ N (µ β , Σ β ), and thus the joint prior distribution of the contrasts (ξ) also follows a multivariate Gaussian distribution.
One interesting aspect of this Bayes factor is that the prior is fully determined from the data, similar to empirical Bayesian estimation of hierarchical models.This property is especially useful here as Bayes factors are known to be sensitive to the choice of the prior.Note that empirical Bayesian approaches to obtain Bayes factors have been proposed for testing regression coefficients using g priors (e.g., see Liang et al. 2008, and the references therein), but (to our knowledge) not for testing variance components as we do here.The advantage of this approach is that no external prior information is required and no other ad-hoc choices for prior distributions are needed.The outcome of the test can be used to quantify the relative evidence in the data of whether a hierarchical structure for the network effects is applicable or not given the observed data.

Test II: Testing the degree of heterogeneity of network effects across sequences
After establishing which network effects are heterogeneous across sequences (using the test from the previous subsection), it is useful to investigate the relative degree of heterogeneity of network effects across sequences, again with the goal to better understand the heterogeneity of social interaction behavior across independent relational event sequences.This comes down to testing whether a specific random effect varies more across sequences than another random effect, or, equivalently, whether one random effect variance is larger than another.When generalizing this further, we would test specific orderings of random effect variances (see also Böing-Messing & Mulder 2020, who consider testing order constraints on variances in a nonhierarchical setting).
Two different tests are proposed for this purpose.First a confirmatory test is proposed of whether an anticipated ordering regarding the degree of heterogeneity of network effects across sequences is present: Second, an exploratory test is proposed to determine which ordering out of all P !possible orderings receives most evidence from the data, i.e., Following the existing Bayesian literature on order constrained hypothesis testing (e.g.Klugkist et al. 2005), we specify the prior under an order-constrained hypothesis as a truncated version of an unconstrained prior truncated in the order-constrained subspace.The Bayes factor of an order-constrained hypothesis against the unconstrained hypothesis, Hu, can then be expressed as the ratio of the posterior and prior probability that the constraints hold (Appendix C): The Bayes factors between constrained hypotheses of interest can then be obtained using the transitive property of the Bayes factors, e.g., B12 = B1u/B2u.
The Bayes factor will not be sensitive to the prior as long as very vague identical priors are specified for the variances.In this case, the prior probability in the denominator will be equal to (P !) −1 , and the posterior probability will be fully determined by the information in the data (similar as in Bayesian estimation).Bartlett's paradox is thus not an issue in Bayesian order or one-sided hypothesis testing (Jeffreys 1961;Klugkist & Hoijtink 2007;Liang et al. 2008;Mulder 2014a).
To compute the posterior probability that the order constraints hold under the unconstrained model, we can simply fit the unconstrained model (using independent vague priors) and compute the proportion of draws satisfying the order constraints.

Test III: Testing common and average network effects over all sequences
In most applications, hypothesis tests are formulated on the common network parameters across clusters (i.e., the fixed effects) and on the average of the random-effects across clusters (i.e., the global random effects means).Since both parameters have the same role, hypothesis tests of these parameters are both discussed here using the same methodology.
First, a common test is of whether a network statistic has no effect, a negative effect, or a positive effect on the interaction behavior, which could be formulated as H0 : ψ1 = 0 versus H1 : ψ1 < 0 versus H2 : ψ1 > 0 for the first fixed effect of the receiver model, for instance.Second, hypotheses are often formulated on combinations of parameters using equality and/or order constraints on multiple parameters of interest based on existing theories or scientific expectations (Hoijtink et al. 2019), e.g., where hypothesis H1 assumes that the first effect (ψ1) is equal to the second effect (ψ2) and larger than the third effect (ψ3), and hypothesis H2 assumes the complement is true.Third, when testing the effects of categorical (dummy) variables on the event rate, scientific expectations may be formulated which categorical variable has the largest impact on the event rate.A challenging testing problem is when the interest is in assessing which categorical variable has the largest impact when it is not known which category of each categorical variable results in the largest event rate.
For example, it may be of interest whether the dichotomous gender variable (0 or 1) has a larger, equal, or smaller impact on the event rate than a dichotomous race variable (0 or 1), but it is not of particular interest which category results in the largest rate.This could be translated to constrained hypotheses on the absolute values of the effects of these two categorical variables according to Finally, it is possible to test parameters between the sender model and the receiver model.This allows researchers to assess whether variable has a larger/equal/smaller to predict the next sender than to predict the next receiver.
The procedure for Test III builds on default Bayes factor methodology where the information in the data is split between a minimal subset, which is used for default prior specification (in combination with a noninformative prior) and a maximal subset which is used for hypothesis testing (O'Hagan 1995; Pérez & Berger 2002, among others).
When using these methods, the resulting Bayes factors do not depend on the undefined normalizing constants of the improper priors (O'Hagan 1995), which practically means that arbitrarily vague priors can be used, as we do in this paper.These Bayes factors can be computed in an automatic fashion and manual prior specification can be avoided 1 .Here we use fractional Bayes factors where a (minimal) fraction of the data, denoted by b, is used to construct a (default) fractional prior, given by and the remaining fraction is used for hypothesis testing.It has been shown that fractional Bayes factors follow many important properties which are not shared by all default Bayes factors, such as large sample consistency, satisfying the likelihood principle, invariance to transformations of the data (O'Hagan 1997), and information consistency (Mulder 2014b).To properly incorporate the model complexity of order-constrained hypotheses, we apply a prior adjustment of the fractional prior (Mulder 2014b), and Gaussian approximations are applied to the posterior and fractional prior to simply the computation (Kim & Ibrahim, 2000;Gu, et al., 2017).The accuracy of the Gaussian approximations will be discussed in the next section.The approximated fractional Bayes factor of a constrained hypothesis (e.g., with any set of equality and/or constraints on certain parameters) against an unconstrained model is then computed as the integral over the unconstrained posterior over the constrained subspace divided by the integral over the unconstrained fractional prior over the constrained subspace as a Savage-Dickey density ratio (Mulder 2014b).Existing functions in R using the 'mvtnorm' package (Genz et al., 2021) can be used to compute the Bayes factors (Mulder et al. 2021).
The choice of the fraction 'b' is based on the recommendation by Mulder & Fox (2019), which implies selecting a minimal sample that is based on the ratio between the total number of parameters (excluding group specific effects) and the total number of observations.Since sender and receiver models are two different processes assumed to be conditionally independent given the past, we have to define two difference fractions, one for each model (see also Hoijtink, Gu, & Mulder 2019).In the sender model, we have P random effects and Q fixed effects, hence the fraction is b Snd = (P (P +1)/2)+(P +Q) and in the receiver model we have V random effects and U fixed effects, resulting in the 4 Estimating and testing classroom dynamics using the multilevel relational event model The data that are considered here were collected by McFarland (2001)

Model specification
We fit the hierarchical actor-oriented relational event model (described in Section 2) to the data.Initially, all effects are considered random, so the need for a hierarchical structure can be evaluated using the Bayes factor presented in Section 3.This first step will allow us to determine whether some effects can be treated as fixed and the model can hence be simplified.
We set the prior parameters to σ 2 φ = σ 2 µ β = σ 2 ψ = σ 2 ζγ = στ = 10, making the priors for the fixed and random effects relatively vague, and η = 2, slightly favoring smaller correlations.For our illustration of the proposed methods, we will include the following covariates into the model: Teacher: A dummy variable that indicates whether the actor is the teacher (one if the actor is the teacher and zero otherwise).
Gender: A dummy variable that indicates the gender of the actor (one if the actor is male and zero otherwise).
Race: A dummy variable that indicates the race of the actor.McFarland (2001) notes that 50 percent of the Magnet High population is Caucasian.Hence, the variable is one if the actor is Caucasian and zero otherwise.
Inertia: Inertia captures the persistence of the communication, where past interaction is likely to be repeated (Leenders, Contractor, & DeChurch 2016).It is computed as the accumulated volume of past communication from a specific sender to a specific receiver.Because this statistic is dyadic, it will be included only in the receiver model.
Participation Shifts: These statistics are used to reflect expectations of adherence of communication norms in small groups (Butts 2008).Gibson (2005) describes the framework for several types of participation shifts.In this application two distinct types of participation shifts will be included.The first belongs to the group of "turn-receiving" and is represented by the event pattern ABBA: an interaction from person A to person B is immediately followed by an interaction from B to A. The second is ABAB, which can be considered as a special case of "turn-continuing": an interaction from person A to B is immediately followed by another interaction event from A to B.
The ABBA and ABAB participation shifts are concerned with dyads, therefore they will be included in the receiver model only.However, we do include adapted versions in the sender model.Here, the statistics become ABB (after A has spoken to B, the next event is B starting the conversation) and ABA (after A has spoken to B, A speaks again).
Activity: This set of covariates captures the effect of actor activity as a sender or as a receiver (Vu et al. 2017).
The first is the outgoingness of a person, defined as the number of events sent by one actor up until a specific point in time.This captures the tendency of the person to start conversations (or just to talk).The second is the popularity of a person, given by the number of events received by one actor up until a specific point in time.This captures the popularity of an individual as a receiver of the conversation.
Relational event models are vulnerable to process explosion, which happens when λ(t) → ∞ (often due to a feedback loop that may be caused by using statistics that are computed as cumulative sums (Aalen, Borgan, & Gjessing 2008).This is particularly a problem for inertia and activity statistics.One way to alleviate this problem is via z-score standardization at every time point, which is defined as z(t) = (x(t) − x(t))/Sx(t), where x(t) is the value of the statistic, x(t) is the sample mean, and Sx(t) is the sample standard deviation at time t.We use standardized statistics in this application.

Exploring social interaction dynamics as class time progresses
By analyzing the relational event sequences as the time during class progresses, we can explore how statistical certainty increases over time using the proposed multilevel relational event model, and how interaction dynamics between children and the teacher evolves as class time progresses.To study this, we analyze the sequences with increasing batches of 20% of the total number of events in each sequence (i.e.20%, 40%, 60%, 80%, and 100%).

Testing for homogeneous social interaction behavior across school classes
The first step is to fit the model with all effects considered random.We ran 2000 MCMC iterations and discarded the first 1000 samples as burn-in.We used R as a convergence diagnostic measure.For all models, the metric values were fairly close to 1 all, falling all below the 1.05 value recommended by Carpenter et al. (2017).As an example, Figure 1 shows several trace plots for the model run with 100% of the events.The chains show convergence.The number of posterior draws is not high because the Hamiltonian dynamics in Stan's algorithm allow for faster convergence than standard MCMC methods (Hoffman & Gelman 2014).To check this we investigated the bulk effective sample size and the tail effective sample size, which fell between 800 and 1100 for all parameters (Carpenter et al. (2017) recommends these values to be both above 100 for the sample to be considered reliable).This confirms the computational efficiency of the sampler.For other sample size recommendations, see Hecht & Zitzmann (2021).Since a more parsimonious model is preferred, we consider an effect to vary across classrooms via a random effect when the posterior probability for the hypothesis that an effect is fixed, H0, is less than 0.25, corresponding to BF01 < 1/3 (which can be interpreted as positive evidence that the effect varies across clusters; Kass & Raftery, 1995).
Figure 2 shows the evolution of the posterior probabilities of H0 as we increase the number of events.In the beginning, when the number of events is small (varying from 17 to 125 events across school classes), the uncertainty in the posterior is relatively large leading to considerable overlap in the posterior distributions of the random effects.
Therefore, the posterior probabilities for a fixed effect tend to start out large and eventually decrease to zero when the full sequences are considered for most effects, providing strong evidence for their heterogeneous nature.This is the case for most effects that are clearly heterogeneous if the sample sizes are large enough.On the other hand, for some covariates the posterior probability gets larger as the sample sizes near 100% of the events.This is the case for ABA, outgoingness and popularity in the sender model and for race in the receiver model, resulting in posterior probabilities for a fixed effect of 0.768, 0.681 and 0.579, 0.777, respectively.Thus, the parameters corresponding to those effects will be fixed in our model.Therefore, we continue with a mixed-effect model, where the sender model has ABA, outgoingness, and popularity as fixed effects, and intercept, teacher , gender , race, and ABB as random effects, and a receiver model where only race is a fixed effect and all other effects are considered random.
To understand the nonmonotonic development of the posterior probabilities where the lines first go down and eventually increase, it is useful to check the corresponding estimates of the random-effects.As an example, the top left panel of Figure 3 shows random-effect estimates of the race covariate in the receiver model, where each line represents a five different classrooms (to keep the plot clear; the other classrooms showed similar patterns).As the lines show we see that the random effects become more heterogeneous until 80% of the events are included in the analysis but when all 100% of the events are included the estimates become more homogeneous which explains the evidence for a fixed effect for the race effect in the receiver model in Figure 2 when all 100% of the events are considered.
In addition, in the bottom right panel of Figure 3, which presents the random-effect estimates for outgoingness in the sender model, a similar pattern is observed.The random effects becomes heterogeneous between 40% and 60% of the events, but they rapidly converge displaying a pattern similar to race in the sender model.This behavior is reflected in the posterior probabilities, with the line going down for 40% and 60% of the events and then getting back up again.On the other hand, the estimates on the left hand side panels of Figure 3, which display popularity in the receiver model and teacher in the sender model, are consistently heterogeneous as class time progresses.This is corroborated by the posterior probabilities with both effects presenting small evidence for the null.

Evaluating the shrinkage behavior of random effects
Given that there is asymmetry in the distribution of the number of events across classrooms (varying from 86 to 628 events over classrooms), there will be a variation regarding the statistical uncertainty of the estimated random effects across classrooms.Due to the multilevel structure of the model, estimated effects with a relatively large uncertainty which also deviate much from the other estimates will then be (slightly) pulled towards the global mean of the effect over all class rooms.This statistical behavior is also known as shrinkage, or borrowing information across groups, which is an intrinsic property of multilevel models.The current sequential analysis where the sample sizes grows allow us to investigate this behavior, since with fewer events we should have less information and, therefore, less accurate classroom specific estimates.
We use the R package remstimate (Arena et al. 2022) to obtain the estimates of the classroom specific effects  First, we observe see that there is considerably more variation of the estimated effects in the separate analyses (top of each panel), which ignore the multilevel structure, in comparison to the multilevel estimates (bottom of each panel).This confirms that the proposed multilevel model shows the anticipated shrinkage behavior.Second, we see that when using subsets of only 20% complete event sequences in the classrooms we obtain more extreme estimates than when using the proposed multilevel model.Finally we see that estimates for classrooms with shorter event sequences (represented with light grey lines) are on average more extreme than estimates for classrooms for longer event sequences.This illustrates that the proposed multilevel relational event models is particularly useful when the event sequences are relatively short to avoid extreme estimates.

Evolution of classroom behavior
Another interesting aspect to look at in the sequential analysis is the potential change of the estimated of the effects as class time progresses.This will give us an insight into the social interaction dynamics in classrooms.Figure 5 shows the evolution of some mean effects in the sender and receiver models.Here we restrict ourselves to some trends of average network behavior across all classrooms.Note that the intervals are relatively wide because the limited lengths of the event sequences in a class period.
First, we see that the teacher effect in the sender model is positive and large in the sender model.This was expected given the dominant role of the teacher during lectures.The plot in Figure 5 (upper left panel) also shows a slight decrease as class time progresses.This suggests that the dominant role of the teacher slightly decreases as the lectures approach the end.This indicates that teachers may loose some of the control towards the end or that more students become engaged or activated towards the end.This is also confirmed by the ABA effect in the sender model (Figure 5, upper right panel)) which is initially positive and decrease to just below zero towards the end of the class period.Thus, it becomes less likely that the sender of the previous messages (typically the teacher because these are lectures) becomes the sender of the next event.
Another interesting trend can be observed for the global inertia effect across classrooms in the receiver model (Figure 5, lower left panel) In the beginning of the class period, inertia starts out negative which indicates that it is less likely that actors who were the receiver of relational events in the past (typically a question or remark by the teacher towards a student because lectures are considered) become a receiver of the next.This indicates that teachers aim to address different students in the beginning of the lecture with the goal to get more students engaged or activated during lectures.The effect however gradually increases as the class time progressed.As the class period nears the end, the inertia effect becomes positive which implies that it has become more likely that a student is addressed by the teacher who was also addressed by the teacher in the past.This suggests that teachers tend to focus more on the same students (possibly the more motivated students) rather than trying to involve other less active students in the discussion.This is also confirmed by the positive and increasing ABAB effect in the receiver model during class time (Figure 5, lower right panel).

Mixed-effects analysis on the full sequences
In this subsection, we carry out the other proposed tests under the mixed-effects model obtained from the homogeneity tests that were carried out in the previous section.We only restrict ourselves to the results for the full event sequences to keep the discussion of the results as concise as possible.We ran 4 chains, each with 5000 MCMC iterations with a burn-in of 2500 draws.Figure 6 shows the trace plots of a few parameters, all chains show convergence.Here we also used the R as a convergence measure, they were essentially 1 in all chains for all parameters.In addition, the bulk and tail effective sample sizes were all close to around 2000 samples.
Moreover, Figures 7 and 8 display the density estimated from each of the 4 chains and a Gaussian approximation on top.As can be seen, the Gaussian approximation seems to be acceptable in the bulk as well as in the tails of the distribution.Therefore, empirically, the approach of conducting the tests using these approximations as proxies for the posterior distributions seems reasonable.From now on, given that the results for the four chains are virtually identical, we proceed to analyze the results of a single chain.
We can see that the proposed testing procedure to determine which effects are random and which are fixed works correctly because a similar fit to the data is obtained in the more parsimonious model (where some effects are fixed) compared to the larger (i.e.all random) model.We can see this by computing the point-wise deviance residuals for both models under the full multilevel model where all effects are random across classrooms and the mixed effects model where certain effects are assumed fixed across classrooms (according to the proposed Bayes factor test).If a similar fit is obtained for both models, the point-wise deviance residuals for both models are similar meaning that none of the two models dominates the other in terms of fit. Figure 9 clearly shows that that is the case: most points lie around the grey diagonal where the values are equal.This indicates that the fit provided by both models is similar.
Details of the point-wise deviance computation are provided in appendix E. Finally, table 1 shows a comparison of the posterior estimates for both models, where the bold ones are effects treated as fixed in the mixed-effects model.
In sum, the mixed effects model results in a comparable fit but with much less parameters to estimate, resulting in a more parsimonious model and a simpler explanation of social interaction behavior across classrooms.

Testing the degree of heterogeneity of network effects across school classes
By testing order constraints on the random effects variance parameters we can assess which characteristics cause most variation in the social interaction behavior across school classes.In particular, we focus on the random effects of the teacher , gender , and race variables in the sender model (which are all dummy variables).The same test can also be applied to the variance parameters in the receiver model in a similar manner.The objective is to determine the amount of evidence regarding the level of heterogeneity among these effects.Let σ 2 teacher , σ 2 gender , and σ 2 race be the between classrooms variances of the teacher , gender , and race effects, respectively, where these parameters are extracted from the diagonal of Σγ.Given the special (dominant) role of the teacher in a classroom, and the fact that different teachers have different teaching styles, it is expected that the variance of the teacher effect is largest.We tested all six possible order hypotheses for the random effects variances: The goal is to determine which hypothesis receives most evidence from the data.Below the Bayes factors for each one of these hypotheses against one another is represented in the form of an evidence matrix: Each cell of the matrix represent the comparison of the hypothesis in the rows against the hypothesis in the columns.
So, for example, the evidence for H1 against H2 is equal to 1.5.The evidence for H2 against H1 is then 1/1.50 = 0.67.
The evidence matrix clearly shows that the evidence for hypotheses H1 and H2 is much larger than the evidence for any of the other hypotheses, which indicates that the teacher effect shows indeed most heterogeneity across the fifteen classrooms.When inspecting the estimates of the variances, this is also confirmed: σ2 teacher = 1.65, σ2 gender = 0.25, and σ2 race = 0.22.Note that the added value of the Bayes factor complementary to eyeballing the estimates is that the effect sizes and their uncertainty are combined in a principled manner to determine which effects are most heterogeneous across classrooms.Finally note that there is no clear evidence that either the gender effect or the race effect is more heterogeneous across classrooms.Thus, we can conclude that teachers have the biggest impact on the variability of social interaction behavior across classrooms.

Testing the impact of nodal characteristics on classroom dynamics
In this subsection, we test the effects of nodal characteristics how they affect classroom dynamics.In the sender model, due to the special role of the teacher to control interaction behavior during lectures, the teacher variable is expected to have the strongest impact on the rate at which an actor starts an interaction when compared to the other two personal-trait variables gender and race.No expectations are formulated about which of these latter two variables has the largest effect in the model, or about which category of these dichotomous variables results in a positive effect.This results in the following hypotheses where constraints are formulated on the absolute values of H4 0.00 0.00 0.00 1.00 .
Again, we translate these Bayes factors to posterior probabilities when assuming that each hypothesis is equally likely a priori, which results in posterior probabilities of 0.235, 0.723, 0.042, and 0.000, for H1, H2, H3, and H4, respectively.Based on these results we can safely rule out the complement hypothesis H4.Moreover, hypothesis H2 (which assumes that the teacher variable has the largest effect, and the gender effect and the race effect are equal) receives most evidence; but only approximately 3 times more evidence than the hypothesis H1 (which assumes that gender plays a larger role than race after the teacher effect).In order to draw more decisive conclusions of whether H1 or H2 is true more data would be required.These results are confirmed when looking at the posterior estimates, i.e., 2.22, 0.289, 0.143, and their posterior standard deviations, i.e., 0.344, 0.131, 0.103, of the absolute values of the teacher, gender, and race effect.Note again that the Bayes factors and posterior probabilities are a principled probabilistic methodology to summarize these findings without requiring subjective eyeballing of the estimates.
Under the receiver model, we consider a more exploratory approach by considering all possible order hypotheses on the absolute effects of these nodal characteristics.Note that even though the teacher still has a dominant role, this may not necessarily be the case in the receiver model.We also display the results in the form of an evidence matrix which yields: and translating these to posterior probabilities result in the following posterior probabilities for the respective order hypotheses: 0.291, 0.006, 0.001, 0.001, 0.068, 0.633.These results indicate that almost all probabilities mass goes to H1 and H6, which suggests that the race variable has the smallest impact.Furthermore, the gender variable is expected to play the largest role in the receiver model among these nodal characteristics.These results are confirmed when looking at the posterior estimates, which equal 0.243, 0.310, 0.070, and posterior standard deviations 0.172, 0.100, and 0.045, for the absolute teacher, gender, and race effect, respective.

Testing the impact of network statistics on classroom dynamics
Next we focus on testing the effects of endogenous (network) characteristics of actors.Specifically, we evaluate whether the fixed effects of popularity or outgoingness makes actors more likely to be the next sender.Let φ popularity be the mean effect of popularity and φoutgoingness be the mean outgoingness effect, both in the sender model.We Assuming equal prior probabilities of the three hypotheses, the posterior probabilities are equal to P (H1|E) = 0.384, P (H2|E) = 0.006, and P (H3|E) = 0.610.The results indicate that H3 is most likely to be true after observing the data, being approximately twice more likely than H1, and approximately 100 times more likely than H2.For this reason, we can with almost complete certainty state that the popularity effect is not smaller than the outgoingness effect.This is confirmed by the posterior estimates of the popularity and the outgoingness effect which are equal to 0.17 and 0.07, respectively; see also the posteriors in Figure 10.More data is required in order to obtain more decisive evidence of whether H1 or H3 is true.
In addition to testing the fixed effects as above, we can also test the global means of random effects.Here, we test the participation shift ABBA is more likely than the participation shift ABAB (both are random-effects) after controlling for the other effects in the receiver model.Let µ abba be the global mean effect of ABBA and µ abab be the global mean effect of ABAB.This comparison captures the tension between an individual's tendency to continue to speak to the same individual or the societal norm of reciprocity where the previous receiver because the sender of the next event and the previous sender because the receiver.In polite conversation norms, we would expect to see more turn-switches between two individuals (A speaks to B and then B responds to A) than turn-continuing (where A directs a comment to B and then continues to speak without giving B the opportunity to respond first).The The evidence favoring H3 is overwhelmingly larger than the evidence for any other hypothesis, strongly supporting the hypothesis that the effect of immediate reciprocity is larger than that of immediate inertia across these fifteen classrooms at Magnet High school.This is also confirmed by the posterior probabilities P (H1|E) = 0.000, P (H2|E) = 0.000, and P (H3|E) = 1.000, when assuming equal prior probabilities.
Finally we showcase a hypothesis test of network effects between the receiver model and the sender model.
Specifically, we test whether the teacher variable has a larger effect on the rate of being a sender or on the rate of being chosen as receiver.Given the dominant role of the teacher as a sender, we expect the teacher effect is larger in the sender model.We formulate the hypotheses as follows: Indeed, the evidence matrix shows convincing evidence that the teacher effect is larger in the sender model than in the receiver model.The posterior probabilities of the three hypotheses are equal to P (H1|E) = 0.000, P (H2|E) = 0.000, and P (H3|E) = 1.000, when assuming equal prior probabilities, which results in the same conclusion.The results are also confirmed by the estimates which are equal to ζteacher = 2.22 and μteacher = 0.20.

Discussion
In this paper we presented a Bayesian actor-oriented multilevel relational event model for studying social interaction behavior from independent relational event sequences.This model allows for inferences at the actor level, thus opening the possibility of unveiling effects that make an actor more prone to send or receive an interaction in the population under study.Our results show that the model is able to capture the effects that have the largest impacts in actors preferences, even when those effects are different in the sender activity and receiver choice rates.The models can be estimated using the Stan programming language, the Stan code for the actor-oriented relational event model is available on Github (link suppressed for peer review).
Furthermore, a flexible set of hypothesis testing procedures was proposed for this class of models, facilitating inferences on population parameters of the estimated effects.These tests can be used for testing whether effects should be treated as fixed or random across sequences, for testing the relative heterogeneity of different network effects across sequences, and for testing equality and order constraints on the fixed effects, even on their absolute values.These tests are useful to get a better understanding about the heterogeneity of social interaction behavior across independent relational event sequences.These novel testing procedures can also be applied to other mixed effects models to better understand the heterogeneity in other types of multilevel or clustered data.
The computational issues originating from the inefficiencies that result from the hierarchical structure of the model are the main limitations in the estimation process.We have taken advantage of the multivariate normal structure of the hierarchical prior to induce a non-centered transformation in the random-effects parameters, which is more efficient in practice for the reasons aforementioned.An issue that we have not directly addressed in our paper is the sparsity of relational event data which may cause the geometry of the posterior distribution to be highly complex (Betancourt & Girolami (2015)).The alien form of the likelihood of the relational event model also presents a challenge to the estimation of hierarchical relational event models.An especially attractive next step would to construct alternative representations of the model so that posterior distributions in closed form could be derived, improving efficiency in the estimation process.
Another promising direction to improve computational feasibility when fitting multilevel relational event models to the large clustered relational event sequences is using meta-analytic approximations (Borenstein et al. 2010).A Gibbs sampler could be derived using the point estimates from the independent relational event sequences as data observations.The question would then be how large the relational event sequences should be in order for these meta-analytic approximation to be accurate enough to make reliable statistical inferences.Our flexible Bayes factor tests could then be built on top of that.Finally, another important direction for future research would be to model time-varying coefficients across sequences in order to discover complex temporal changes of the network drivers of social interactions over time.

D Bayes factor for testing random-effect structures
In this appendix we show details on the derivation of the parameters in the Bayes factor for testing random-effect structures.Assuming we have K social networks and each of them has its own β k , for i = 1, . . ., K, regression parameter.Where β k ∼ N (µ, σ 2 ).Thus, let μ, σ2 , βk be posterior estimates for µ, σ 2 and β k , ∀ k.Also, let τ 2 k be the point estimate for the variance of the posterior distribution of β k .Then we can write posterior: prior: β k ∼ N (μ, σ2 ).The parameters in the prior are derived in the same way with E(ξj) = 0, Var(ξj) = 2σ 2 and Cov(ξj, ξ j+h ) = −σ 2 , if h = 1, and zero otherwise.

E Point-wise deviance residuals
After estimation, it is common practice to check how the proposed model fits the data under study.In DuBois et al.
(2013), they used the posterior draws to evaluate the log likelihood at every data point to determine the dynamic ( This measure is usually called residual deviance (Collett 2015), and it is useful to compare models and see which one fits better each data point by having smaller values of Resm.
in a study to investigate student rebellion in the classroom.The data feature observations of interactions among high-school students in two different schools in the United States.For this illustration, we consider 15 independent classrooms from Magnet High School during the 1996-1997 school year.The student body of this high school can be considered academically homogeneous.The data were collected through classroom observations in which the conversations within the classroom were coded.In each of these fifteen classes, a teacher is present in addition to the students.The number of events (the number of times one person said something to another person) ranges from 86 to 628, and the number of persons (the students plus the teacher) ranges between 19 and 30 across the event sequences.The conversations happened in an orderly lecture-like fashion, so only one person was speaking at each time.The aim of this application is two-fold.At first, we present a sequential analysis that is used as proof of concept for the proposed model and statistical testing procedures.Secondly, we aim to provide insights into classroom dynamics by showing the evolution of network effects over time and extensively applying our tests to showcase multiple substantively interesting hypothesis tests.Below we first discuss the actor-oriented multilevel relational event model that we use to analyze the data.

Figure 1 :
Figure 1: Trace plots of some of the random-effect parameters.The number 2 in the subscript shows that the parameters belong to cluster number two.The top (bottom) line shows parameters in the sender (receiver) model.
which ignores the multilevel structure of the data.Figures 4 shows the shrinkage for several effects in the sender and receiver models.The lines were plotted using a grey scale where lighter shades mean smaller sample sizes.In each panel, we have the estimates based on the estimation of the separate classrooms at the top, and the multilevel estimates at the bottom.

Figure 3 :
Figure 3: Random effect estimates for receiver (top) and sender (bottom) models for 5 groups.

Figure 5 :Figure 6 :
Figure 5: The development of different drivers of interaction behavior as class time progresses from 20% to 100% of the observed events in each classroom.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Estimated posterior density and Gaussian approximation for a some fixed effects in the sender model.

Evidence for the H 0 − sender model
Posterior probabilities of the fixed effect hypothesis H 0 in the homogeneity test for different effects.

Table 1 :
Parameter estimates of the mean effects.Numbers in bold are effects that were treated as fixed, either φ or ψ, in the mixed-effects model.Not-bolded numbers represent random-effect means, µ.Numbers between parentheses are 95% credible intervals.Dashes indicate that the covariate was not included in the respective model.
estimates in Table1suggests that this is indeed the case since the mean effect of ABBA is almost four times as large as the mean effect of ABAB, 4.27 and 1.18, respectively.Here we formally test this by considering the following