{'fraction_non_alphanumeric': 0.035664714047098124, 'fraction_numerical': 0.009833674453619403, 'mean_word_length': 4.666977736301827, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Marked Self-Exciting Point Process Modelling of Information Diffusion on Twitter

Information diffusion occurs on microblogging platforms like Twitter as retweet cascades. When a tweet is posted, it may be retweeted and henceforth further retweeted, and the retweeting process continues iteratively and indefinitely. A natural measure of the popularity of a tweet is the number of retweets it generates. Accurate predictions of tweet popularity can assist Twitter to rank contents more effectively and facilitate the assessment of potential for marketing and campaigning strategies. In this paper, we propose a model called the Marked Self-Exciting Process with Time-Dependent Excitation Function, or MaSEPTiDE for short, to model the retweeting dynamics and to predict the tweet popularity. Our model does not require expensive feature engineering but is capable of leveraging the observed dynamics to accurately predict the future evolution of retweet cascades. We apply our proposed methodology on a large amount of Twitter data and report substantial improvement in prediction performance over existing approaches in the literature.


Introduction
The advancement of technology has dramatically changed the ways people connect to each other over the past few years. This contributes to the increasing popularity of microblogging platforms, which integrate the features of instant messaging and blogging, enabling users to conveniently share contents like short sentences, images or videos. Twitter is a microblogging service that allows the users to share information in the form of 140-character messages called tweets. As a tweet is posted by a user, it may be shared by the followers of the original poster through an action known as retweeting, which explicitly refers to the original tweet via its unique identification number, or a retweet. This retweeting process can iterate indefinitely, resulting in a cascade of retweets.
Information diffusion modelling in Twitter has been an active field of research. [28] considered Twitter as a platform used for political deliberation and analyzed tweet sentiments by machine learning to forecast the results of elections. By learning the features, one can also find the likelihood of retweets based on the interestingness of contents [23]. [21] proposed a model to account for the rise and fall of influence propagation whereas [3] recently modelled random series of events prevalent in Twitter by Poissonian and self-feeding processes. These studies, however, do not emphasize on popularity prediction, which is our primary concern.
The popularity of a tweet is naturally measured by the size of the retweet cascade, or the number of retweets it generates. The predictions of tweet popularity are important as they can assist Twitter to rank contents more effectively and facilitate the assessment of potential for marketing and campaigning strategies. As such, many models have been proposed to capture the retweeting dynamics and to predict the popularity. One noticeable work is by [31], who proposed a Bayesian approach which predicts tweet popularity using network information. Models based on the theory of point processes [16,32], which do not require network information, were also shown to have good prediction performances. The model proposed by [22], on the other hand, combines the point process models with feature based approaches to predict the popularity. Other models like the growth-adoption model by [20], the spatial-temporal heterogeneous Bass model by [30], and the concept drift model by [19], were all proposed for the purpose of popularity predictions.
Some models employed on other online social networks (OSNs) with the same purpose of popularity predictions are also relevant as the proposed methodologies may be applicable to predictions on Twitter network. Notably, [1] proposed a dynamic linear regression model to predict the clickthrough rate for Today Module on Yahoo! Front Page. Activities on platforms like Youtube and Digg were also modelled, for instance using linear regression models [26] and classification models [2]. Other related works include the reinforced Poisson model [13] applied on Sina Weibo and the model by [29] that incorporates temporality and seasonality, applied on Flickr image data.
Recently [32] proposed a model termed the SEISMIC (Self-Exciting Model of Information Cascades) to model the retweeting dynamics on Twitter. The model describes the retweet intensity of a tweet, or the expected number of retweets per unit time, as a product of the infectivity of the original tweet and the accumulated excitation effect of all previous retweets. [32] estimated the infectivity as a function of time using a kernel smoothing estimator, and the excitation function, or memory kernel, using a graphical approach under the assumption that some retweeting processes follow an inhomogeneous Poisson process with the excitation function as its intensity function. They also proposed to predict the future popularity of a tweet based on calculating the expected number of future retweets by first assuming that the infectivity remains constant since the censoring time, and a subsequent ad hoc adjustment to the expectation to incorporate the decaying trend of the infectivity. They reported that the predictions of tweet popularity using their approach outperformed those based on competing approaches, under several performance measures. [16] proposed another model termed the TiDeH (Time-Dependent Hawkes) model, which models the retweet intensity similar to the SEISMIC, and estimated the infectivity and memory kernel using similar nonparametric kernel smoothing estimators. [16] fitted a circadian rhythmic function to the nonparametrically estimated infectivity function up to the censoring time, and extrapolated it beyond the censoring time to predict the number of future retweets. With certain choices of the smoothing parameters, the tweet popularity predictions based on the TiDeH model are superior to those based on the SEISMIC, especially on long cascades. However, Kobayashi and Lambiotte's approach requires sufficiently long observation time on a retweet sequence to have reliable estimation of the infectivity function, and the prediction performance depends critically on the window size parameter used in the estimation of the infectivity function.
In this work, we propose a marked self-exciting point process model to capture the retweeting dynamics and to predict tweet popularity. Our model is motivated by the SEISMIC and the TiDeH model, and bears some similarities to them. However, our model has some important advantages. First, the intensity process in our model has a linear form similar to that of the original self-exciting process of [15], and therefore the resulting point process is interpretable as a cluster Poisson process, which means our model can be simulated using a cascading algorithm similar to that used for the efficient simulation of Hawkes processes. Second, the estimation of our model and the assessment of the goodness-of-fit can be implemented using principled approaches from point process theory, and the predictions based on our model can also be done properly by exploiting the probabilistic properties of the model, without resorting to ad hoc assumptions such as those needed by the SEISMIC. Moreover, our model is found to be able to capture the retweeting dynamics and make accurate popularity predictions based on much shorter observation times than those required by the TiDeH model.
The rest of the paper is organized as follows. In Section 2 we describe the tweet data and previous models for tweet popularity predictions which motivated our work, and show how our model is built, how we estimate the parameters and evaluate the goodness-of-fit, and finally how to make predictions. In Section 3 we apply the proposed model to the tweet data and compare the prediction performance of our model with that of the SEISMIC and the TiDeH model, to show the superior performance of our model. Finally, in Section 4, we conclude the paper with a discussion.
2. The data, the model, and the methodology

Twitter data
The data * which motivated our work and which will be used to demonstrate our modelling and prediction methodology is that used recently by [32]. The data contains a total of 166,069 reasonably popular tweets published from October 7 to November 7, 2011, each with at least 49 retweets within seven days of publishing. For each tweet, the data includes the Twitter ID of the original tweeter, the posting times of the original tweet and all the retweets within seven days, and the numbers of followers of the original poster and of the retweeters. Following [32], we use data on the 71,815 tweets published in the first seven days of the study period as training data and the remaining 94,254 tweets published in the next eight days as test data. See Figure 1 for 5 randomly selected retweet cascades from the training data set.  Note that the data lacks the complete Twitter network information, that is, for a retweet, the data only has its publishing time and the number of followers of the retweeting account, without information on whether the original tweet or any previous retweet is being retweeted. This implies that the methodology of [31], which assumes the complete Twitter network information, does not apply here. The contents of the original tweets and of the retweets are also not included in the data, and therefore methodologies depending on features of the tweet contents or features of the posters (other than their numbers of followers) do not apply either.
From Figure 1 we observe that the retweets tend to occur in clusters or bursts. This suggests that self-exciting processes like the SEISMIC of [32] are potentially suitable for such data. In the next subsection we present our model for the retweeting dynamics, which is a marked self-exciting process model similar to that of the SEISMIC.

Model formulation
Let (τ i , m i ), i = 0, 1, . . . be a marked point process where 0 = τ 0 < τ 1 < . . . denote the event times and m 0 , m 1 , . . . denote the respective event marks. In the context of information diffusion modelling on Twitter, the event times shall refer to the retweeting times, except that τ 0 = 0 denotes the posting time of the original tweet, and the event marks refer to the numbers of followers of the retweeting (or tweeting, in the case of m 0 ) accounts. Let N (t) = ∞ i=1 I {τ i t} , t 0 be the associated counting process of retweets, and F = {F t ; t 0}, with F t = σ {N (t), m 0 , (τ j , m j ), j = 1, . . . , N (t)}, be the natural filtration of the marked point process. The (conditional) intensity process of N relative to the filtration F is an F-predictable process λ(t), t 0, such that M (t) = N (t) − t 0 λ(s) ds, t 0 is a mean zero F-martingale. In an informal but intuitive notation, the intensity can be written as λ(t) = E [ dN (t)| F t− ]/ dt, from which we note that the intensity at any time point is the expected number of events per unit time given the history of the process prior to that time point.
As the evolution of a point process over time is fully determined by its intensity process, a commonly used approach to specify a point process model is to specify the form of the dependence of its intensity process on the prior-t history of the process F t− . In particular, the SEISMIC of [32] assumes that the intensity of the retweeting process N (t) takes this form, where p(t) is an unspecified positive function, called the infectivity function, which typically decreases in t, and φ(·) is a positive function called the memory kernel. [32] proposed to estimate the infectivity function p(·) nonparametrically using a kernel smoothing estimator with a triangular kernel. To estimate the memory kernel, they assumed that it is of a power law decaying form and that 15 "carefully chosen" retweet cascades follow inhomogeneous Poisson processes with intensity functions proportional to the memory kernel. They then estimated the parameters using histogram and complementary cumulative distribution function plots of the retweeting times in those 15 cascades. As for the TiDeH model of [16], it assumed an intensity process of the same form as (1), except the further assumption that the infectivity function p(·) is also parametric, and takes a dampened circadian oscillation form. To estimate the infectivity function p(·), [16] proposed a twostep approach where a preliminary estimatep(·) was obtained first using a kernel method and then the parametric form of p(·) was fitted to the preliminary estimate by a least squares method.
The point process model we propose in the current work for the purpose of retweeting dynamics modelling is given by where ν(·) is the baseline intensity function, with ν(t) denoting the part of the event intensity at time t that is due to the initial event at time 0; and ω(·, ·, ·) is the excitation function, with ω(τ, m, t − τ ) denoting the impact of an event at time τ with mark m on the event intensity at time t, where t is the time since the initial tweet was posted. Furthermore, both the baseline intensity and the excitation functions are time-dependent and take multiplicatively separable forms as follows, Here α > 0 is a constant giving the direct excitation effect of the original tweet, that is, how many retweets it is expected to generate directly. The function φ(·) is called the memory kernel function, which describes how the excitation effect due to the original tweet or a retweet is distributed over time. Similar to [32], we require φ(·) to be a probability density function, so that φ(·) 0 and ∞ 0 φ(t) dt = 1. The function p(·) indicates how the "infectivity" of a retweet varies over time and is also called the infectivity function, although its influence on the intensity process λ(t) is different than that of the infectivity function p(·) in (1). For identifiability, we assume that p(0) = 1. The function r(·) is called the impact function, and describes the total excitation effect of a retweet attributed to the number of followers of the retweeter. Note, we do not require α = r(m 0 ), to allow for the potentially different influences of the original tweet and of the retweets.
More specifically, the functions in (3) are assumed to take the following parametric forms, for parameters β 0, γ 0, δ 1 > 1, and δ 2 > 0. Here we have adopted an exponential decay form for the infectivity function, based on the intuition that the infectivity, or newsworthiness of a retweet should decay very quickly over time. We further assume that the impact function is linear in the number of followers on a log scale, rather than on the original scale as in [32], because of the high degree of right skewness for the distribution of the number of followers [4,6,17]. Our choice of the power law decay form for the memory kernel is motivated by [32] and the empirical findings of the heavy tailed distributions for the human response time in social networks, reported in the literature [5,10,31]. Finally, as in [32] and [16], we also assume the event marks m i are i.i.d. with a common density function f (·) relative to a suitable reference measure on the space M of event marks, and moreover, m i is independent of τ i and F τi− for all i. As the excitation function associated with an event is allowed to depend on the time of that event, the model will be called a Marked Self-Exciting Process with Time-Dependent Excitation Function, or MaSEPTiDE for short.
At this point we emphasize an important difference between the MaSEPTiDE we propose in this work and the SEISMIC of [32]. From (2), we note that, unlike the SEISMIC, the MaSEPTiDE has an intensity process that is of a linear form similar to the self-exciting process of [15], whose intensity process takes the form, In fact, if we choose p(τ ) ≡ 1 and r(m) ≡ r for a constant r in (3), then (2) reduces to the timevarying version of the Hawkes process considered by [7,8]. The linear structure of the intensity process implies that the MaSEPTiDE can also be interpreted as a Poisson cluster process, as the original Hawkes process or the generalized version with a time-varying background intensity. By this interpretation, immigrants arrive according to a marked inhomogeneous Poisson process with its intensity function equal to the baseline intensity function ν(·), and event mark distributed according to the density function f (·). Once an immigrant with mark m arrives at τ , it starts to independently produce children according to a marked inhomogeneous Poisson process with intensity function ω(τ, m, ·) = p(τ )r(m)φ(·) and event marks distributed according to f (·), so that the total number of children is Poisson distributed with mean ∞ 0 ω(τ, m, t) dt = p(τ )r(m), and given the total number of children, the waiting times to births of the children are i.i.d. with common density function φ(·) if the order of births is ignored. Moreover, once an offspring of any generation is born, say at time τ and with mark m , it starts to independently produce children of its own according to a similar marked inhomogeneous Poisson process with intensity function ω(τ , m , ·) and event marks distributed according to f (·). The events of the MaSEPTiDE process by time t consist of all immigrants and offspring of any generation that have arrived by time t. This Poisson cluster process interpretation implies an efficient recursive cascading algorithm to simulate the MaSEPTiDE process as described in Section 2.5, which has important implications for simulation based predictions by the MaSEPTiDE process.
Because of the Poisson cluster interpretation, the memory kernel function φ(·) in the MaSEP-TiDE can also be called the offspring density function, and the function p(·)r(·) might be interpreted as the branching ratio function which specifies how the branching ratio, that is, the average number of direct (or generation 1) offspring from an individual (be it an immigrant or an offspring), depends on the birth time and event mark of the individual. In contrast, the functions p(·) and φ(·) in the SEISMIC or the TiDeH model do not permit such a neat interpretation.
It might also be of interest to note the difference between the treatments of the background intensity in the MaSEPTiDE model and in the Hawkes process model with a general time-varying background intensity. In the former model, we require the baseline intensity function to be proportional to the memory kernel φ(·), while in the latter, the background intensity and the memory kernel can take different shapes. The advantage of our treatment is that it leads to a more parsimonious model, while the time-varying background intensity model can easily accommodate nonstationarity, such as that due to the diurnal patterns of human activity levels.

Parameter estimation
Before we can use the MaSEPTiDE model for future events prediction, we need to first estimate the model parameters. Since the event marks are assumed to be i.i.d., their distribution can simply be estimated by the empirical distribution of m i , for i = 1, . . . , N (T ). The main estimation problem is to estimate the parameter vector θ = (α, β, γ, δ 1 , δ 2 ) . To this end, we shall use the maximum likelihood (ML) approach. By the point process theory [11,Proposition 7.3.III], the likelihood of the MaSEPTiDE process based on observations over the interval [0, T ], where T denotes the censoring time, takes the following form where λ(·) depends on the parameters through (2)-(4), and f (·) denotes the event mark density, which is assumed to be free of the parameters θ.
To compute the ML estimator of the parameter vector θ using general-purpose numerical optimization routines, the efficient evaluation of the likelihood function or its logarithm is very important. For this purpose, we need to be able to evaluate the definite integral of the intensity function in (5) efficiently. Due to the linear structure of the intensity function, the integral of the intensity function can be shown to take an explicit form similar to the intensity function itself, and therefore can be exactly computed without resorting to numerical quadrature routines. To show this, it is convenient to use the random measure interpretation of a marked point process. That is, we interpret as a random measure on [0, ∞) × M, so that the intensity in (2) can be written as Therefore, by Fubini's theorem, a change of variables, and the assumed forms of the functions ν, ω and φ, we have where From the separable form of the likelihood function in (5) and the assumption that the event mark distribution does not depend on the parameter vector θ, the ML estimation of the parameters θ can be based on maximizing the logarithm of the part of the likelihood that does not involve f (·), that is In practice, the maximization can be done numerically using various general-purpose optimization routines. In our numerical experiments, we have used the downhill simplex method of [24], which is the default method used by the function optim in the R software environment for statistical computing [27].

Goodness-of-fit assessment
The assessment of the goodness-of-fit of models to historical data can guide us to seek models that can describe the observed data well and therefore serves as the basis of predictions for future observations. To assess the goodness-of-fit of the MaSEPTiDE, we shall use the residual point process approach based on Papangelou's random time change theorem [11,Theorem 7.4.I]. By the time change theorem, with Λ(t) = t 0 λ(s) ds denoting the cumulative intensity process, the transformed process N (Λ −1 (t)) is a Poisson process with unit rate or equivalently, the random times Λ(τ i ), i = 1, 2, . . . , will be the event times of a unit rate Poisson process. Therefore, if the MaSEPTiDE with the parameters θ set to their ML estimatesθ is a sufficient model for the observed event times up to the censoring time T , then the transformed event times,Λ(τ i ), i = 1, . . . , N (T ) should be approximately equal in distribution to the event times of a unit rate Poisson process up to timeΛ(T ), whereΛ(t), t > 0 is the plugin estimate of the cumulative intensity Λ(t; θ) = t 0 λ(s; θ) ds, that is,Λ with p(·) and r(·) defined in (4), and Φ(·) defined as in (6). As the conditional distribution of the event times of a Poisson process in a fixed interval given the total number of events in the interval is equal in distribution to the order statistics of the same number of i.i.d. random variables uniformly distributed in the interval, to assess the goodness-of-fit of the MaSEPTiDE (or any point process model specified via the intensity process), we can assess the uniformity of the transformed event timesΛ(τ i ), i = 1, . . . , N (T ), in the interval (0,Λ(T )], informally using graphical approaches such as the histogram or the Q-Q (quantile-quantile) plots, or formally using statistical tests like the K-S (Kolmogorov-Smirnov) test. A similar residual analysis was performed by [25] to assess the goodness-of-fit of point process models on earthquake data.

Predicting future number of events
Given observations up to T , to predict the number of events from T to a future time pointT > T , one commonly uses its conditional expectation or its conditional median, which is optimal relative to the mean squared error or the mean absolute error accordingly [14]. To obtain the conditional expectation, we can use either a solve-the-equation approach or a simulation based approach. The former approach involves deriving a functional equation satisfied by the conditional expectation as a function of a future time point, solving the equation, and evaluating the solution function at the desired time point. The latter approach involves simulating the sample path of the MaSEPTiDE on the time interval (T,T ] conditional on the observations up to time T for a large number of times, counting the number of events on each simulated sample path, and using the average of the simulated event counts to approximate its expectation. While the first approach is computationally less expensive, the solution of the functional equation is not always easy to obtain. For the second approach, although it is relatively less efficient, especially if the process to be simulated has a large expected number of events, it is more robust than the first approach. To obtain the conditional median, the only option seems to be a simulation based approach, which involves simulating the conditional sample path of the MaSEPTiDE process a large number of times and extracting the median of the resultant empirical distribution of the number of events in the prediction interval. Both solve-the-equation approach and simulation based approach rely on the observation that, conditional on the history of the MaSEPTiDE process up to time T , its future evolution is the same as that of another MaSEPTiDE process with a different baseline intensity function and a similar excitation function. LetÑ (t) = N (T + t) − N (T ), for t 0,τ j = τ N (T )+j − T ,m j = m N (T )+j for j = 1, 2, . . . , andF t = F T +t , t 0. Then, theF-intensity process ofÑ (t) is given bỹ whereν(·) denotes the functioñ andω(·, ·, ·) denotes the functioñ ω(τ, m, t) = ω(T + τ, m, t) = p(T + τ )r(m)φ(t) ≡p(τ )r(m)φ(t).
Therefore,Ñ (t), t 0 is a MaSEPTiDE process with baseline intensity functionν and excitation functionω given as above in (7) and (8) respectively. The excitation functionω has a similar separable form as ω, with r and φ the same as before, and infectivity functionp equal to a time shift of the previous infectivity function, that is,p(τ ) = p(T + τ ).
To calculate the expected number of events E[N (T )−N (T ) | F T ] without resorting to simulations, we first note from the definition of the conditional intensity that, withλ(s) = E[λ(s) | F T ] denoting the mean intensity function ofÑ (t) given F T . By the independence between event marks and previous event times, we havē where we have also usedÑ ( dτ, dm) to denote the associated random measure again, and F denotes the distribution of the i.i.d. event marks, while R = E [r(m i )] = M r(m) dF (m) is the expected total excitation effect due to an event.
In general, we need to solve the integral equation in (10) numerically to obtainλ(t) on [0,T − T ] and use it in finding the conditional expectation of the number of events in (9). One method to solve (10) is to approximateλ(t) by a flexible parametric function and identify the parameters by requiring both sides of the equation to be equal or approximately equal at sufficiently many points in the interval [0,T − T ]. Examples of the flexible parametric functions to approximateλ(t) include a B-spline function with a specified order and knot sequence, or a truncated Fourier series. In both cases, the unknown parameters of the approximating function can be obtained by solving a linear equation of the unknown parameters. In practice, we would try approximating functions with increasing flexibility until convergence in the solution is achieved. We select the B-spline function in this work as a method to findλ(t) for its ease of implementation and computational stability. Specifically, we let B(t) = (B 1 (t), . . . , B k (t)) denote the set of B-spline basis functions of a certain order on the interval (0,T − T ], and assume thatλ(t) ≈ B(t) η for a k-vector η. Plugging this into (10), we have the following equation of η, To solve (11) for η, we evaluate both sides of (11) at sufficiently many ( k) t values in the interval (0,T − T ], and solve the resulting overdetermined linear system using the method of least squares to get η. Once η is obtained, the predicted value is calculated as In evaluating the integrals in (11), we often need to use numerical quadrature routines. In our case, we have used the R function integrate for this purpose.
To simulate the MaSEPTiDE processÑ over the interval (0,T − T ], we can use the following cascading algorithm, which is a generalization of that used for the simulation of non-stationary self-exciting point processes [7,8]. A similar algorithm has been used by [9] to simulate renewal Hawkes processes. The algorithm shown above requires the simulation of inhomogeneous Poisson processes, which can be achieved using the thinning algorithm of [18]. In our numerical experiments, we have used the R implementation simPois from the IHSEP package. Our implementation of the above cascading algorithm is based on a simple modification of the function simHawkes1 from the R package IHSEP. To predict the number of events in the interval (T,T ], we simulate the sample path of the process N (t) over the interval (0,T − T ] for a large number (say 100) of times, and count the number of events on each simulated sample path. The mean or median of these simulated event numbers will then be our point prediction of the number of events of the MaSEPTiDE process N in the interval (T,T ]. In practice, when we use the fitted model to make predictions, whether by using the solvethe-equation approach or by using the simulation based approach, the unknown functionsν and ω, and the event mark distribution F need to be replaced by their respective estimators. In our numerical experiments, we use the plugin estimatorsν(t;θ) andω(·, ·, ·;θ) forν andω, and the empirical distribution functionF of the event marks m 1 , . . . , m N (T ) for F . One implication is that the constant R in (10)

Application to the tweet data
In this section, we report the results of applying the proposed model and inference methodologies to the tweet data. The performance of our prediction methods is also compared to those of the SEISMIC and the TiDeH model.

The model fit
We fitted the MaSEPTiDE model to the 71,815 retweet cascades in the training data set described in Section 2.1 with different censoring times, using the maximum likelihood method described in Section 2.3. The estimated parameter values with censoring time of seven days are highly skewed, with the median estimates of α, β, γ, δ 1 , and δ 2 equal to 48.349, 0.072, 7.209, 1.416, and 0.007 respectively. To have some idea about the typical parameter values found in practice, we display in Table 1 the estimated parameter values for the five sample cascades shown in Figure 1, together with their final popularity.
The estimated values of the parameter β suggest very fast decays of infectivity, with the times taken for the infectivity to drop to 1% of the initial levels vary from about 19 seconds (log(100)/0.246 Table 1 Fitted parameter values on the sample cascades shown in Figure 1 Table 2. From this table we note that the percentage of cascades from which the estimated model passes the test decreases when the censoring time increases. This is to be expected as the amount of data increases with the censoring time, implying that the difficulty of finding a fitting model also increases. At significance level of 0.01, when fitted to the complete cascade data, that is, with the censoring time of 168 hours (seven days), the MaSEPTiDE model passes the goodness-of-fit test on roughly 75% of the cascades. By the censoring time of 12 hours, the MaSEPTiDE model passes the goodness-of-fit test on the majority (approximately 82%) of the cascades. Given that the majority of the retweets, or 80% on average, have already happened within 12 hours of the posting of the original tweets, we conclude that the MaSEPTiDE model is able to describe the retweeting dynamics reasonably well.

Popularity prediction
For each of the 94,254 tweets in the test data set, we applied the fitted MaSEPTiDE model with the retweet cascades censored at different times to predict their final popularity, using the prediction methods discussed in Section 2.5. For the purpose of comparison, we also obtained the predictions based on the SEISMIC of [32] and the TiDeH model of [16]. We only report the results of comparison with these two methods, because they were found to outperform other methods in the literature, such as those reported in [10], [1], [26], and [13], both in our own numerical experiments and in the works of [32] and [16].
Our point prediction of the "final" popularity of a tweet, or the total number of retweets by timeT = 7 days, using the MaSEPTiDE model estimated with the retweet cascade observed up to the censoring time T , is given by N (T ) pred = N (T ) + (N (T ) − N (T )) pred , where (N (T ) − N (T )) pred is obtained either as the conditional expectation using the solve-the-equation approach or the simulation based approach, or as the conditional median using the simulation based approach. Our numerical experiments have confirmed that in the case of conditional expectation, the two approaches produce identical predictions up to a negligible numerical error, as expected. When using the simulation based approach, we note that, for some very popular tweets, the retweeting cascades are very long and the numbers of retweeting events to be simulated in the prediction intervals are very large, and therefore simulations can take a long time to complete. This issue was also noted by [16]. A trick we used to mitigate this issue is to simulate the processÑ with a smaller baseline intensity function, sayν(·)/S with S = 100 or larger, and inflate the simulated event numbers by the factor S. For the majority of the cascades, a moderately large number of simulation replications, such as 100 or even 50, was enough to produce a prediction consistent with that by the solve-the-equation approach. The same set of simulations were used to calculate the median based prediction.
To assess the performance of the conditional mean based predictions, we first follow the recent literature [16,32] and use the Absolute Percentage Error (APE), to compare the accuracy of prediction by different models. Each prediction method under evaluation was applied to each of the 94,254 retweet cascades in the test data set with censoring times T = 2, 4, . . . , 12 hours as most retweeting events would have occurred within the first few hours, based on our analysis on the training data set. For each censoring time, we calculated the APEs of the conditional mean predictions by the proposed model and by the two competing models. The predictions by the SEISMIC approach were calculated using the R package seismic. The predictions by the TiDeH model approach were calculated using the algorithm described in [16] with the window size parameter ∆ obs [cf. 16, p. 194] in the estimation step set to one hour. Due to the lack of a principled approach to select the window size, we chose this value based on experimenting with different values and selecting the one that produced reasonable estimates of the infectivity function by visual inspection. The boxplots of the APEs of conditional mean predictions by the three models (MaSEPTiDE, SEISMIC, and TiDeH) with different censoring times are shown in Figure 2. In each boxplot, the horizontal thick bar indicates the median APE, and the circular point indicates the mean APE.
The actual values of the median and mean APE are given in Table 3. From Figure 2 and Table 3, the MaSEPTiDE prediction has consistently smaller median APE and mean APE for each T than the SEISMIC prediction. Compared to the TiDeH prediction, the MaSEPTiDE prediction has clearly better performances when T = 2, 4, 6 hours, either by the median APE or by the mean APE. The performances of these two models are comparable when T = 8 hours, but the MaSEPTiDE is at a slight disadvantage when T = 10, 12 hours. Despite its widespread use in evaluating prediction performance, the APE as a prediction error measure is not consistent with the feature of the predictive distribution used as a point prediction here, which is the expectation/mean; see [14] for a systematic discussion of this issue. A more appropriate prediction error measure when the conditional mean is used as the point prediction is the squared error. Therefore, we also calculated the squared errors of the predictions by the three models. The boxplots of the squared errors by the three models at different censoring times are shown in the left panel of Figure 3. As the models can occasionally produce extremely large predictions, and even infinity in the case of SEISMIC, the outlying values were not shown in the boxplot for better visualization. The mean squared prediction errors (MSEs) at different censoring times are indicated by the circular points in the boxplots, and their squared roots, that is, the root mean squared errors (RMSEs), are shown in Table 4. From Figure 3 and Table 4, we note that, using the RMSE as the performance measure, the MaSEPTiDE prediction outperforms the SEISMIC prediction at all the censoring times, and similar to the conclusion drawn based on the median APEs, the MaSEPTiDE again, outperforms the TiDeH model when T = 2, 4, 6 hours but slightly underperforms when T = 8, 10, 12 hours. In comparison, the SEISMIC only outperforms the TiDeH model when T = 2, 4 hours. To compare the prediction performances when using the conditional median by different models, we use the mean absolute error (MAE) as the criterion of comparison, as advised by [14]. The conditional median prediction by the MaSEPTiDE model was calculated by the simulation based approach described in Section 2.5. The conditional median prediction by the TiDeH model was similarly calculated using a simulation based approach, although the simulation of the TiDeH model was achieved by a less efficient method where the events have to be simulated serially one after another, using the rejective method of [18]. The conditional median prediction by the SEISMIC is not included in this comparison because this model does not specify the form of the intensity process beyond the censoring time, and therefore we cannot calculate the conditional median using the simulation based approach. The right panel of Figure 3 shows the absolute errors of the conditional median predictions by the MaSEPTiDE and the TiDeH models at different censoring times, where, as before, the circular points indicate the MAEs of the predictions at the corresponding censoring times. See also Table 4 for MAE values. By the MAE, the conditional median prediction by the MaSEPTiDE model is clearly superior to that by the TiDeH model at the censoring times T = 2, 4, 6 hours, and is comparable albeit slightly inferior at the larger censoring times T = 8, 10, 12 hours.
By all the performance evaluation criteria considered, the prediction by the MaSEPTiDE model is clearly more accurate than by the two competing models, especially when prediction needs to be made based on shorter observation times, for example, within six hours or shorter of the posting of the original tweet.

Conclusion and discussion
In this work, we have proposed a marked self-exciting point process model, called the MaSEPTiDE, to model the retweeting dynamics and to predict tweet popularity. The MaSEPTiDE is able to model a large number of retweet cascades adequately, and its prediction performance is superior to those of the competing models and approaches in the literature that require the same input.
When the prediction is based on observing a cascade for a long period of time, the approach based on the TiDeH model by [16] is found to outperform our model by a small margin. However, considering the fact that this small advantage of the TiDeH model is not realized until the retweet cascades are observed for eight hours or longer, when the majority of the retweets would have already happened, its practical significance is rather limited. On the contrary, the approach based on the MaSEPTiDE is able to provide accurate predictions of the final popularity based on observations within two hours of posting of the original tweet. Another issue with the TiDeH model is that the nonparametric estimation step to obtain the initial raw estimate of the infectivity curve needs a large amount of data to work well. In fact, in their numerical experimentation, [16] only verified the superior performance of their prediction approach relative to the SEISMIC on 738 very long cascades (containing 2,000 or more retweets), which account for less than 0.5% of all the cascades. In contrast, the approach based on the MaSEPTiDE is fully parametric and does not require as much data to estimate.
The specific parametric forms of the functions in the MaSEPTiDE model have been selected from a class of candidate models by comparing their goodness-of-fit on the retweet cascades in the training data set and identifying the model that can fit most of the cascades. In the class of candidate models, we have considered other parametric forms of the component functions, such as infectivity functions that decay at polynomial rate, and memory kernel functions that decay exponentially fast. The model with the specific forms of the component functions reported herein has the best goodness-of-fit on the training data.
To further improve the MaSEPTiDE model, more complex models, such as those that incorporate the calendar time effects [12,16] are worth considering. Another aspect of our approach that can be improved is that our approach still requires the observation of the retweet cascade for a substantial amount of time to accumulate enough data to identify the model, even though the required observation time is much less compared to approaches based on other models such as the TiDeH model. If we make stronger assumptions on the model parameters across the cascades, then parameter estimation might be achieved using only training data, which will allow us to predict its final popularity as soon as a tweet is published, or even before it is published.
Finally, an important limitation of the data considered in our work, originally collected by [32], is that the data contains only cascades with at least 49 retweets. Such data is by no means representative of all tweets published by Twitter users, as vast majority of the tweets do not get even a single retweet. Therefore, models developed based on such data are only useful for popularity predictions of reasonably popular tweets. To develop models suitable for the predictions of the popularity of average tweets, one would need to collect suitable random samples of tweets and their retweet cascades, and build models accordingly.