Generalized Poisson Dynamic Network Models

Generalized Poisson Dynamic Network Models

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 {

} \fi

\if01 {

} \fi

{\it Keywords:} Overdispersion, Concentration Inequalities, Bayesian Inference, Latent Space, Identifiability.

Introduction

Temporal networks have attracted considerable attention across a wide range of disciplines holme2012temporal, including biology pastor2015, neuroscience betzel2017generative,vavsa2022null, economics jackson2002evolution, friel2016interlocking, and the social sciences barbera2015birds, casarin2025media. The distinctive feature of temporal networks is that connections between nodes are not static, but change over time.

In recent years, an increasing number of studies have focused on count-weighted temporal networks, that is, dynamic networks in which the weight associated with each edge represents a count sewell2016latent. Examples include online media interaction networks casarin2025media, email communication networks yin2017local, transportation networks he2025semiparametric, and brain connectivity networks zhang2019nonparametric. It is well established that count-weighted network data frequently exhibit overdispersion zheng2006many, corsini2022dealing and underdispersion lord2010statistical,lux2020distribution. Nevertheless, many modeling approaches often overlook this feature, which can lead to biased estimates and misleading inferences.

The contribution of this work is manifold. First, we propose a flexible model for weighted networks with integer-valued weights based on the Generalized Poisson (GP) distribution family and derive theoretical properties of the GP network model, including expected strength and node centrality. Second, we introduce three dynamic specifications for our Generalized Poisson network model. Third, a Bayesian inference framework is provided together with an efficient posterior approximation procedure. The proposed approximation procedure is evaluated through an extensive simulation study, followed by two empirical applications. The first application focuses on bike-sharing data from New York City citibike2019, while the second examines the evolution of media interaction networks casarin2025media. The analysis of these data can have a significant impact across many fields and aspects of society.

Regarding the weight distribution, we assume a GP. This family, introduced by ConJai1973,consul1973some, can capture both over- and underdispersion, excess kurtosis, and includes the Poisson distribution as a special case. See ConFam2006 and FamLee2020 for an introduction. The GP distribution is well studied in the literature and has found applications across diverse fields, including health and epidemiology Zam2016, sport statistics ShaMoy2016, as well as economics, finance and insurance WanFam1997,Ambagaspitiya_Balakrishnan_1994,FemCon1995,FamEtal2004,Lin2004,finner2015some.

GP has been successfully used in count time-series analysis zhu2012modeling,ShaMoy2014, CheLee2016,CARALLO20241359 and naturally emerges in random graph theory from branching processes aldous1998tree,aldous2004tractable,bertoin2012fires. Unlike other distribution families that allow both overdispersion and underdispersion, it preserves tractability in deriving relevant random network properties. Other popular models for count data, such as the Negative Binomial, Poisson-gamma, and Poisson-lognormal, have limitations: they allow only overdispersion and struggle to handle low sample means and small samples. The Conway-Maxwell-Poisson is also appealing, since it includes the Poisson as a special case, belongs to the exponential family, and admits conjugate distributions boatwright2006conjugate, shmueli2005useful. However, its moments are not tractable, which prevents its use in the analytical study of the network model properties.

Regarding the dynamic specification of our model, we propose three alternative formulations that capture temporal dependence in distinct ways. The first specification introduces a common dynamic latent factor that affects all edges simultaneously BrauningKoopman2020. The second specification is more parsimonious and incorporates a lagged measure of network strength. In this setting, past global network features influence current link intensities, producing autoregressive dynamics in the network tsikerdekis2021network,jiang2023autoregressive. Finally, the third specification assumes a Latent Space (LS) model hoff2002latent with time-varying latent node coordinates that drive the connectivity dynamics sewell2016latent, rastelli2016properties, artico2023dynamic, casarin2025media. Some theoretical properties of the models, such as node centrality, are derived by building on concentration inequalities vershynin2018high, which show how the dispersion parameter affects the random network's connectivity.

We adopt a Bayesian inference framework, as it provides greater flexibility for working with nonlinear and latent-variable models via the data-augmentation principle and the complete-data likelihood robert2007bayesian and accommodates uncertainty in the prediction (mccabe2011efficient). We provide sufficient conditions for latent variable identifiability and an efficient Markov Chain Monte Carlo sampler for the posterior distribution RobCas2013.

The Citibike citibike2019 and Media network datasets schmidt2018polarization provide two examples of unequal dispersion and edge persistence in temporal networks. These two applications show that GP network models fit better than Poisson models and capture salient global and local features of the networks considered.

The work is organized as follows. Section (ref) presents the Citibike and Media Network datasets, and illustrates the dispersion features in the network weights. Section (ref) introduces the proposed Generalized Poisson network model, provides some model properties, and a Bayesian inference procedure. Section (ref) presents the results of a simulation study. Section (ref) illustrates the model performance on the two network datasets. Finally, Section 5 concludes.

Overdispersion in Dynamic Networks

We consider two datasets relevant to many fields and with impact on society: the Citibike dataset citibike2019 and the Media network dataset. schmidt2018polarization. While bike-sharing data, as the former, are of interest in transportation noland2016bikeshare, geography an2019weather, urban planning yu2018cost, and sustainability chen2022environmental, media networks, as the latter, are relevant to media and communication studies xu2022evolution, political science eady2025news, and computational social science del2017mapping, cinelli2021echo. Both datasets have recently attracted the attention of the network statistical modeling literature he2025semiparametric, casarin2025media. The Citibike dataset contains information on rides between any two Citibike stations in New York City and allows obtaining count-weighted symmetric temporal networks with nodes representing the neighborhoods of New York City and edges representing the count of bike-sharing connections between neighborhoods. The Media network dataset contains information on the Facebook activity of national and local news outlets in France, Germany, Italy, and Spain, and allows obtaining temporal networks in which nodes represent the news outlets' pages and edge weights denote the number of unique users who commented on a pair of outlets within a given time interval.

The left panel in Figure (ref) reports the comparison between log-mean and log-variance for each edge of the Citibike and Media network dataset for Italy for two sub-periods. Compared with the 45-degree line, we observe some heterogeneity: a few edge weights exhibit underdispersion, while most exhibit overdispersion. Moreover, there is evidence of a relationship between log-dispersion and autocorrelation at lag 1 for the edges (right panel). This behavior is consistent across sub-periods. The preliminary results indicate the need to use a conditional distribution for network edges that accommodates unequal dispersion and dynamic features. We report in the Supplement (Section (ref)) further results which corroborate this evidence, not only for the Citybike and Italian Media network datasets but also for the Media networks of France, Germany, and Spain.

GP Dynamic Network Models

The GP is defined, and some properties of the GP networks are derived. GP model classes are introduced together with a suitable inference procedure. The proofs of the results in this section are in the Supplement (Section (ref)).

Generalized Poisson networks

Let $\mathcal{G}=\{\mathcal{G}_{t},\, t=1,2,\ldots, T\}$ be an undirected and count--weighted temporal network with $\mathcal{G}_{t}=(V,E_{t}, Y_{t})$ where $V$ denotes the vertex set, $E_{t}$ the edge set and $Y_t = \{Y_{ijt}\}_{i,j= 1}^N$ the corresponding $N \times N$ count--weighted adjacency matrix at time $t$. We consider a GP specification for modeling the count of interactions $Y_{ijt}$ between nodes $i$ and $j$ at time $t$, and we assume conditional independence given a set of latent parameters. The conditional distribution of $Y_{ijt}\in\mathbb{N}$ follows the GP: $$ p(y_{ijt} \mid \lambda_{ijt}, \theta) = \frac{\lambda_{ijt} \left( \lambda_{ijt} + \theta y_{ijt} \right)^{y_{ijt}-1} \exp(-(\lambda_{ijt} + \theta y_{ijt}))}{y_{ijt}!}, \,\, y_{ijt} = 0,1,2,\ldots, $$ where $\lambda_{ijt} > 0$ governs the mean intensity of the process and $\theta \in (-1,1)$ controls the degree of dispersion ConJai1973. We write $Y_{ijt}\sim\mathcal{GP}(\lambda_{ijt},\theta)$. For $\theta = 0$ the model reduces to the standard Poisson distribution, while $\theta > 0$ ($\theta < 0$) allows for overdispersion (underdispersion). Let $\mathcal{F}_{t}=\sigma(\{Y_{ijs},i=1,\ldots,N,j=1,\ldots,N,j\neq i, s=1,\ldots,t\})$ denote the $\sigma$-algebra generated by the collection of edge weights up to time $t$, with $Y_{ijt}\sim\mathcal{GP}(\lambda_{ijt},\theta)$. Following Ambagaspitiya_Balakrishnan_1994 the conditional moment generating function (mgf) of $Y_{ijt}$, defined as $M_t(u)=\mathbb{E}(Y_{ijt}|\mathcal{F}_{t-1})$, is given by

$u < u_0(\theta)$, where $ z(u) = -\theta \exp(u-\theta)$, and $W(z)$ is the Lambert's W function mezo2022lambert. The conditional expected value and variance of $Y_{ijt}$ are given by \[ \mathbb{E}[Y_{ijt}|\mathcal{F}_{t-1}] = \frac{\lambda_{ijt}}{1 - \theta}, \qquad \mathbb{V}(Y_{ijt}|\mathcal{F}_{t-1}) = \frac{\lambda_{ijt}}{(1 - \theta)^3}. \] Since the probability of a missing edge between two nodes $\mathbb{P}(\{Y_{ijt}=0\}|\mathcal{F}_{t-1})=1-\exp(-\lambda)$ does not depend on $\theta$, then the unweighted random graph, where there is an edge between $(i,j)$ if and only if $Y_{ijt}>0$, has the same properties as an Erd\"os-R\'enyi graph with edge probability $p=1-\exp(-\lambda)$. The dispersion parameter $\theta$ plays a role when weights are used in graph metrics. For instance, since the GP satisfies the convolution property Con1989, the total strength of the network follows a GP with dispersion parameter $\theta$.

From the previous proposition, one can deduce that larger values of $\theta$ provide an increase in the expected strength and consequently in the average node centrality. In the underdispersion setting, that is $\theta\in(-1,0)$, the node strength is smaller than in the oversdisperion case $\theta\in(0,1)$. Also, since the GP is a sub-exponential distribution, Bernstein concentration inequalities can be applied to study how indirect global connectivity is affected by dispersion features vershynin2018high.

The following proposition states the sub-exponential property of a GP variable and provides a useful bound for the logarithm of the moment generating function.

Before stating the main result on the concentration for the spectral radius $\varrho(Y_t)$ of the random matrix $Y_t$, we introduce some notation and a preliminary result. Denote with $||X||_{op}=\sup\{||Xz||_2:||z||_2=1\}$ the operator norm of a $(N\times N)$ matrix $X$, with $||X||_2$ the square norm, and with $||X||_{\Psi_1}=\inf\{u>0|\mathbb{E}(\exp(||X||_{op}/u))\leq 2\}$ the Orlicz norm. Note that if $X$ is symmetric then $||X||_{op}=\max\{|\varrho_{i}(X)|,\,i=1,\ldots,N\}=\varrho(X)$, where $\varrho_{i}(X)$ are the eigenvalues of $X$ and $\varrho(X)$ its the spectral radius. The following propositions establish some bounds for the Orlicz norm $||\cdot||_{\psi_1}$ of GP variables.

This bound says that the spectral radius of the random matrix $Y_{t}$ is sharply concentrated around its expectation’s eigenvalue $\varrho(\Lambda_t)$, implying that some graph properties, such as diffusion rate and eigenvector centrality, can be predicted using the expected connectivity $\Lambda_t$. In addition if $\varrho(\Lambda_t)>> \sqrt{v_t^2 \log N}+K_t \log N$ than $\varrho(Y_t)/\varrho(\Lambda_t)\overset{p}{\rightarrow} 1$ as $N\rightarrow+\infty$. The larger the value of $\theta$, the larger are the expected spectral radius and the absolute fluctuation bound. An increase in $\varrho(\Lambda_t)$ implies highly connected nodes become more central, and every node’s potential to propagate influence (in the eigenvector sense) increases.

The overdispersion parameter is not only central to the properties of the random graph model at each point in time, but also crucial for making inferences about network dynamic properties. As shown in the following sections, assuming equal dispersion, $\theta=0$, a priori, can produce significant misspecification bias.

Dynamic specifications

For a better interpretability, we reparametrize the Generalized Poisson model in terms of location and dispersion parameters, that is $\mu_{ijt} = \mathbb{E}[Y_{ijt}]$ and $\rho= \mathbb{V}(Y_{ijt})/\mathbb{E}[Y_{ijt}] = (1-\theta)^{-2}$ thus $\lambda_{ijt}=\mu_{ijt}\rho^{-1/2}$ and $\theta=1-\rho^{-1/2}$ CARALLO20241359. Since $\rho = (1-\theta)^{-2} \in (\frac{1}{4}, \infty)$ for $\theta \in (-1,1)$, we reparametrize further the model and assume the dispersion ratio to be guided by the parameter $\zeta$, i.e. $\rho = \tfrac{1}{4} + \exp(\zeta)$.

The GP parameter is then given by $\theta = 1 - 1/\sqrt{\rho}$, so that overdispersion relative to the Poisson model ($\theta = 0$) corresponds to $\zeta > \log(3/4)$. When the goal is to capture the dynamics in temporal network data, additional features can be incorporated to improve predictions. In the following, we introduce three model classes.

In the first model class, denoted $\mathcal{M}_1$, latent factors are used to explain the formation of edge weights. As an example, consider a dynamic factor $f_{ijt}$ in the log-mean parameter: $$ \log \mu_{ijt} = \alpha_i + \alpha_j + f_{ijt}, $$ where $\alpha_i$ and $\alpha_j$ are individual-specific node effects that account for systematic heterogeneity in nodes’ propensity to form ties -- i.e., baseline differences in connectivity that translate into higher or lower expected degree. These effects capture node centrality, and are often interpreted as node popularity in the network literature, since larger values correspond to nodes that are more likely to attract links across the network friel2016interlocking. A common latent trend is useful because many networks are driven not only by node-specific heterogeneity and pairwise structure, but also by system-wide forces that affect all nodes simultaneously (e.g., overall growth or contraction in activity, macroeconomic conditions, regulatory changes, platform-wide shocks). If such global variation is ignored, the model may attempt to absorb broad changes in link density into the node effects $\alpha_i$ or the latent positions, thereby complicating interpretation and potentially biasing inference. Introducing a common time component separates global intensity from relational structure, improving fit and identifiability.

To capture time variation shared across the whole network, we therefore introduce a common latent local--level component $f_t$ in the link-formation process. This term represents unobserved time-specific conditions that increase or decrease the propensity to form ties for all node pairs. Accordingly, we assume $f_{ijt}=f_t$ for all $i,j$, so that pairwise link probabilities co-move over time through a single scalar factor. We model $f_t$ as a random-walk (local--level) process, \[ f_t = f_{t-1} + \varepsilon_t,\qquad \varepsilon_t \sim \mathcal{N}(0,\sigma_\varepsilon^2), \] which yields a flexible evolution over time: the level can drift gradually, while $\sigma_\varepsilon^2 > 0$ controls the smoothness of the trend (small values imply a slowly varying baseline density, whereas larger values allow sharper shifts). Models in this class require a careful inference procedure with a non-negligible computational cost. This specification relates to state-space models commonly used in time-series analysis xing2010statespace, mazzarisi2020dynamic, BrauningKoopman2020, BuccheriMazzarisi2024.

In the second model class, denoted $\mathcal{M}_2$, past edge weights affect the edge-formation process. We consider a parsimonious autoregressive formulation with $f_{ijt}=f_t$, where $f_t = \sum_{\ell=1}^{k}\delta_\ell \log(\tilde{y}_{t-\ell}),\,\, \tilde{y}_{t-\ell} = \frac{\bar{S}_{t-\ell}}{(N-1)} $ with $\delta_\ell$, $\ell=1,\ldots,k$ autoregressive coefficients and $\bar{S}_{t-\ell} = N^{-1}\sum_{i=1}^{N}\sum_{j\neq i}^N y_{ij,t-\ell}$ is the average network strength at time $t-\ell$, $\ell=1,\ldots,k$. This specification shares points with time-varying ERGMs, in which global network features are employed to predict the network's average behavior shang2009exponential,tsikerdekis2021network,jiang2023autoregressive. The framework can be readily extended to accommodate general INGARCH specifications zhu2012modeling, in which lagged $f_{t}$ values enter the dynamics.

The third class of models, denoted by $\mathcal{M}_3$, is obtained by augmenting the edge specification with time-varying node-specific features (or latent coordinates) $\mathbf{x}_{it}\in\mathbb{R}^d$, $i=1,\ldots,N$. In this setting, the pairwise propensity to form a link is allowed to depend on how “close” two nodes are in the feature space, namely \[ f_{ijt} = f_t - d(\mathbf{x}_{it},\mathbf{x}_{jt}), \] where $d(\cdot,\cdot)$ is a user-chosen function that quantifies similarity or dissimilarity between $\mathbf{x}_{it}$ and $\mathbf{x}_{jt}$. The choice of $d(\cdot,\cdot)$ directly controls how differences in node features translate into link formation, and hence the flexibility of the model in capturing latent clustering and other geometric network effects. Different distance measures are used because the notion of “closeness” between node features depends on the application and on the type of latent representation one wishes to impose.

For example, one may adopt a similarity measure such as the inner product $d(\mathbf{x},\mathbf{y})=-\mathbf{x}^\top\mathbf{y}$, $\mathbf{x},\mathbf{y}\in\mathbb{R}^d$, so that more aligned feature vectors increase the propensity to link. Similarity-based measures emphasize alignment in direction and allow for assortative patterns, in which ties are more likely between nodes with concordant feature profiles. Alternatively, one may use a dissimilarity measure, such as those based on Euclidean distances, such as $d(\mathbf{x},\mathbf{y})=\|\mathbf{x}-\mathbf{y}\|_2$ or $d(\mathbf{x},\mathbf{y})=\|\mathbf{x}-\mathbf{y}\|_2^2$ with $\|\cdot\|_2$ denoting the Euclidean norm, so that nodes that are closer in the latent space are more likely to be connected. Euclidean distances encode homophily by penalizing separation in the latent space, whereas squared distances further strengthen the penalty for large separations and are often computationally convenient (they are smooth everywhere and avoid the singularity at zero distance), whereas the Euclidean norm yields a more direct geometric interpretation in terms of straight-line distances. More generally, the chosen metric determines the geometry of the latent space and thus the class of network structures that the model can represent (e.g., clustering, transitivity, or the presence of hubs). In the following, we focus on the squared Euclidean distance $d(\mathbf{x},\mathbf{y})=||\mathbf{x}-\mathbf{y}||_2^2$, nevertheless our modeling strategy and inference procedure can be readily extended to other distance functions.

The features can be observable or latent, node-specific or common, with a dynamic specification. In the following, we assume a random walk for the common factor $f_t = f_{t-1} + \varepsilon_t$ with $\varepsilon_t \sim \mathcal{N}(0, \sigma_\varepsilon^2)$, and a random walk specification of the latent positions, that is \[ \mathbf{x}_{it}=\mathbf{x}_{it-1}+\boldsymbol{\nu}_{it},\, \boldsymbol{\nu}_{it}\overset{iid}{\sim}\mathcal{N}_d(\mathbf{0},\Sigma_x), \] with $\Sigma_x = \sigma^2_xI_d$. This allows the network to be represented in a latent space - an approach introduced in the latent position model of hoff2002latent, where nodes are embedded as points in a low-dimensional Euclidean space and tie probabilities depend on distances between points, and later extended by sewell2016latent to accommodate time-varying latent coordinates for dynamic networks. Although the modeling strategies discussed above could be further combined, we leave this extension to future research, as our goal is to demonstrate the contribution of GP dispersion to network model properties across classes.

Prior distribution and model properties

We propose a Bayesian inference framework that naturally accounts for parameter uncertainty. Our model choice is to assume independent Gaussian prior distributions on $\alpha_i$ and $\delta_\ell$, that is

as it is standard in relational data models, see hoff2005bilinear. Following CARALLO20241359, a Gaussian prior distribution is assumed for the dispersion parameter $\zeta$, which induces a weakly informative prior on the dispersion ratio $\rho$ and, consequently, on $\theta$

As is common in random network literature rastelli2016properties, casarin2024dynamic, casarin2025media, we investigate the impact of node random effects on the conditional strength distribution, its expected value, and dispersion index. Let $\bar{S}_t$ be the network average strength at time $t$ and $\mathcal{F}_{t}^{\ast}=\mathcal{F}_{t}\vee \{f_{s},i=1,\ldots,N,j=1,\ldots,N,j\neq i, s=1,\ldots,t\}$ the augmented sigma algebra. Note that for model class $\mathcal{M}_1$, $f_t,f_{t-1}\in\mathcal{F}_{t-1}$ and for $\mathcal{M}_{2}$ and $\mathcal{M}_3$, $f_t\notin\mathcal{F}^{\ast}_{t-1}$ and $f_{t-1}\in\mathcal{F}^{\ast}_{t-1}$.

Proposition (ref) allows deriving the following properties for $\mathcal{M}_1$ and $\mathcal{M}_2$; similar properties can be derived for $\mathcal{M}_3$.

The average strength can also be derived for model $\mathcal{M}_3$ following similar arguments as in the proof of Prop . 1 of casarin2024dynamic and of Proposition 2.1 and Cor. 2.1 of casarin2025media.

The individual effect variance, the factor level, and the variance of the factor level can all contribute to overdispersion. As $\sigma^2_\alpha$ increases, one obtains larger overdispersion for $\theta\in(0,1)$ and reduced overdispersion for $\theta\in(-1,0)$. One obtains equal dispersion when $\theta=0$ and $\sigma^2_\alpha$ goes to zero. If $\theta\in(-1,0)$, the model can also account for underdispersion as illustrated in Figure (ref). In the left plot, each line shows the values of $\theta$ and $\sigma_\alpha^2$ such that the dispersion index is 0.5 (solid lines) or 1.5 (dashed lines). Also, negative values of $\theta\in(-1,0)$ are compatible with overdispersion (dashed lines) for larger values of the individual effect variance for different network sizes (left) and factor levels (right). The left plot also shows that underdispersion cannot always be achieved if the value of $\sigma^2_\alpha$ is too large, since one needs $\theta>-1$ to have a well-defined GP distribution. The right plot shows that for larger values of $\sigma^2_\alpha$, $\theta$ must be small to guarantee a given level of average strength across different levels of the factor.

Parameters Identifiability

We note that the three proposed models face the same identification issues as other models in the same class. As concerns $\mathcal{M}_1$, the likelihood remains invariant when we sum and subtract the same quantity from $\mu_{ijt}$. This affects the identification of $\alpha_i$ for $i = 1, \ldots, N$ and $f_t$ for $t = 1, \ldots, T$. This issue is per-se mitigated by choosing a prior centered at zero for either $\alpha_i$ or $f_t$, but can be solved by imposing a zero-sum restriction for either of the two sets of parameters. In the following result, we show how a zero-sum restriction on the vector of parameters $\boldsymbol{\alpha}$ is sufficient to guarantee identifiability.

Regarding $\mathcal{M}_2$, we notice that a further assumption is required on the design matrix $W =(\boldsymbol{w}_1',\ldots,\boldsymbol{w}_T')' \in \mathbb{R}^{T\times k}$ with $\boldsymbol{w}_t = (\tilde{y}_{t-1}, \ldots, \tilde{y}_{t-k})$. One needs $W$ to be pre-determined and full-ranked.

In addition to the aforementioned issues, $\mathcal{M}_3$ suffers from well-known indeterminacy problems related to the likelihood invariance of the latent coordinates in terms of translation, reflection and rotation (see hoff2002latent), which can be solved via a combination of dimension-specific zero-sum restrictions and with the use of Procrustes transformation as stated in the following.

Posterior Approximation

We derive the full conditional distributions used in the Gibbs sampling posterior approximation and conduct simulation exercises to assess the efficiency and effectiveness of our algorithm. Moreover, we demonstrate numerically that substantial misspecification bias and prediction errors arise when a Poisson model is applied to data with unequal dispersion.

Gibbs sampler

Since the posterior distribution is not tractable, Bayesian inference is performed through Markov-Chain-Monte Carlo (MCMC) sampling, in particular via Metropolis-within-Gibbs steps chib1995understanding. Posterior inference is conducted via a Metropolis-within-Gibbs algorithm, which alternates between draws from the full conditional distributions of the model parameters. To simply the notation, we define the following collections of variables and parameters: $\mathbf{Y} = \{Y_1, \ldots, Y_T \}$ $\mathbf{X} = \{X_1, \ldots, X_T \}$ with $X_t = (\mathbf{x}'_{1t}, \ldots, \mathbf{x}'_{Nt})' \in \mathbb{R}^{N \times d}$, $\boldsymbol{\alpha} = (\alpha_1,\ldots, \alpha_N)$, $\boldsymbol{\delta} = (\delta_1,\ldots, \delta_k)$, and $\mathbf{f} = (f_1, \ldots, f_T)$. At each iteration $h=1,\dots,H$, the following updates are performed.

In the first block of updates, sample the individual effects in $\boldsymbol{\alpha}$ by drawing each $\alpha_i$ from its full conditional distribution:

using a random-walk Metropolis-Hastings update.

In the second block of the Gibbs sampler, the parameters and the latent variables of the dynamic components are sampled. Depending on the chosen specification $\mathcal{M}_j$.

1) In the autoregressive specification $\mathcal{M}_1$, update the autoregressive coefficients $\boldsymbol{\delta}=(\delta_1,\ldots,\delta_K)$ from

typically via a multivariate random-walk Metropolis-Hastings step.

2) In the latent factor specification $\mathcal{M}_2$, sample the latent factor sequence $\mathbf{f}$ from its state-space posterior

assuming $f_0 = 0$ and $\sigma_{f_0}^2$ large, while using a random-walk Metropolis-Hastings step. Sample $\sigma^2_\epsilon$ exactly from its inverse gamma posterior: $\sigma^2_{\epsilon}\mid \mathbf{f} \sim \mathcal{IG}(\overline{a}, \overline{b})$, where $\overline{a} = \underline{a} + (T-1)/2$ and $\overline{b} =\underline{b} + \sum_{t = 2}^T(f_{t}-f_{t-1})/2$.

3) In the dynamic latent position specification $\mathcal{M}_3$, first, draw the latent factor sequence $\mathbf{f}$ from its state-space posterior $p(\mathbf{f}|\mathbf{Y},\boldsymbol{\alpha}, \boldsymbol{X},\zeta,\sigma_\varepsilon^2)$ similarly to what was done in $\mathcal{M}_2$, and secondly draw the latent coordinates $\mathbf{x}_{it}$ $i=1,\ldots,N$ $t = 1,\ldots, T$ from the joint full conditional distribution

assuming $\mathbf{x}_{i,0}$ equal to the zero vector. We use a Metropolis-Hastings algorithm with a proposal derived from a log-Taylor expansion of the likelihood.

At a given Gibbs iteration $j$, when sampling $\mathbf{x}_{it}^{(j)}$, the log-Taylor approximation is centered at the previous iteration value of the coordinate, that is $\tilde{\mathbf{x}}=\mathbf{x}_{it}^{(j-1)}$.

In the last block of the Gibbs sampler, the dispersion parameters $(\zeta, \rho, \theta)$ are sampled conditionally given the first and second block, and the other observations. First $\zeta$ is updated (and hence $\rho = \tfrac{1}{4}+\exp(\zeta)$ and $\theta = 1 - 1/\sqrt{\rho}$) by sampling from \[ p(\zeta\mid \mathbf{Y},\boldsymbol{\alpha},\boldsymbol{\delta},\mathbf{f},\mathbf{X}) \propto \mathcal{N}(\zeta;\mu_\zeta,\sigma_\zeta^2)\prod_{t=1}^T\prod_{j > i}\mathcal{GP}(Y_{ijt}\mid \mu_{ijt},\theta(\zeta)). \]

A further step can be added to provide intensity posterior point forecasts and network distribution forecasts via an MCMC sampling scheme, following CARALLO20241359 and CheLee2016.

We implement our MCMC algorithm in C++, leveraging the Rcpp package eddelbuettel2013seamless for easy integration with the R environment. For the random number generator and likelihood function of the GP distribution, we adapted the code developed in the VGAM package yee2015vector.

Approximation efficiency and misspecification bias

We assess the performance of the proposed estimation procedure through Monte Carlo simulations. Data are generated from three alternative specifications ($\mathcal{M}_1$--$\mathcal{M}_3$) under a Generalized Poisson framework, incorporating dynamic latent factors, autoregressive dynamics, and time-varying latent coordinates. Across all designs, the MCMC algorithm (5,000 iterations, 2,000 burn-in, thinning every 5 draws) shows good mixing and convergence, as confirmed by graphical diagnostics, effective sample sizes, and Geweke tests. Posterior distributions accurately recover the true structural and dispersion parameters, with only minor bias in the most complex latent-coordinate specification. Overall, the results indicate that the proposed methodology reliably recovers the underlying data-generating parameters. See Section (ref) in the Supplement for full details on the simulation exercise.

Omitting underdispersion or overdispersion features in the model can lead to significant estimation bias in the parameters and imprecise uncertainty quantification. Figure (ref) reports the posterior distribution (light blue) and the true value of some characteristic parameters for $\mathcal{M}_1$, $\mathcal{M}_2$ and $\mathcal{M}_3$ under Poisson (upper panel) and GP specification (bottom panel) when the data generating process assumes a GP likelihood. The bias in the estimates induced by the misspecified Poisson model clearly emerges (top vis-a-vis bottom panels).

We compare the Deviance Information Criterion gelman2014understanding of the correctly specified Generalized Poisson models and the misspecified Poisson models under the three specifications $\mathcal{M}_1$, $\mathcal{M}_2$, and $\mathcal{M}_3$. Correctly specified models exhibit lower DIC values $DIC(\mathcal{M}_1)=120851.9$, $DIC(\mathcal{M}_2)=23774.6$, $DIC(\mathcal{M}_3)=25953.6$ than those of the misspecified models $DIC(\mathcal{M}_1)=336530.6$, $DIC(\mathcal{M}_2)=345649.9$, $DIC(\mathcal{M}_3)=52027.5$. Failing to account for overdispersion leads to higher DIC values.

Applications

We illustrate the relevance of the proposed models through two empirical applications exhibiting overdispersion. In both cases, we estimate and compare the three model specifications, $\mathcal{M}_1$, $\mathcal{M}_2$, and $\mathcal{M}_3$, under both the Poisson and the Generalized Poisson (GP) likelihoods. Model $\mathcal{M}_2$ is estimated by fixing $k=2$, while Model $\mathcal{M}_3$ adopts a latent space representation with dimension $d=2$.

Citibike dataset

In the first application, we consider the Citibike dataset citibike2019, which contains information on rides between any two stations of the Citibike bike-sharing service in New York City. Part of this dataset -- one weekday in 2019 -- has been previously studied in he2025semiparametric where the authors proposed a semi-parametric dynamic Poisson LS network model. In this application, we consider bike-sharing activity for the entire year 2019 at a monthly frequency. We aggregate ride counts between any two stalls at the Neighborhood--Tabulation--Area (NTA) level. We thus obtain a count-weighted symmetric network with 61 nodes representing the 61 NTA areas and edges representing the bike-sharing service.

The top panel of Table (ref) reports the overdispersion parameter $\hat\rho$ for for model $\mathcal{M}_1$, $\mathcal{M}_2$, and $\mathcal{M}_3$ under the GP likelihood assumption. In all three model specifications, we observe strong overdispersion. Model $\mathcal{M}_3$ exhibits the lowest $\hat\rho$ as overdispersion is also captured by the latent coordinates. The bottom panel of Table (ref) reports the model-sections results in terms of DIC for $\mathcal{M}_1$, $\mathcal{M}_2$, and $\mathcal{M}_3$ both under Poisson and GP likelihood. Overall, the GP specification is preferred to its Poisson counterpart. As $\mathcal{M}_3$ is the preferred model under both likelihood assumptions, we present a description and comparative analysis of the results under this model. Further estimates are reported in Section (ref) of the Supplement.

The top-left panel in Figure (ref) shows the posterior mean of the $\alpha_i$ parameter associated with each of the 61 New York City neighborhoods. We observe that prominent and touristic NTAs, e.g., those located on the island of Manhattan (Murray Hill-Kips Bay, East Village, and Union Square), exhibit a higher $\alpha_i$ coherently with the interpretation of these parameters in terms of node centrality. The top-right panel shows the posterior predictive (PP) of the average expected strength over time. The model correctly captures the network's yearly seasonality. Higher riding activity is observed in the spring and summer, while a reduction in the activity is observed in the fall and winter.

In the bottom panels, we present a comparison of latent space estimation with $d=2$ under the Poisson and GP specifications for April 2019. Further estimates are reported in the Supplement (Section (ref)). We observe that the latent-space representation enables the retrieval of the relative geographical positions of the nodes (compare with the left panel). Stalls located in the Manhattan Island cluster together, and so do stalls located in the Bronx and Queens. Latent coordinates are more scattered under the Poisson specification since the Poisson model struggles to capture overdispersion, and a higher posterior variance of the latent coordinates partially mitigates this misspecification. This result is consistent with the misspecification bias for the latent coordinates discussed in Section (ref) (see also Figure (ref)).

Media network dataset

As a second application, we consider the time-varying media network dataset (schmidt2018polarization,casarin2025media). The dataset enables the construction of timely networks in which nodes are news outlets, and the weight of each edge is the count of unique commenters who interact with a pair of news outlets within the time window under consideration. In this application, we consider media networks aggregated monthly, totaling 24 networks spanning 2015 and 2016.

The top panel of Table (ref) reports the overdispersion parameter $\hat\rho$ for for model $\mathcal{M}_1$, $\mathcal{M}_2$, and $\mathcal{M}_3$ under the GP likelihood assumption. In all three model specifications and for all countries, we observe the presence of overdispersion. Model $\mathcal{M}_3$ exhibits the lowest $\hat\rho$ as overdispersion is also partially captured by the latent coordinates. The bottom panel reports the comparison between Poisson and GP assumptions in terms of DIC. The improvement in DIC when accounting for overdipersion is evident in all the specifications considered.

As the GP LS model ($\mathcal{M}_3$) is the best performing in terms of DIC, we will report and comment on the estimation results under this specification. For illustration purposes, we show in the main text the results for France in Figure (ref). The left panel shows a spatial representation of the news outlets (located at their headquarters), with point size and color proportional to $\alpha_i$. As expected, we find that national news outlets, often located in the most prominent cities of the four countries, are more central within the networks, whereas local news outlets tend to be more peripheral at the geographical level. We also observe differences in the geographical distributions of the most central news outlets. From a cross-country comparison (see Figures (ref)-(ref) in the Supplement), France and Spain are more centralized, with prominent news outlets located in their respective capital cities, whereas Germany and Italy exhibit a more heterogeneous geographical distribution of central news outlets. This may be due to very well-known historical reasons related to a monocentric vs policentric development of the four countries. See tezis2008european, newman2015reuters for a complete description of the different four media environments.

The top-right panel in Figure (ref) shows the latent-space representation of news outlets in January 2016. There are two forms of clustering. First, national news outlets with similar editorial lines are close in the latent space. Second, news outlets that are geographically closer tend to cluster. This effect is stronger for local news outlets. Similar evidence is found for the other countries in the dataset (see top-plots of Figures (ref)-(ref) in Supplement).

Regarding the coordinate dynamics, there is evidence of heterogeneous trending behavior across countries (see the latent factor estimates in the bottom panels of Figures (ref)-(ref) in Supplement). For both France and Germany, there is evidence of an upward trend. In Italy, the factor exhibits a pronounced upward trend at the beginning of the sample, peaking in July 2015, followed by a gradual decline that reaches its minimum around March 2016. Subsequently, the latent factor increases again, suggesting a sharpening in the overall connectivity. The time-varying latent factor for Spain is relatively stable with peaks in July 2015, December 2015, February and March 2016, and October 2016. Overall, the trend behavior confirms the results obtained in casarin2025media, which fitted a Poisson Markov switching LS model with a large number of hidden states.

We conclude our analysis with an out-of-sample prediction exercise. As LS models reveal themselves to be powerful for edge weight imputation, we inject 3 NAs at random in the last adjacency matrix of each country. We then test the ability of the estimated LS model to impute the missing values under the Poisson and GP specifications.

Table (ref) reports the prediction metrics by country under the two likelihood assumptions averaged over 10 simulations. Dealing with point-wise prediction MAE, MSE, RMSE), we notice that the results are mixed (GP outperforms Poisson for France and Germany, vice versa for Italy and Spain). Dealing with the distributional forecast, we notice that the GP model consistently shows a more satisfactory coverage ($>$ 90.0%). As concerns tail probabilities, the GP model appears much better calibrated (mean and variance of the tail probability, respectively, close to 1/2 and 1/12) at the cost of a higher interval width in all cases. To sum up, the Poisson model, despite providing competitive point-wise predictions, is also overconfident, while the GP clearly wins in terms of uncertainty quantification.

Conclusion

New model classes for count--weighted temporal networks are proposed. The Generalized Poisson (GP) distribution captures edge--weight underdispersion and overdispersion, while lagged network features or latent variables drive the dynamics of edge weights. Sufficient conditions for latent variable identification are derived, and a Bayesian inference framework is developed together with an efficient posterior sampling procedure. Numerical illustrations provide evidence of misspecification bias and significant posterior predictive errors when weight overdispersion is neglected. Bike-sharing and media temporal network applications demonstrated the GP latent space model's ability to accommodate unequal dispersion, as well as other features such as seasonality, trend, and spatial effects.

\setcounter{equation}{0} \setcounter{proposition}{0}