Seasonality in Dynamic Stochastic Block Models

TL;DR

This paper introduces a new statistical model for dynamic networks with seasonal time dependencies, enabling the recovery of underlying seasonal processes from data, which is crucial for understanding time-varying sociotechnological and geospatial systems.

Contribution

It proposes a novel dynamic stochastic block model incorporating seasonal effects and provides an inference method to recover these seasonal processes from data.

Findings

Successful recovery of latent seasonal processes in synthetic networks

Model captures time-varying edge formation driven by seasonal patterns

Demonstrates effectiveness in understanding dynamic network structures

Abstract

Sociotechnological and geospatial processes exhibit time varying structure that make insight discovery challenging. This paper proposes a new statistical model for such systems, modeled as dynamic networks, to address this challenge. It assumes that vertices fall into one of k types and that the probability of edge formation at a particular time depends on the types of the incident nodes and the current time. The time dependencies are driven by unique seasonal processes, which many systems exhibit (e.g., predictable spikes in geospatial or web traffic each day). The paper defines the model as a generative process and an inference procedure to recover the seasonal processes from data when they are unknown. Evaluation with synthetic dynamic networks show the recovery of the latent seasonal processes that drive its formation.

Figures

Equations

$$c_t^{(a,b)}=m_t^{(a,b)}+s_t^{(a,b)}$$

$$s_t^{(a,b)}=-\sum _{i=1}^{d-1}{s}_{t-i}^{(a,b)}+{\delta }_{s_t^{(a,b)}}$$

$$E_t^{(a,b)}=c_t^{(a,b)}+{ϵ}_{E_t^{(a,b)}}$$

$$[A_t^{(a,b)}{]}_{ij}\sim Bernoulli(E_t^{(a,b)})$$

$${\mathbf{x}}_t^{(a,b)}={\left[\begin{array}{ccccc}{m}_t^{(a,b)}& s_t^{(a,b)}& s_{t-1}^{(a,b)}& \dots & s_{t-d+2}^{(a,b)}\end{array}\right]}^T$$

$$\mathbf{G}=\left[\begin{array}{cccccc}1& 0& 0& \dots & 0& 0\\ 0& -1& -1& \dots & -1& -1\\ 0& 1& 0& \dots & 0& 0\\ 0& 0& 1& \dots & 0& 0\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& \dots & 1& 0\end{array}\right]$$

$$n^{(a,b)}=\{\begin{array}{ll}a=b& \frac{|a|(|a|-1)}{2}\\ a\ne b& |a||b|\end{array}$$

$$p_t^{(a,b)}=n^{(a,b)}{E}_t^{(a,b)}+{\omega }_{p_t^{(a,b)}}$$

$$\mathbf{h}=\left[\begin{array}{ccccc}{n}^{(a,b)}& n^{(a,b)}& 0& \dots & 0\end{array}\right]$$

Full text

Dynamic Networks, Kalman Filter, Structural Time Series, State Space Model

††copyright: rightsretained††doi: ††isbn: ††conference: 4th Workshop on Complex Methods for Data and Web Mining; August 2017; Leipzig, Germany††journalyear: 2017††journalyear: 2017††copyright: acmlicensed††conference: WI ’17; August 23-26, 2017; Leipzig, Germany††price: 15.00††doi: 10.1145/3106426.3109424††isbn: 978-1-4503-4951-2/17/08

1. Introduction & Motivation

Many complex systems exhibit regular, time dependent, seasonal patterns. For example, human movement patterns are driven by the time of day (Song et al., 2010), and vehicle traffic densities exhibit predictable increases at certain hours causing rush hours and decreases at night (Li et al., 2009). This same ‘seasonal’, time dependent effect occurs when monitoring network bandwidth usage (Yoo and Sim, 2016) or when counting the number of clicks per day on a web page(Gonçalves and Ramasco, 2008).

We look to bring the notion of seasonality to statistical network modeling with a new kind of dynamic stochastic block model (DSBM). A DSBM asserts that system components (nodes) are grouped into several types, and the probability of observing a component relation or interaction (edges) are determined by the types of the incident nodes and time. Different kinds of DSBM consider different assumptions about the network formation process (Xu and Hero III, 2013; Xu, 2015; Yang et al., 2011; Xing et al., 2010; Ho et al., 2011), but none consider seasonality. A conceptual overview of the model we propose is given in Figure 1. It fuses structural time series (plates on the right of Figure 1) with a generative network model. We call this a seasonal DSBM (SDSBM). The model is transformed to a state-space model for scalable fitting to data by Kalman Filters, and expectation-maximization is used for parameter learning.

2. Model Specification

We first specify the model of the seasonal processes controlling edge dynamics. We assume that time is discrete, with the current time $t$ representing a time period of some resolution. We also assume the node types are provided. For each pair of node types $a$ and $b$, we consider a structural time series with a bias $m_{t}^{(a,b)}$ establishing an anchor for values of the time series and a seasonal offset $s_{t}^{(a,b)}$ that shifts the bias by the current seasonality position. The process at time $t$, denoted $c^{(a,b)}_{t}$, is

[TABLE]

with bias $m^{(a,b)}_{t}$ described by $m^{(a,b)}_{t}=m^{(a,b)}_{t-1}+\delta_{m^{(a,b)}_{t}}$ where $\delta_{m^{(a,b)}_{t}}\sim\mathcal{N}(0,q^{(a,b)}_{m})$ models possible process noise. The set of seasonal offsets are stored in a vector $\mathbf{s}^{(a,b)}=(s^{(a,b)}_{1},s^{(a,b)}_{2},...,s^{(a,b)}_{d})$ having $d$ components. $d$ reflects either the length or the resolution of a seasonal process (e.g., $d=60$ to model per minute changes over a process that cycles per hour) and is assumed to be provided by the user of the model. The components are:

[TABLE]

where $\delta_{s^{(a,b)}_{t}}\sim\mathcal{N}(0,q^{(a,b)}_{s})$. This form enforces a zero-sum constraint to increase identifiability (Murphy, 2012). It should be emphasized that $q^{(a,b)}_{m}$ and $q^{(a,b)}_{s}$ control the noise of the underlying seasonal process, where low noise represents a stable seasonality that does not vary significantly between realized periods.

To model how the seasonal processes govern the shape of a dynamic network, we define a random variable $E^{(a,b)}_{t}\in[0,1]$ as the expected density of edges spanning node types $(a,b)$ at time $t$:

[TABLE]

where $\epsilon_{E^{(a,b)}_{t}}\sim\mathcal{N}(0,r^{(a,b)})$ models possible measurement noise. Next we define an adjacency matrix $A_{t}$, where $[A_{t}]_{ij}=1$ if there exists an edge between nodes $i$ and $j$ at time $t$ and $[A_{t}]_{ij}=0$ otherwise. Denote $A_{t}^{(a,b)}$ as the submatrix of $A_{t}$ only containing the rows and columns representing type $a$ and type $b$ nodes. Then $A_{t}^{(a,b)}$ is defined by the random variable:

[TABLE]

Repeating this process for all blocks $(a,b)$ and time steps $t$ will generate a desired dynamic network $\mathcal{A}=\{A_{1},A_{2},...,A_{t}\}$.

3. Model Fitting

We now describe an inference procedure to fit the model to an observed $\mathcal{A}$. Since seasonal processes are latent in data, the task is to estimate a posterior distribution on each $m^{(a,b)}_{t}$ and $s^{(a,b)}_{t}$ for each pair of node types $(a,b)$. Kalman filters (Kalman et al., 1960) are an appropriate tool for this task, but requires transforming the generative model into a state-space model (SSM). A SSM is a time series model with hidden and observed variables (Murphy, 2012); here we define $\mathbf{x}_{t}$ as hidden and $w_{t}$ as observed variables, respectively. We can notice the bias and seasonal offsets from before were hidden, while the adjacency matrix is observed, foreshadowing the structure to be defined. A SSM creates observations at time $t$ by two linear models: An observation model $w_{t}=\mathbf{h}\mathbf{x}_{t}+\epsilon_{t}$ and a transition model $\mathbf{x}_{t}=\mathbf{G}\mathbf{x}_{t-1}+\Delta_{t}$. Here, observations $w_{t}$ are generated by a transformation of the output (defined by $\mathbf{h}$) of the underlying transition model. The transition model describes transformations within a hidden state space where transitions from time $t-1$ to time $t$ are defined by the matrix $\mathbf{G}$. Observations and transitions are to be subject to time dependent random noise, which are modeled by Gaussian distributions $\epsilon_{t}\sim\mathcal{N}(0,R_{t})$ and $\Delta_{t}\sim\mathcal{N}(0,\mathbf{Q}_{t})$ With $R_{t}$ and $\mathbf{Q}_{t}$ controlling the amount of observation and transition noise, respectively. Assuming parameters $\theta_{t}=\{\mathbf{h},\mathbf{G},R_{t},\mathbf{Q}_{t}\}$ are known, a Kalman Filter can be used to derive the exact posterior $Pr(\mathbf{x}_{t}|\theta_{t},w_{1},w_{2},...,w_{t})$, i.e., the probability of the hidden state value at time $t$ given the observations up to and including time $t$ (Murphy, 2012).

Now we transform the model specification into a state space to define the transition model $\mathbf{x}_{t}=\mathbf{G}\mathbf{x}_{t-1}+\Delta_{t}$. As we are assuming edges of different vertex types $(a,b)$ are independent of each other, we will formulate the inference in terms of a pairing $(a,b)$. The full inference is completed by repeating the process for all pairs $(a,b)$. First we transform Equations 1 and 2 to define the hidden state variable $\mathbf{x}_{t}$ and state transition $\mathbf{G}$. The hidden state will be composed of the bias and vector of seasonal offsets as a $d\times 1$ seasonal state vector for a period of length $d$:

[TABLE]

Note that all the seasonal offsets from $\mathbf{s}^{(a,b)}$ are maintained in the state for a given $t$, with the $d$th seasonal offset implicitly defined based on the zero-sum constraint. Now to perform the state transition from time $t-1$ to time $t$ we define a $d\times d$ matrix $\mathbf{G}$:

[TABLE]

In $\mathbf{G}$, we see that multiplication of the first row of $\mathbf{G}$ by $\mathbf{x}^{(a,b)}_{t}$ yields Equation 1 without noise being added, as only the bias term is updated. Multiplication of the second row of $\mathbf{G}$ by $\mathbf{x}^{(a,b)}_{t}$ will update a single seasonal offset as shown in Equation 2. The remaining rows of $\mathbf{G}$ serve to permute the remaining seasonal offsets, such that each offset $s^{(a,b)}_{i}$ is updated after a full period of $d$ time steps. Each time step will update the most current seasonal offset, and shift the remaining offsets right one index in the state vector. Next we define the $d\times 1$ noise vector $\Delta_{t}=[\delta_{m^{(a,b)}_{t}},\delta_{s^{(a,b)}_{t}},0,\dots,0]^{T}$ where the first element is the bias noise in Equation 1, the second element is seasonal noise in Equation 2, and the remaining elements are all 0 as there is no additional noise for permuting the seasonal offsets. These noise values are sampled from a zero mean Gaussian with $d\times d$ covariance matrix $\mathbf{Q}=\texttt{diag}[q_{m}^{(a,b)},q_{s}^{(a,b)},0,\dots,0]$. Assuming the bias noise $\delta_{m^{(a,b)}_{t}}$ and seasonal offset noise $\delta_{s^{(a,b)}_{t}}$ are independent, the off-diagonal elements of $\mathbf{Q}$ are zero. The remaining elements are all zero, reflecting the lack of noise for the permutation operations. We can assume $\mathbf{Q}$ is stationary, and drop the dependence on $t$. This complete formulation of the transition model is not new, and has been completed by other researchers such as in (Davis and Keck, 2005).

Our next task is to transform Equations 3 and 4 into the observation model $w_{t}=\mathbf{h}\mathbf{x}_{t}+\epsilon_{t}$. To do this, we will need to define some additional variables, and take advantage of a result of the central limit theorem for a large number of vertices with types $(a,b)$ to create an approximate Gaussian transformation. First we need a count of the number of possible edges in block $(a,b)$, so if there are $|a|$ nodes of type $a$ and $|b|$ nodes of type $b$ then define:

[TABLE]

Also define the random variable $p_{t}^{(a,b)}\sim Binomial(E^{(a,b)}_{t},n^{(a,b)})$ as the number of formed edges in block $(a,b)$ at time $t$, where $E^{(a,b)}_{t}$ is the expected edge density as determined by Equation 3. $p_{t}^{(a,b)}$ is simply a more mathematically convenient way to define the edge generation process from Equation 4 and does not change the overall model. Having a binomial random variable, for large enough $n^{(a,b)}E_{t}^{(a,b)}$ we can apply the central limit theorem to approximate the distribution of $p_{t}^{(a,b)}$ as Gaussian:

[TABLE]

where $\omega_{p^{(a,b)}_{t}}\sim\mathcal{N}(0,n^{(a,b)}E_{t}^{(a,b)}(1-E_{t}^{(a,b)}))$ is observation noise that is time dependent on $E_{t}^{(a,b)}$. This represents variation in the repeated binary decision process of whether pairs of complex system components will interact. For example, in a geospatial context, a seasonal process may dictate that people travel from home to work at 8am, yet people individually decide the precise time they leave for work. In this example, locations of the geospace are system components (vertices), and movement between these locations are interactions (edges). In the SDSBM, each system component will have a type, such that each individual home is a vertex and all these homes can have the same type of ‘residence’. The sum of individual departure time decisions creates variation captured in $\omega_{p^{(a,b)}_{t}}$.

Now that we have successfully transformed our edge sampling procedure from Equation 4 to an approximately Gaussian formulation in Equation 8, we can return to defining observation model parameters $w_{t}$ and $\mathbf{h}$. We will define the observed variable $w^{(a,b)}_{t}$ as the number of formed edges $p_{t}^{(a,b)}$. To create $\mathbf{h}$, define a transformation which takes as input a seasonal state vector $\mathbf{x}_{t}^{(a,b)}$ and produces as output the number of formed edges $w^{(a,b)}_{t}$. By combining the operations of Equations 3 and 8 we define:

[TABLE]

Examining $w^{(a,b)}_{t}=\mathbf{h}\mathbf{x}^{(a,b)}_{t}$ closer, we can see this multiplication both sums the bias $m_{t}^{(a,b)}$ and first seasonal offset $s^{(a,b)}_{t}$ following Equation 3, and multiplies by $n^{(a,b)}$ to match the repeated Bernoulli trials from Equation 8. Finally, we model the variance of the measurement noise $\epsilon_{t}$ by summing measurement noise $\epsilon_{E_{t}^{(a,b)}}$ and observation noise $\omega_{p^{(a,b)}_{t}}$. This definition allowing for more flexible modeling of many complex systems. $\epsilon_{t}$ is sampled from a zero mean Gaussian distribution with time dependent variance $R^{(a,b)}_{t}=n^{(a,b)}E_{t}^{(a,b)}(1-E_{t}^{(a,b)})+(n^{(a,b)})^{2}r^{(a,b)}$.

We have now transformed the generative procedure to the suitable SSM to allow easy inference via the Kalman Filter. Given an initial Gaussian belief state $Pr(\mathbf{x}^{(a,b)}_{0})$ with mean $\mathbf{x}_{0}^{(a,b)}$ and variance $\Sigma_{0}^{(a,b)}$, all subsequent belief states will be Gaussian as well. The closed-form updates for the posterior distribution

$Pr(\mathbf{x}^{(a,b)}_{t}|\theta^{(a,b)}_{t},w^{(a,b)}_{1},w^{(a,b)}_{2},...,w^{(a,b)}_{t})$ are not defined in this paper due to space constraints but are available in (Murphy, 2012).

To estimate the unknown noise parameters of the SSM and Kalman Filter $\theta_{t}=\{R_{t},\mathbf{Q}\}$ from data, we derive an expectation-maximization routine that iteratively converges to locally optimal point estimates. The update equations for $\mathbf{Q}$ can be found in other research such as (Harvey and Peters, 1990). The updates for $R_{t}$ is a new formulation, which deviates from the conventional EM routine for Kalman Filters, due to our separation of noise parameter $R^{(a,b)}_{t}$ into a combination of $E_{t}^{(a,b)}$, $n^{(a,b)}$, and $r^{(a,b)}$. $n^{(a,b)}$ is assumed fixed and $E_{t}^{(a,b)}$ can be estimated using the prediction step of the Kalman Filter as described in (Xu and Hero III, 2013). To estimate $r^{(a,b)}$, we set up an optimization routine, which maximizes the log-likelihood of the complete joint distribution $Pr(\mathbf{x}^{(a,b)}_{1},...,\mathbf{x}^{(a,b)}_{t},w^{(a,b)}_{1},...,w^{(a,b)}_{t})$. This optimization is completed each iteration of EM, until a locally optimal estimate for both $\mathbf{Q}$ and $R_{t}$ is found.

4. Results

To demonstrate the fidelity of the SDSBM, we perform several experiments on a synthetic dynamic network $\mathcal{A}=\{A_{1},A_{2},...,A_{T}\}$ generated following the process defined in Section 2. The network is defined with $k=3$ vertex types, creating $6$ unique pairs $(a,b)$. There are approximately $32$ nodes of each type, creating the number of possibles edges for each $(a,b)$ as $n^{(a,b)}=1000$. All pairings $(a,b)$ have the same underlying true seasonality. We start with an initial bias $m^{(a,b)}_{0}$ and initial vector of seasonal offsets $\mathbf{s}^{(a,b)}=(s_{1}^{(a,b)},s_{2}^{(a,b)},...,s_{d}^{(a,b)})$. We set $m_{0}^{(a,b)}=0.5$, and the seasonal offsets as samples from a sine wave with amplitude $0.3$. This will create expected block densities $E_{t}^{(a,b)}\in[.2,.8]$, where noisy samples may cover the full range $[0,1]$. For example, in some of the experiments, we set the length of seasonality as $8$, creating $\mathbf{s}^{(a,b)}=(.2,.3,.2,0,-.2,-.3,-.2,0)$. To determine the noise level of the system, we set bias variance $q^{(a,b)}_{m}=10^{-8}$ and seasonal offset variance $q^{(a,b)}_{s}=10^{-8}$. These are set very low to reflect the expectation of near constant seasonality. If we set the number of possible edges of types $(a,b)$ as $n^{(a,b)}=1000$, we define a ‘medium’ amount of noise at $r^{(a,b)}=5.5*10^{-3}$, which will result in a standard deviation of the number of formed edges of $75$. Repeating the operations in Equations 1-4 $T$ times will form the desired synthetic dynamic network $\mathcal{A}$.

To infer the underlying seasonal processes, thus fitting the model to $\mathcal{A}$, the user needs to provide values for seven hyperparameters. We assume the user knows both the node labels and desired length of seasonality $d$. Initial guesses need to be provided for the remaining hyperparameters of initial state mean $\mathbf{x}_{0}^{(a,b)}$, state variance $\Sigma_{0}^{(a,b)}$, measurement variance $r^{(a,b)}$, bias variance $q_{m}^{(a,b)}$, and seasonal offset variance $q_{s}^{(a,b)}$. These can be defaulted to initial values such as $1$, without having a significant consequence on the fidelity performance. A more formal exploration of the sensitivity of the model to these hyperparameters is left as future work. We present the results of three experiments. First we simply infer the hidden seasonal state vectors $\mathbf{x}_{t}^{(x,y)}$ for a specific pairing $(x,y)$ and for all $t$, to qualitatively demonstrate the goodness of fit of the algorithm. A visual of the seasonal state vector transformed to expected edge counts is given in Figure 2. The underlying seasonality is estimated well, as the estimates shown as blue squares, closely follow the true seasonality of the solid red line. In the second experiment, we look to demonstrate the model’s ability to recover the seasonal state vector with differing numbers of observations as controlled by $T$, the length of the dynamic network. We start with $T=2d$ and increase the number of observed periods until $T=10d$. To assess performance we calculated the mean-squared error (MSE) between the true $\mathbf{x}_{t}^{(a,b)}$ and estimated $\hat{\mathbf{x}}_{t}^{(a,b)}$ seasonal vectors for all $t$. In Figure 3 a negative linear dependence between MSE and number of observations is demonstrated in all pairs $(a,b)$. The model becomes increasingly accurate for larger datasets. In the third experiment, we again evaluate the MSE between truth and estimated seasonality, but this time we vary the amount of measurement noise $r^{(a,b)}$. We start with a low $r^{(a,b)}=5*10^{-4}$ and increasing a half order of magnitude (e.g., $5*10^{-4}$,$1*10^{-3}$, $5*10^{-3}$, etc…) until a high level of noise at $5*10^{-2}$. In Figure 4 we see the exponential dependence on the noise level. The performance of this model on real data may thus be sensitive to this measurement noise value.

5. Conclusion

This paper proposed a new statistical model for dynamic networks, leveraging the benefits of structural time series and stochastic block models. The generative specification and inference procedure are defined. We demonstrate the capabilities of the model on a synthetic dataset, showing some properties of the model such as the negative linear dependence to the number of observed seasonal periods.

Acknowledgements.

This work is supported by industry and government partners at the National Science Foundation’s I/UCRC Center for Surveillance Research and the Air Force Research Laboratory.