Multivariate count time series segmentation with “sums and shares” and Poisson lognormal mixture models: a comparative study using pedestrian flows within a multimodal transport hub

Abstract

This paper deals with a clustering approach based on mixture models to analyze multidimensional mobility count time-series data within a multimodal transport hub. These time series are very likely to evolve depending on various periods characterized by strikes, maintenance works, or health measures against the Covid19 pandemic. In addition, exogenous one-off factors, such as concerts and transport disruptions, can also impact mobility. Our approach flexibly detects time segments within which the very noisy count data is synthesized into regular spatio-temporal mobility profiles. At the upper level of the modeling, evolving mixing weights are designed to detect segments properly. At the lower level, segment-specific count regression models take into account correlations between series and overdispersion as well as the impact of exogenous factors. For this purpose, we set up and compare two promising strategies that can address this issue, namely the “sums and shares” and “Poisson log-normal” models. The proposed methodologies are applied to actual data collected within a multimodal transport hub in the Paris region. Ticketing logs and pedestrian counts provided by stereo cameras are considered here. Experiments are carried out to show the ability of the statistical models to highlight mobility patterns within the transport hub. One model is chosen based on its ability to detect the most continuous segments possible while fitting the count time series well. An in-depth analysis of the time segmentation, mobility patterns, and impact of exogenous factors obtained with the chosen model is finally performed.

Appendices

Appendix A: Mixture of sums and shares model estimation

Given $z_{j,s}=1$ and $\textbf{x}_{j,t}$, the series $\textbf{y}_{j,t}$ are distributed according to the following mixture model:

$$\begin{aligned} p(\textbf{y}_{j,t};{\varvec{{\varvec{\alpha }}}}, {\varvec{{\varvec{\gamma }}}}, r, {\varvec{{\varvec{\xi }}}}) = \sum _{s=1}^S \pi _s(j;{\varvec{{\varvec{\alpha }}}}) g(v_{j,t}|\textbf{x}_{j,t}, {\varvec{{\varvec{\gamma }}}}_s, r_s) h(\textbf{y}_{j,t}|v_{j,t},\textbf{x}_{j,t}, {\varvec{{\varvec{\xi }}}}_s), \end{aligned}$$

(A1)

with ${\varvec{{\varvec{\gamma }}}} = ({\varvec{{\varvec{\gamma }}}}_s)_{s=1,...,S}$, $r = (r_s)_{s=1,...,S}$ and ${\varvec{{\varvec{\xi }}}} = ({\varvec{{\varvec{\xi }}}}_s)_{s=1,...,S}$. The parameters of the model are estimated with the Expectation Maximization (EM) algorithm (Dempster et al. 1977) which requires a complete data log-likelihood maximization. The complete data log-likelihood can be written:

$$\begin{aligned} {\mathcal {L}}_c({\varvec{{\varvec{\alpha }}}}, {\varvec{{\varvec{\gamma }}}}, r, {\varvec{{\varvec{\xi }}}} ) = \sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T z_{j,s} \log \left( \pi _s(j;{\varvec{{\varvec{\alpha }}}}) g(v_{j,t}|\textbf{x}_{j,t}, {\varvec{{\varvec{\gamma }}}}_s, r_s) h(\textbf{y}_{j,t}|v_{j,t}, \textbf{x}_{j,t}, {\varvec{{\varvec{\xi }}}}_s) \right) . \end{aligned}$$

(A2)

Given the initial value of the parameters ${\varvec{{\varvec{\xi }}}}^{(0)}$, ${\varvec{{\varvec{\gamma }}}}^{(0)}$, $r^{(0)}$ and ${\varvec{{\varvec{\alpha }}}}^{(0)}$, the following two steps are repeated until convergence.

• Expectation step (E) The expectation of the completed log-likelihood is evaluated knowing the observed data Y and the set of current parameters: ${\varvec{{\varvec{\xi }}}}^{(c)}$, ${\varvec{{\varvec{\gamma }}}}^{(c)}$, $r^{(c)}$ and ${\varvec{{\varvec{\alpha }}}}^{(c)}$.

  $$\begin{aligned} Q({\varvec{{\varvec{\alpha }}}}^{(c)}, {\varvec{{\varvec{\gamma }}}}^{(c)}, r^{(c)}, {\varvec{{\varvec{\xi }}}}^{(c)}) =&\sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T E_{{\varvec{\xi }}^{(c)}, {\varvec{\gamma }}^{(c)}, r^{(c)}, {\varvec{\alpha }}^{(c)}}[z_{j,s}|Y] \end{aligned}$$

  (A3)
  $$\begin{aligned}&log\left( \pi _s(j;{\varvec{\alpha }}^{(c)}) g(v_{j,t}|\textbf{x}_{j,t}, {\varvec{\gamma }}^{(c)}_s, r^{(c)}_s) h(\textbf{y}_{j,t}|v_{j,t},\textbf{x}_{j,t}, {\varvec{\xi }}^{(c)}_s) \right) , \end{aligned}$$

  (A4)

  where

  $$\begin{aligned} E_{{\varvec{\xi }}^{(c)}, {\varvec{\gamma }}^{(c)}, r^{(c)}, {\varvec{\alpha }}^{(c)}}[z_{j,s}|Y]&= \tau ^{(c)}_{j,s} \end{aligned}$$

  (A5)
  $$\begin{aligned}&= \frac{\pi _s(j;{\varvec{\alpha }}^{(c)}) \prod _T g(v_{j,t}|{\varvec{x}}_{j,t}, {\varvec{\gamma }}^{(c)}_s, r^{(c)}_s) h(\textbf{y}_{j,t}|\textbf{x}_{j,t}, v_{j,t}, {\varvec{\xi }}^{(c)}_s)}{ \sum _{s'}\pi _{s'}(j;{\varvec{\alpha }}^{(c)}) \prod _T g(v_{j,t}|\textbf{x}_{j,t}, {\varvec{\gamma }}^{(c)}_s, r^{(c)}_s) h(\textbf{y}_{j,t}|\textbf{x}_{j,t}, v_{j,t}, {\varvec{\xi }}^{(c)}_s)}. \end{aligned}$$

  (A6)

  The a posteriori probabilities that each day j belongs to segment s, $\tau ^{(c)}_{j,s}$, are updated at each iteration of step E.

• Maximization step (M) Parameters ${\varvec{\xi }}^{(c+1)}$, ${\varvec{\gamma }}^{(c+1)}$, $r^{(c+1)}$ and ${\varvec{\alpha }}^{(c+1)}$ that maximize $Q({\varvec{{\varvec{\alpha }}}}^{(c)}, {\varvec{{\varvec{\gamma }}}}^{(c)},r^{(c)}, {\varvec{{\varvec{\xi }}}}^{(c)})$ are calculated. It is possible to rewrite this quantity as:

  $$\begin{aligned} Q({\varvec{{\varvec{\alpha }}}}^{(c)}, {\varvec{{\varvec{\gamma }}}}^{(c)}, r^{(c)}, {\varvec{{\varvec{\xi }}}}^{(c)}) = Q_1({\varvec{\alpha }}^{(c)}) + Q_2({\varvec{\gamma }}^{(c)},r^{(c)}) + Q_3({\varvec{\xi }}^{(c)}) \end{aligned}$$

  (A7)

  where

  $$\begin{aligned} Q_1({\varvec{\alpha }}^{(c)}) = \sum _{s=1}^S\sum _{j=1}^J \tau ^{(c)}_{j,s} log(\pi _s(j;{\varvec{\alpha }}^{(c)})) \end{aligned}$$

  (A8)
  $$\begin{aligned} Q_2({\varvec{\gamma }}^{(c)},r^{(c)}) = \sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T \tau ^{(c)}_{j,s} log(g(v_{j,t}|\textbf{x}_{j,t}, {\varvec{\gamma }}^{(c)}_s, r^{(c)}_s)) \end{aligned}$$

  (A9)
  $$\begin{aligned} Q_3({\varvec{\xi }}^{(c)}) = \sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T\tau ^{(c)}_{j,s} log(h(\textbf{y}_{j,t}|\textbf{x}_{j,t}, v_{j,t}, {\varvec{\xi }}^{(c)}_s)). \end{aligned}$$

  (A10)

  The maximisation of $Q_1$ consists in solving a weighted multinomial logistic regression. New values of ${\varvec{\alpha }}$ can be found using iterative procedures such as iteratively reweighted least squares (IRLS) (Holland and Welsch 1977). This problem is solved with with the multinom function of the nnet package (Ripley et al. 2016). $Q_2$ is the log-likelihood corresponding to a negative binomial generalized linear model. Its maximisation is solved through an alternating iteration process provided by the glm.nb function in the MASS package (Ripley et al. 2013). Within each segment s, for a given value of $r^{(c)}_s$ the linear model is fitted using an IRLS method. Next for fixed found ${\varvec{\gamma }}^{(c)}_s$ parameters, the $r^{(c)}_s$ parameter is estimated with score and information iterations. The two steps are alterned until convergence and ${\varvec{\gamma }}^{(c+1)}_s$ and $r^{(c+1)}_s$ are found. Note that $\tau ^{(c)}_{j,s}$ are here used as prior weights in the fitting process. The criterion $Q_3$, which is associated to a weighted Dirichlet multinomial regression model, is solved with the MGLM package (Kim et al. 2018). Because Dirichlet multinomial distribution does not belong to the exponential family, IRLS method is not used, as the expected information matrix is difficult to calculate. The method used here combine the minorization-maximization (MM) (Lange et al. 2000) algorithm and the Newton’s method. MM and Newton updates are computed at each iteration and the one with the higher log-likelihood is chosen.

Appendix B: Poisson log-normal mixture model estimation

The series $\textbf{y}_{j,t}$ are distributed according to the following mixture model:

$$\begin{aligned} p(\textbf{y}_{j,t};{\varvec{\alpha }}, {\varvec{\rho }}, {\varvec{\Sigma }}) = \sum _{s=1}^S \pi _s(j;{\varvec{\alpha }}) \int _{R^P}\left[ \prod _{p=1}^P g(\textbf{y}_{j,t,p}|\theta _{j,t,p})\right] m(\theta _{j,t}|{\varvec{\rho }}_s, {\varvec{\Sigma }}_s)d\theta _{j,t}, \end{aligned}$$

(B11)

with ${\varvec{\rho }} = ({\varvec{\rho }}_s)_{s=1,...,S}$ and ${\varvec{\Sigma }} = ({\varvec{\Sigma }}_s)_{s=1,...,S}$. g is a Poisson distribution and m a Gaussian distribution function. The EM algorithm can be used for parameter estimation but finding the expected value of the complete data log-likelihood requires estimating the conditional expectations $\mathop {\mathbb {E}}(Z_{js}\theta _{j,t}|\textbf{y}_{j,t},{\varvec{\rho }}_s,{\varvec{\Sigma }}_s)$ and $\mathop {\mathbb {E}}(Z_{js}\theta _{j,t}\theta '_{j,t}|\textbf{y}_{j,t},{\varvec{\rho }}_s,{\varvec{\Sigma }}_s)$ which are intractable. These conditional expectations can be calculated with an EM algorithm coupled with a Markov chain Monte Carlo (MCMC-EM) algorithm as presented in the work of Silva et al. (2019), which however comes with a heavy calculation load. We refer instead to the work presented by Chiquet et al. (2019) that uses variational approximation which is an approximate inference technique. The idea behind variational inference is to use Gaussian densities and approximate complex posterior distributions by minimizing the Kullback–Leibler divergence between the true $p(\theta )$ and approximating densities $q(\theta )$. The marginal log-likelihood for $\textbf{y}_{j,t}$ can be written as

$$\begin{aligned} \log p(\textbf{y}_{j,t}) = F(q(\theta _{j,t}), \textbf{y}_{j,t}) + D_{KL}(q(\theta _{j,t})|p(\theta _{j,t})), \end{aligned}$$

(B12)

with $D_{KL}(q(\theta _{j,t})|p(\theta _{j,t}))$ the Kullback-Leibler divergence between $p(\theta _{j,t})$ and $q(\theta _{j,t})$. $F(q(\theta _{j,t}), \textbf{y}_{j,t})$ is the expression of the variational lower bound of the log-likelihood. This is the criterion that we aim to maximize in the parameter estimation process. In the case of the Poisson-Lognormal model, q is assumed to be a Gaussian distribution:

$$\begin{aligned} q(\theta _{j,t}; \textbf{m}_{j,t}, \textbf{S}_{j,t}) = {\mathcal {N}}(\theta _{j,t}; \textbf{m}_{j,t}, \textbf{S}_{j,t}), \end{aligned}$$

(B13)

with $\textbf{m}_{j,t}$ and $\textbf{S}_{j,t} = diag(\textbf{S}_{j,t})$ the variational parameters associated with sample $\textbf{y}_{j,t}$ at day j and time slot t. To minimize the Kullback–Leibler divergence, the variational lower bound has to be maximized. The complete data log-likelihood can be written as follows:

$$\begin{aligned} \begin{aligned} {\mathcal {L}}_c({\varvec{\alpha }}, {\varvec{\rho }}, {\varvec{\Sigma }}, {\varvec{m}}, {\varvec{S}})&= \sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T z_{j,s} \log (\pi _s(j;{\varvec{\alpha }}))+\\ {}&\sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T z_{j,s}[F(q^{(s)}(\theta _{j,t}), \textbf{y}_{j,t}) + D_{KL}(q^{(s)}(\theta _{j,t})|p^{(s)}(\theta _{j,t}))], \end{aligned} \end{aligned}$$

(B14)

where $D_{KL}(q^{(s)}(\theta _{j,t})|p^{(s)}(\theta _{j,t}))$ is the Kullback–Leibler divergence between $p(\theta _{j,t}|\textbf{y}_{j,t}, z_{j}=s)$ and $q^{(s)}(\theta _{j,t})$ with $q^{(s)}(\theta _{j,t}) = {\mathcal {N}}(\textbf{m}_{j,t}^{(s)}, \textbf{S}_{j,t}^{(s)})$. And the variational lower bound of the log-likelihood for each observation $\textbf{y}_{j,t}$ is

$$\begin{aligned} \begin{aligned} & F(q^{(s)}(\theta _{j,t}), \textbf{y}_{j,t})= \frac{1}{2}\log |\textbf{S}_{j,t}^{(s)}|- \frac{1}{2}(\textbf{m}_{j,t}^{(s)} - \textbf{x}_{j,t}^T{\varvec{\rho }}_{s})'{\varvec{\Sigma }}_s^{-1}(\textbf{m}_{j,t}^{(s)} - \textbf{x}_{j,t}^T{\varvec{\rho }}_{s}) - tr({\varvec{\Sigma }}_s^{-1}\textbf{S}_{j,t}^{(s)}) - \\&\quad \frac{1}{2}\log |{\varvec{\Sigma }}_s|- \frac{P}{2} + (\textbf{m}^{(s)})^{'}_{j,t}\textbf{y}_{j,t} - \sum _{p=1}^P(\exp (m_{j,t,p}^{(s)} + \frac{1}{2}s_{j,t,p}^{(s)}) + \log (y_{j,t,p}!)). \end{aligned} \end{aligned}$$

(B15)

The EM algorithm is used to estimate the parameters and the following two steps are repeated until convergence.

• Expectation step (E) The expectation of the completed log-likelihood is evaluated knowing the observed data Y, the set of current parameters ${\varvec{\rho }}^{(c)}$, ${\varvec{\Sigma }}^{(c)}$ and ${\varvec{\alpha }}^{(c)}$ and variational parameters $\textbf{m}^{(c)}_{j,t}$, $\textbf{S}^{(c)}_{j,t}$.

  $$\begin{aligned} \begin{aligned} Q({\varvec{\rho }}^{(c)},{\varvec{\Sigma }}^{(c)},{\varvec{\alpha }}^{(c)},\textbf{m}^{(c)},\textbf{S}^{(c)})&= \sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T \tau _{j,s}^{(c)} \log (\pi _s(j;{\varvec{\alpha }}^{(c)})) +\\&\sum _{s=1}^S\sum _{j=1}^J\sum _{t=1}^T \tau _{j,s}^{(c)} E_{{\varvec{\rho }}^{(c)}, {\varvec{\Sigma }}^{(c)}, {\varvec{\alpha }}^{(c)}, m^{(c)}_{j,t}, S^{(c)}_{j,t}}[F(q^{(s)}(\theta _{j,t}), \textbf{y}_{j,t}) +\\&D_{KL}(q^{(s)}(\theta _{j,t})|p^{(s)}(\theta _{j,t}))], \end{aligned} \end{aligned}$$

  (B16)

  with $\tau _{j,s}^{(c)} = E_{{\varvec{\rho }}^{(c)}, {\varvec{\Sigma }}^{(c)}, {\varvec{\alpha }}^{(c)}, m^{(c)}_{j,t}, S^{(c)}_{j,t}}[z_{j,s}|Y]$. The variational lower bound of the log-likelihood is used to approximate $\tau _{j,s}^{(c)}$:

  $$\begin{aligned} \tau _{j,s}^{(c)} = \frac{\pi _s(j;{\varvec{\alpha }}^{(c)})\prod _{t=1}^T \exp (F(q^{(s)}(\theta _{j,t}), \textbf{y}_{j,t}))}{\sum _{h=1}^S\pi _h(j;{\varvec{\alpha }}^{(c)})\prod _{t=1}^T \exp (F(q^{(h)}(\theta _{j,t}), \textbf{y}_{j,t}))}. \end{aligned}$$

  (B17)

  Note that this approximation is used in the R package PLNmodels.

• Maximization step (M) The maximization step is divided into two parts:

  • Conditionally on ${\varvec{\rho }}_s$ and ${\varvec{\Sigma }}_s$ and given $\tau _{j,s}$, variational parameters $\textbf{m}^{(c)}_{j,t}$ and $\textbf{S}^{(c)}_{j,t}$ are updated. Because $F(q^{(s)}(\theta _{j,t}), \textbf{y}_{j,t})$ is strictly concave with respect to $\textbf{m}^{(c)}_{j,t}$ and $\textbf{S}^{(c)}_{j,t}$, it is possible to obtain $\textbf{S}^{(c+1)}_{j,t}$ with the fixed-point method and $\textbf{m}^{(c+1)}_{j,t}$ with Newton’s method.

  • Knowing $\tau _{j,s}^{(c)}$, $\textbf{m}^{(c+1)}_{j,t}$ and $\textbf{S}^{(c+1)}_{j,t}$ parameters ${\varvec{\rho }}^{(c+1)}$, ${\varvec{\Sigma }}^{(c+1)}$ and ${\varvec{\alpha }}^{(c+1)}$ are obtained.

Appendix C: Description of the time segments

See Table 5.

Table 5 Time segmentation