Supervised Multiscale Dimension Reduction for Spatial Interaction Networks

TL;DR

This paper presents spinlets, a multiscale supervised dimension reduction method for spatial interaction network data, using an empirical Bayes approach with a tree-guided partitioning to improve interpretability and relevance to responses.

Contribution

It introduces a novel multiscale Bayesian dimension reduction technique for spatial interaction networks, incorporating a tree-structured prior and an inverse Poisson model.

Findings

Effective in producing compact, interpretable representations

Applied successfully to soccer analytics data

Enhances understanding of complex spatial networks

Abstract

We introduce a multiscale supervised dimension reduction method for SPatial Interaction Network (SPIN) data, which consist of a collection of spatially coordinated interactions. This type of predictor arises when the sampling unit of data is composed of a collection of primitive variables, each of them being essentially unique, so that it becomes necessary to group the variables in order to simplify the representation and enhance interpretability. In this paper, we introduce an empirical Bayes approach called spinlets, which first constructs a partitioning tree to guide the reduction over multiple spatial granularities, and then refines the representation of predictors according to the relevance to the response. We consider an inverse Poisson regression model and propose a new multiscale generalized double Pareto prior, which is induced via a tree-structured parameter expansion scheme. Our approach is motivated by an application in soccer analytics, in which we obtain compact vectorial representations and readily interpretable visualizations of the complex network objects, supervised by the response of interest.

Tables

Table 1. Table 1: Mean and standard deviation (subscripts) of MAEs.

	Lasso $\left(\alpha =1\right)$		Elastic-net $\left(\alpha =0.5\right)$		Ridge $\left(\alpha =0\right)$
Method	Training	Testing	Training	Testing	Training	Testing
rMETIS (h=7)	${0.851}_{0.046}$	${0.908}_{0.183}$	${0.845}_{0.058}$	${0.910}_{0.184}$	${0.756}_{0.063}$	${0.899}_{0.182}$
rMETIS (h=8)	${0.861}_{0.068}$	${0.924}_{0.183}$	${0.857}_{0.073}$	${0.934}_{0.192}$	${0.848}_{0.085}$	${0.921}_{0.177}$
rMETIS (h=9)	${0.712}_{0.112}$	${0.913}_{0.189}$	${0.702}_{0.120}$	${0.909}_{0.188}$	${0.580}_{0.114}$	${0.901}_{0.180}$
S3O	${0.699}_{0.128}$	${0.919}_{0.208}$	${0.694}_{0.116}$	${0.910}_{0.190}$	${0.476}_{0.031}$	${0.879}_{0.190}$
Spinlets	${0.574}_{0.078}$	${0.909}_{0.208}$	${0.564}_{0.071}$	${0.899}_{0.205}$	${0.611}_{0.057}$	${0.867}_{0.198}$

Equations

$$(X_{i,j}|Y_i=y_i)\sim \text{Poisson}({\lambda }_{i,j}),{\eta }_{i,j}(y_i)=\ln ({\lambda }_{i,j})={\alpha }_j+{\mu }_i+y_i{\beta }_j,$$

$$f_j(x_j|Y=y)$$

$$a_j({\eta }_j(y))$$

$$f(\mathbf{x}|y)=g({\mathbf{\beta }}^T\mathbf{x},y)h(\mathbf{x}),$$

$$R_{\mathbf{\beta }}({\mathbf{x}}_i)={\mathbf{\beta }}^T{\mathbf{x}}_i=\sum _{j=1}^m{x}_{i,j}{\beta }_j=\sum _{j=1}^m{x}_{i,j}\sum _{(s,\ell )\in A_j}{\gamma }_{s,\ell }=\sum _{s=0}^h\sum _{\ell =1}^{2^s}{z}_{i,s,\ell }{\gamma }_{s,\ell }={\mathbf{\gamma }}^T{\mathbf{z}}_i=R_{\mathbf{\gamma }}({\mathbf{x}}_i).$$

$$\mathbf{D}=\text{setcounter}MaxMatrixCols15\left[\begin{array}{ccccccccccccccc}1& 1& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0\\ 1& 1& 0& 1& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 1& 0& 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 1& 1& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 1& 0& 0& 1& 0& 0& 0& 0& 0& 1& 0& 0& 0\\ 1& 0& 1& 0& 0& 1& 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 1& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 1& 0\\ 1& 0& 1& 0& 0& 0& 1& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right].$$

$${\gamma }_{s,\ell }$$

$${\gamma }_{s,\ell }\sim \mathcal{N}(0,{\tau }_{s,\ell }),{\gamma }_{s^{\prime }+1,2\ell -1}-{\gamma }_{s^{\prime }+1,2\ell }\sim \mathcal{N}(0,{\varphi }_{s^{\prime }+1,2\ell -1,2\ell }).$$

$$\mathbf{\Lambda }(\mathbf{\tau },\mathbf{\varphi })=\mathrm{blockdiag}[1/{\tau }_{0,1};{\mathbf{\Omega }}_{0,1};{\mathbf{\Omega }}_{1,1};{\mathbf{\Omega }}_{1,2};\dots ;{\mathbf{\Omega }}_{h-1,1},\dots ,{\mathbf{\Omega }}_{h-1,m/2}],$$

$${\mathbf{\Omega }}_{s^{\prime },\ell }=\left[\begin{array}{cc}\frac{1}{{\tau }_{s^{\prime }+1,2\ell -1}}& 0\\ 0& \frac{1}{{\tau }_{s^{\prime }+1,2\ell }}\end{array}\right]+\frac{1}{{\varphi }_{s^{\prime }+1,2\ell -1,\ell }}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right],(s^{\prime },\ell )\in {\mathcal{I}}_{h-1}.$$

$$\ln p(\mathbf{\gamma })$$

$$\displaystyle+\sum_{(s^{\prime},\ell)\in\mathcal{I}_{h-1}}\bigg{[}-\ln{(2\xi_{2})}-(\alpha_{2}+1)\ln{\bigg{(}1+\frac{|{\gamma}_{s^{\prime}+1,2\ell-1}-{\gamma}_{s^{\prime}+1,2\ell}|}{\alpha_{2}\xi_{2}}}\bigg{)}\bigg{]},$$

$$x_{i,j}\sim \text{Poisson}({\mu }_{i,j}),{\mu }_{i,j}=t_i{e}^{{\eta }_{i,j}},{\eta }_{i,j}=a+b_i+c_j+y_i{\beta }_j,$$

$$a\sim \mathcal{N}(0,{\omega }_a),b_i\sim \mathcal{N}(0,{\omega }_b),c_j\sim \mathcal{N}(0,{\omega }_c),i=2,\dots ,n,j=2,\dots ,m,$$

$$\ln p(\mathbf{X}|\mathbf{y},\mathbf{\gamma },\mathbf{\omega })=\ln \int p(\mathbf{X}|\mathbf{y},\mathbf{\gamma },a,\mathbf{b},\mathbf{c})p(a|{\omega }_a)p(\mathbf{b}|{\omega }_b)p(\mathbf{c}|{\omega }_c)\text{d}a\text{d}\mathbf{b}\text{d}\mathbf{c},$$

$$\mathcal{F}(q,\mathbf{\gamma },\mathbf{\omega })=⟨\log P(\mathbf{X};\mathbf{\gamma },\mathbf{\Psi },\mathbf{U}|\mathbf{y},\mathbf{\omega }){⟩}_{q(\mathbf{\Psi },\mathbf{U})}+H[q(\mathbf{\Psi },\mathbf{U})],$$

$$q^{(t)}(\mathbf{\Psi }):=\underset{q(\mathbf{\Psi })}{\mathrm{argmax}}\mathcal{F}(q(\mathbf{\Psi }),q^{(t-1)}(\mathbf{U}),{\mathbf{\gamma }}^{(t-1)},{\mathbf{\omega }}^{(t-1)}).$$

$$\mathcal{F}(q^{(t)}(\mathbf{\Psi }),q^{(t)}(\mathbf{U}),{\mathbf{\gamma }}^{(t-1)},{\mathbf{\omega }}^{(t-1)})\ge \mathcal{F}(q^{(t)}(\mathbf{\Psi }),q^{(t-1)}(\mathbf{U}),{\mathbf{\gamma }}^{(t-1)},{\mathbf{\omega }}^{(t-1)}).$$

$$\mathcal{F}(q^{(t)}(\mathbf{\Psi }),q^{(t)}(\mathbf{U}),{\mathbf{\gamma }}^{(t)},{\mathbf{\omega }}^{(t)})\ge \mathcal{F}(q^{(t)}(\mathbf{\Psi }),q^{(t)}(\mathbf{U}),{\mathbf{\gamma }}^{(t-1)},{\mathbf{\omega }}^{(t-1)}).$$

$$\mathcal{F}(q,\mathbf{\gamma },\mathbf{\omega })=⟨{\ell }_1(\mathbf{\Psi };\mathbf{\gamma }){⟩}_{q(\mathbf{\Psi })}+⟨{\ell }_2(\mathbf{U};\mathbf{\gamma },\mathbf{\omega }){⟩}_{q(\mathbf{U})},$$

$$⟨{\ell }_1(\mathbf{\Psi };\mathbf{\gamma }){⟩}_{q(\mathbf{\Psi })}={\mathbb{E}}_{p(\mathbf{\Psi }|{\mathbf{\gamma }}^{(t)},{\mathbf{\omega }}^{(t)},\mathbf{X},\mathbf{y})}[{\ell }_1(\mathbf{\Psi },\mathbf{\gamma })],$$

$${\ell }_1(\mathbf{\Psi },\mathbf{\gamma })$$

$$\displaystyle+\sum_{(s^{\prime},\ell)\in\mathcal{I}_{h-1}}\bigg{[}\ln p(\delta_{s^{\prime}+1,2\ell-1,2\ell},\phi_{s^{\prime}+1,2\ell-1,2\ell},\nu_{s^{\prime}+1,2\ell-1,2\ell}|\alpha_{2},\eta_{2})\bigg{]}.$$

$$\ln p({\gamma }_{s,\ell },{\tau }_{s,\ell },{\lambda }_{s,\ell }|{\alpha }_1,{\eta }_1)=\ln p({\gamma }_{s,\ell }|{\tau }_{s,\ell })+\ln p({\tau }_{s,\ell }|{\lambda }_{s,\ell })+\ln p({\lambda }_{s,\ell }|{\alpha }_1,{\eta }_1),$$

$$\ln p({\delta }_{r,2\ell -1,2\ell },{\varphi }_{r,2\ell -1,2\ell },{\nu }_{r,2\ell -1,2\ell }|{\alpha }_2,{\eta }_2)$$

$$+\ln p({\nu }_{r,2\ell -1,2\ell }|{\alpha }_2,{\eta }_2).$$

$$p(\mathbf{\tau },\mathbf{\lambda }|-)=\prod _{l\in {\mathcal{I}}_{h-1}\bigcup \mathcal{L}}p({\tau }_l,{\lambda }_l|-),p({\tau }_l,{\lambda }_l|-)=p({\tau }_l|{\lambda }_l,-)p({\lambda }_l|-),$$

$$⟨{\ell }_1(\mathbf{\Psi };\mathbf{\gamma }){⟩}_{q(\mathbf{\Psi })}=⟨{\ell }_1(\mathbf{\Psi };\mathbf{\gamma }){⟩}_{p(\mathbf{\Psi }|{\mathbf{\gamma }}^{(t)},-)}=-\sum _{l=1}^L\frac{{\gamma }_l^2}{2}⟨{\tau }_l^{-1}⟩-\sum _{(u,w)\in {\mathcal{C}}_{s^{\prime },\ell }}\frac{{\delta }_{u,w}^2}{2}⟨{\varphi }_{u,w}^{-1}⟩.$$

$${\mathbb{E}}_{p({\rho }_l|{\lambda }_l,-)}[{\rho }_l]$$

$$⟨{\rho }_l⟩$$

$$⟨{\upsilon }_{u,w}⟩=⟨{\varphi }_{u,w}^{-1}⟩=\frac{({\alpha }_2+1)}{|{{\delta }_{u,w}}^{(t)}|[|{{\delta }_{u,w}}^{(t)}|+{\eta }_2]}.$$

$$\tilde{\mathbf{\Lambda }}=\mathrm{blockdiag}[⟨{\rho }_{0,1}⟩;{\tilde{\mathbf{\Omega }}}_{0,1};{\tilde{\mathbf{\Omega }}}_{1,1};{\tilde{\mathbf{\Omega }}}_{1,2};\dots ;{\tilde{\mathbf{\Omega }}}_{h-1,1},\dots ,{\tilde{\mathbf{\Omega }}}_{h-1,m/2}],$$

$${\tilde{\mathbf{\Omega }}}_{s^{\prime },\ell }=\left[\begin{array}{cc}⟨{\rho }_{s^{\prime }+1,2\ell -1}⟩& 0\\ 0& ⟨{\rho }_{s^{\prime }+1,2\ell }⟩\end{array}\right]+⟨{\upsilon }_{s^{\prime }+1,2\ell -1,\ell }⟩\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right],(s^{\prime },\ell )\in {\mathcal{I}}_{h-1}.$$

Full text

Keywords: dimension reduction, generalized linear mixed model, multiresolution method, object data analysis, parameter expansion, structured sparsity

1 Introduction

In modern applications, we frequently encounter complex object-type data, such as functions (Ramsay and Silverman, 2006), trees (Wang and Marron, 2007), shapes (Srivastava et al., 2011), and networks (Durante et al., 2017). In many instances, such data are collected repeatedly under different conditions, with an additional response variable of interest available for each replicate. This has motivated an increasingly rich literature on generalizing regression on vector predictors to settings involving more elaborate object-type predictors with special characteristics, such as functions (James, 2002), manifolds (Nilsson et al., 2007), tensors (Zhou et al., 2013), and undirected networks (Guha and Rodriguez, 2018).

Complex objects are often built from simpler parts. In this article, we consider supervised dimension reduction in which the sampling unit of data is in the form of composite objects (CO), composed of a collection of primitive objects (POs). Many data types can be seen as instances of the CO family, such as networks with spatially coordinated edges, pathology images of cells with a variety of shaped nuclei, or microbiome samples with species related on a phylogenetic tree. As shown in Figure 1, in the conventional regression formulation with a vector predictor, each replicate corresponds to a single point in $\mathbb{R}^{p}$. In regression with a composite object-valued predictor, the data could lie in a relatively lower dimensional $\mathbb{R}^{d}$ space, but each replicate corresponds to a collection of points. The component POs in CO-type data can be enormous and mostly distinctive from one another across replicates, presenting new challenges for data exploration, analysis, and visualization.

We are motivated by the 2018 FIFA World Cup data collected by StatsBomb (https://statsbomb.com/), which contain information on the association between coordinated interactions between players and team performance. For each replicate $i\in\{1,\ldots,n\}$, we observe SPatial Interaction Networks (SPIN) data $\mathcal{E}_{i}:=\{e_{k}:k=1,\ldots,q_{i}\}$, which contains a collection of $q_{i}$ completed passes. As illustrated in Figure 2, $\mathcal{E}_{i}$ is in the CO form, with every constituent pass $e_{k}$ viewed as a PO logged with the spatial locations of the passer and receiver. The whole dataset contains $50,159$ completed passes with $49,988$ unique locations of origin-destination in the $64$ matches. The associated response $y_{i}$ may refer to a team performance metric such as goals scored or conceded, the number of shots on target, or other situational game factors.

We consider regression modeling with $n$ observations of a scalar response ${y}_{i}$ and SPIN as a CO-valued predictor $\mathcal{E}_{i}$, $i=1,\ldots,n$. A primary challenge is to represent the CO data in a malleable form that facilitates multivariate analysis. One common practice is to divide $\mathcal{P}:=\bigcup_{i=1}^{n}\mathcal{E}_{i}$, the complete set of POs, into $p$ non-overlapping subsets $\mathcal{P}={\bigcup}_{j=1}^{m}\Pi_{j}$ (with $\Pi_{j}\bigcap\Pi_{j^{\prime}}=\emptyset$ for $j,j^{\prime}\in\{1,\ldots,p\}$ and $j\neq j^{\prime}$) through a predefined partitioning scheme $\boldsymbol{\Pi}=[\Pi_{1},\ldots,\Pi_{p}]$. Under the partition $\boldsymbol{\Pi}$, the $\mathcal{E}_{i}$ can be represented as a $p$-dimensional count vector $\boldsymbol{x}_{i}\in\mathbb{Z}^{p}$ through a bag-of-words representation (Blei et al., 2003; Taddy, 2013), where $x_{i,j}=N(\mathcal{E}_{i}\bigcap{\Pi}_{j})$ counts the occurrences of POs appearing in $\mathcal{E}_{i}$ and belonging to subset $\Pi_{j}$.

For example, Miller et al. (2014) discretizes the basketball court uniformly into tiles and counts the number of shots located in each tile. Durante et al. (2017) parcellates the brain into regions and investigates the network connectivity between pairs of regions. In Figure 3, we plot the resulting tile-based network representation based on a nodal partitioning scheme. Unfortunately, such a non-data adaptive scheme throws away many relevant details and produces an overly sparse representation of the data.

Alternatively, one can directly divide POs into disjoint groups according to their similarity. This partition scheme inherently assumes the equivalence of the POs falling within the same group, and focuses on the variabilities in abundance across groups. The choice of partitioning scheme and its scale will have a critical influence on inference. Ideally, the induced vector representation should promote the interpretability of the CO-type predictor and preserve the relevance to the response. However, such approaches are underdeveloped in the current literature.

In this article, we propose spinlets—a supervised multiscale dimension reduction method for spatial interaction networks and CO-type data more broadly. Spinlets takes into account both the predictor structure and the label supervision: it uses the similarity among POs for deriving a bag-of-POs representation accompanied by a partitioning tree, and then refine the tree structure according to the relevance to the response. Spinlets serves as a tool for both supervised representation learning and visualization of complex SPIN data. We first review the relevant literature and then describe our spinlets method.

1.1 Relevant Literature

There is a rich literature on supervised dimension reduction, covering LASSO (Tibshirani, 1996), supervised PCA (Bair et al., 2006; Barshan et al., 2011), and sufficient dimension reduction (SDR), see Cook (2007), Adragni and Cook (2009) and the references cited therein. The reduction of complexity is typically achieved through variable selection or combination. Such approaches can accommodate vector predictors and perform well in high-dimensional settings; however, our application involves predictors with complex structures. SDR methods have been generalized to handle functional predictors (Ferré and Yao, 2003, 2005), matrix- or array-valued predictors (Li et al., 2010), and irregularly measured longitudinal predictors (Jiang et al., 2014). In this article, we center our focus on spatial networks as an instance of a composite data object.

There is a separate literature on multiscale geometric data representation, including diffusion maps (Lafon and Lee, 2006) and GMRA (Allard et al., 2012; Petralia et al., 2013). These approaches seek a reductive representation that reflects the intrinsic geometry of the high-dimensional data by partitioning the similarity graph of $n$ data observations. In contrast, our spinlets method partitions the similarity graph of $q=\sum_{i=1}^{n}q_{i}$ variables, with a different goal of identifying predictive groups of variables. Spinlets is similar in spirit with the treelets method (Lee et al., 2008), which organizes variables on a hierarchical cluster tree with multiscale bases; however, treelets is an unsupervised approach utilizing the sample correlation to construct the tree with a single cutoff height. Our spinlets approach departs from treelets by incorporating external proximity to construct the tree, and determining non-uniform heights (Meinshausen and Bühlmann, 2008) with reference to the response.

In regression with CO-type predictors, the total number of unique POs is massive, while only a limited number of them are sparsely observed within each replicate. It is advantageous to form groups of POs that are spatially contiguous, such that meaningful analysis can be conducted at a lower level of resolution. In many other applications, predictors are highly correlated, or collectively associated with the response, or domain knowledge exists suggesting the functional similarity among a group of variables. This has motivated a line of research on supervised clustering of predictors in forward regression settings. Examples include Hastie et al. (2001); Jörnsten and Yu (2003) and Dettling and Bühlmann (2004). The averaging operator on the predictors often leads to lower variance (Park et al., 2006).

Regularization methods such as elastic net (Zou and Hastie, 2005) or OSCAR (Bondell and Reich, 2008) can mitigate the multicollinearity issue and encourage grouping effects. Along this thread, Wang and Zhao (2017) recently proposed two tree-guide fused lasso (TFL) penalties, which effectively encode the topology of a phylogenetic tree in determining the taxonomic levels of microbes associated with a given phenotype. However, this approach does not model the variability in the predictors, while we model the conditional distributions of the predictors given the response through inverse regression, with possibilities of alleviating the effects of collinearity (Cook and Li, 2009). Moreover, Lasso-based penalties tend to over-shrink signals not close to zero (Armagan et al., 2011). We introduce a new multiscale prior that induces a locally adaptive shrinkage rule on each scale. This prior executes two special operators on the tree within our proposed spinlets method.

In Section 2, we introduce the Poisson inverse regression model and the reductive operators. Section 3 presents a tree-structured PX scheme and our new multiscale shrinkage prior. A variational expectation-maximization (EM) algorithm for estimation is outlined in Section 4. In Section 5, we evaluate the performance of our approach with simulated data and demonstrate the practical utility through applications to soccer analytics. The implementation of spinlets will be made available on Github.

2 Tree-Guided Supervised Dimension Reduction

As a preprocessing step, we build a partition tree via recursively applying the METIS partitioning algorithm (Karypis and Kumar, 1998) on the similarity graph $\mathcal{G}$ of POs, detailed in Appendix A. The multiscale proximity information of POs is encoded in a binary tree $\mathcal{T}_{h}$ of height $h$, which partitions $\mathcal{P}$ into $m=2^{h}$ groups on the finest scale and yields a primary vectorial representation $\boldsymbol{X}\in\mathbb{Z}^{m}$ of COs, aligned across replicates.

According to Cook (2007), a sufficient reduction of the random vector $\boldsymbol{X}$, denoted as $R(\boldsymbol{X})$, satisfies one of the three equivalent statements: $(i)$ $\boldsymbol{X}|(Y,R(\boldsymbol{X}))\sim\boldsymbol{X}|R(\boldsymbol{X})$; $(ii)$ $Y|\boldsymbol{X}\sim Y|R(\boldsymbol{X})$; $(iii)$ \boldsymbol{X}\rotatebox[origin={c}]{90.0}{\models}Y|R(\boldsymbol{X}), where $\sim$ indicates equivalence in distribution and

$\models$

denotes independence. Our main goal is to determine a reductive rule $R(\boldsymbol{X})$ by pruning the tree $\mathcal{T}_{h}$, such that the resulting representation is interpretable and retains the relevance to the response variable $Y$.

2.1 Poisson Inverse Regression Model

For replicate $i$ with response variable $Y_{i}$, we attach a random variable $X_{i,j}\in\mathbb{Z}$ to each leaf node $j$, counting the number of POs appearing in CO $\mathcal{E}_{i}$ that fall in the $j$-th leaf group. In order to explicitly model the variabilities in occurrences, we adopt an inverse regression formulation. This approach is motivated by poor performance we observed in implementing usual regression due to extreme multicollinearity issues. Sufficiency is guaranteed within our proposed Poisson inverse regression (PIR) model for $X_{i,j}$ conditionally on $Y_{i}$, $i=1,\ldots,n$, $j=1,\ldots,m$,

[TABLE]

where $\alpha_{j}$ is the intercept for predictor $j$, $\mu_{i}$ is the baseline effect for replicate $i$, and ${\beta}_{j}$ is the regression coefficient for predictor $j$. The linear sufficient reduction $R_{\boldsymbol{\beta}}(\boldsymbol{X})$ parameterized by $\boldsymbol{\beta}$ is derived as follows (the replicate index $i$ is omitted for clarity).

Proposition 1

Letting $R_{\boldsymbol{\beta}}(\boldsymbol{X})=\boldsymbol{\beta}^{T}\boldsymbol{X}$, under the inverse Poisson regression model (1), the distribution of $Y|\boldsymbol{X}$ is the same as the distribution of $Y|R_{\boldsymbol{\beta}}(\boldsymbol{X})$ for all values of $\boldsymbol{X}$.

Cook (2007) proves the sufficiency of $R_{\boldsymbol{\beta}}(\boldsymbol{X})=\boldsymbol{\beta}^{T}\boldsymbol{X}$ for one-parameter exponential family models. Within this family, the PIR model (1) can be written in the following form,

[TABLE]

Accordingly, the joint probability mass function $f(\boldsymbol{x}|y)$ of $\boldsymbol{X}|(Y=y)$ can be written as

[TABLE]

where $g(\boldsymbol{\beta}^{T}\boldsymbol{x},y)=\exp{(y\boldsymbol{\beta}^{T}\boldsymbol{x})}\prod_{j=1}^{m}a_{j}(\eta_{j}(y))$ and $h(\boldsymbol{x})=\prod_{j=1}^{m}[b_{j}(x_{j})\exp{({x}_{j}{\alpha}_{j})}]$. Thus, the sufficiency of linear reduction holds according to the Fisher-Neyman factorization theorem for sufficient statistics (Bickel and Doksum, 2015, Theorem 1.5.1, p. 43). The PIR model has a close connection with the multinomial inverse regression (MNIR) model (Taddy, 2013), though, the vector Poisson likelihood departs from multinomial likelihood by accounting for the variability of total number of POs in each replicate.

2.2 Reductive Operators: Deletion and Fusion

One reductive operator enabled by the $\boldsymbol{\beta}$ parameterization from the PIR model (1) is the deletion of irrelevant leaf groups. One can see that $\beta_{j}=0$ implies that $f(x_{j}|Y=y)\equiv f(x_{j})$, that is, the number of POs in the $j$-th leaf group is independent of the value of $Y$. Another reductive operator on the tree is fusion. We observe that if $\beta_{j}=\beta_{j^{\prime}}\neq 0$, $\forall j,j^{\prime}\in\mathcal{D}$, then $R_{\boldsymbol{\beta}}(\boldsymbol{x})$ is a function depending on the predictors only through $\sum_{j\in\mathcal{D}}{X}_{j}$, which is the total number of POs falling into the set $\mathcal{D}$. Therefore, the practitioner can construct a lower resolution vectorial representation by merging the involved leaf sets into one set $\mathcal{D}$ without loss of relevant information. The relevant signals are then captured on coarser scales. To ensure spatial contiguity, we require all leaves contained in $\mathcal{D}$ to share at least one common ancestor node.

The grouping of highly correlated predictors in high-dimensional regression can be incorporated via fusion penalties (Tibshirani et al., 2005; Bondell and Reich, 2008), which encourage sparsity in the differences of coefficients. There exist several generalized fusion schemes that can take into account graph structure (She, 2010; Tibshirani and Taylor, 2011). However, these methods do not support multiscale nested grouping of predictors in accordance with the tree structure. For example, applying the pairwise fused lasso (She, 2010) to all pairs of variables tends to incorrectly encourage merging all the variables together with equal strengths. The TFL penalties (Wang and Zhao, 2017) require careful tuning of the regularization parameters across multiple scales. There exists no shrinkage prior that supports the prior belief of sparsity and piece-wise smoothness across scales and locations along the tree.

3 Multiscale Shrinkage with Parameter Expansion

Parameter expansion (PX) (Liu et al., 1998) has been found useful not only for accelerating computations (Liu and Wu, 1999), but also for inducing new families of prior distributions (Gelman, 2004). In this section, we propose a new tree-structured PX scheme, which induces a multiscale shrinkage prior on $\boldsymbol{\beta}$.

3.1 Tree Structured Parameter Expansion

For $j=1,\ldots,m$, we denote by $\mathcal{A}_{j}$ the “vertical” root-to-leaf path set in the tree $\mathcal{T}_{h}$ connected to the $j$-th leaf group, $j=1,\ldots,m$. The path set $\mathcal{A}_{j}$ includes all the nodes that it visits from the root node to the leaf node $j$. Meanwhile, we denote the “horizontal” descendant set $\mathcal{D}_{s,\ell}$ as the set of leaf nodes who share a most recent common ancestor $(s,\ell)$, $(s,\ell)\in\mathcal{I}_{h-1}$, where $\mathcal{I}_{h-1}:=\{(s,\ell):0\leq s\leq h-1,1\leq\ell\leq 2^{s}\}$ is the internal nodes in $\mathcal{T}_{h}$, in which each parent node $(s,\ell)$ has two child nodes $(s+1,2\ell-2+t)$, $t\in\{1,2\}$. Specifically, for $s=h$, we set $\mathcal{D}_{h,j}\equiv\mathcal{L}_{j}$, for $j=1,\ldots,m$, for ease of notation. Figure 4 illustrates an example of a path set $\mathcal{A}_{8}$ and an example of a descendant set $\mathcal{D}_{2,4}$ with the common ancestor node $(2,4)$.

There exists a dual relationship between random variables and coefficients based on these two notations. We can attach a random variable $Z_{i,s,\ell}\in\mathbb{Z}$ to each node $(s,\ell)$, where $Z_{i,s,\ell}=N(\mathcal{E}_{i}\bigcap\mathcal{D}_{s,\ell})$ counts the appearances of POs in descendant set $\mathcal{D}_{s,\ell}$ and ${Z}_{i,s,\ell}=\sum_{j\in\mathcal{D}_{s,\ell}}{X}_{i,j}$. Letting $\gamma_{s,\ell}$ be the coefficient for node $(s,\ell)$ visited by the path $\mathcal{A}_{j}$, we have ${\beta}_{j}=\sum_{(s,\ell)\in{A_{j}}}{\gamma}_{s,\ell}$. Therefore, the sufficient reduction score $R_{\boldsymbol{\beta}}(\boldsymbol{x}_{i})$ introduced in Section 2 can be re-expressed under the parameterization $\boldsymbol{\gamma}$ on the partition tree $\mathcal{T}_{h}$,

[TABLE]

Clearly, $R_{\boldsymbol{\gamma}}(\boldsymbol{x}_{i})$ is also a linear sufficient reduction. The reparameterization changes neither the data likelihood of the PIR model (1), nor the sufficient reduction score.

The reparameterization by the tree-structured PX scheme can be represented in matrix form as $\boldsymbol{\beta}=\boldsymbol{D}\boldsymbol{\gamma}$, where $\boldsymbol{D}$ is a $m\times L$ design matrix with binary entries, $m<L$. Each column in $\boldsymbol{D}$ can be interpreted as a basis function that encodes the piecewise smoothness at a different location and scale. For example, assuming $h=3$, the number of leaves $m=2^{3}=8$, $L=15$, we have

[TABLE]

3.2 Fused Generalized Double Pareto Prior

We denote by $\mathcal{C}_{s,\ell}$ the set of child nodes of $(s,\ell)$ and the set $\mathcal{F}_{s,\ell}$ all the nodes on the sub-branch rooted from node $(s,\ell)$. The tree $\mathcal{T}_{h}$ is originated from the root node $(0,1)$ with leaf nodes indexed by $\mathcal{L}:=\{(h,\ell):1\leq\ell\leq 2^{h}\}$, $\mathcal{F}_{0,1}=\mathcal{I}_{h-1}\bigcup\mathcal{L}$. In Section 2.2, we introduced two reductive operations: deletion (if $\beta_{j}=0$) and fusion (if $\beta_{j}=\beta_{j^{\prime}}$, $\forall j,j^{\prime}\in\mathcal{D}_{s,\ell}$) on the leaf partitions along the tree. However, exhaustive search for all possible schemes is prohibitive. Even for a binary tree $\mathcal{T}_{4}$, these two operations in combination result in $458,330$ different schemes. Alternatively, effective execution of the following two operations can be induced by regularization in the $\boldsymbol{\gamma}$ parameterization:

1. (i)

Deletion: if ${\gamma}_{s,\ell}=0$, $\forall(s,\ell)\in\mathcal{A}_{j}$, then $\beta_{j}=\sum_{(s,\ell)\in\mathcal{A}_{j}}{\gamma}_{s,\ell}=0$, which implies that the contributions of leaf predictor $j$ across all the scales are pruned out. 2. (ii)

Fusion: If $\forall(s^{\prime},\ell)\in\mathcal{F}_{s,\ell}$, their child nodes satisfy ${\gamma}_{s^{\prime}+1,2\ell-1}={\gamma}_{s^{\prime}+1,2\ell}$, then $\beta_{j}=\beta_{j^{\prime}}$, $\forall j,j^{\prime}\in\mathcal{D}_{s,\ell}$, the leaf variables within $\mathcal{D}_{s,\ell}$ can be condensed into one variable.

Note that the $\boldsymbol{\gamma}$ parameterization is redundant; both conditions above are sufficient but not necessary. Based on the above observations, we impose generalized double Pareto (GDP) priors (Armagan et al., 2013) on $\boldsymbol{\gamma}$ and the pairwise differences between sibling nodes,

[TABLE]

where $(s^{\prime}+1,2\ell-1),(s^{\prime}+1,2\ell)\in\mathcal{C}_{s^{\prime},\ell}$, $(s^{\prime},\ell)\in\mathcal{I}_{h-1}$ and $(s,\ell)\in\mathcal{F}_{0,1}$. The first prior encourages sparsity on the individual coefficients and the second prior promotes sparsity on the differences between pairs of siblings with a common parent node $(s^{\prime},\ell)$. These priors lead to a generalized fused lasso-type penalty (She, 2010; Tibshirani and Taylor, 2011), but the GDP prior corresponds to a reweighted $\ell_{1}$ penalty instead of $\ell_{1}$ (as will be seen in (6)), which better approximates the $\ell_{0}$-like criterion (Candes et al., 2008).

For $\mathcal{T}_{h}$, the number of expanded parameters in $\boldsymbol{\gamma}$ is $L=2^{h+1}-1$. A natural question is whether there exists a multivariate prior on $\boldsymbol{\gamma}$ that could justify the compatibility of these two GDP priors. To see this, we use the latent variable representation of the GDP prior introduced in Armagan et al. (2013). The first level of the hierarchy is written as,

[TABLE]

We postulate a multivariate normal prior for the $L$-dimensional vector $\boldsymbol{\gamma}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Lambda}^{-1})$ with $L\times L$ precision matrix $\boldsymbol{\Lambda}$, whose log marginal density is different than that of priors in (3) only up to a constant. The entries in $\boldsymbol{\Lambda}$ as a function of $(\boldsymbol{\tau},\boldsymbol{\phi})$ can be found by square completing. It takes a block-diagonal form as follows,

[TABLE]

where

[TABLE]

To complete the hierarchy of the multivariate GDP prior with $\boldsymbol{\gamma}\sim\mathcal{N}(\boldsymbol{0},[\boldsymbol{\Lambda}(\boldsymbol{\tau},\boldsymbol{\phi})]^{-1})$, we put $\tau_{s,\ell}\sim\mathrm{Exp}(\lambda_{s,\ell}^{2}/2)$, $\lambda_{s,\ell}\sim\mathrm{Ga}(\alpha_{1},\eta_{1})$, $(s,\ell)\in\mathcal{F}_{0,1}$, and $\phi_{s^{\prime}+1,2\ell-1,2\ell}\sim\mathrm{Exp}(\nu_{s^{\prime}+1,2\ell-1,2\ell}^{2}/2)$, $\nu_{s^{\prime}+1,2\ell-1,2\ell}\sim\mathrm{Ga}(\alpha_{2},\eta_{2})$, $(s^{\prime},\ell)\in\mathcal{I}_{h-1}$. So now we have a fused generalized double Pareto (fGDP) prior, which promotes the desired form of structured sparsity in $\boldsymbol{\gamma}$, and enjoys a latent variable representation that makes the parameter estimation straightforward. Integrating out the latent variables $\boldsymbol{\Psi}$, we obtain the marginal density of the fGDP prior, denoted by $\textrm{fGDP}(\boldsymbol{\gamma};\mathcal{T}_{h},\alpha_{1},\eta_{1},\alpha_{2},\eta_{2})$, whose logarithm takes the following form,

[TABLE]

where $\xi_{1}=\eta_{1}/\alpha_{1}$, $\xi_{2}=\eta_{2}/\alpha_{2}$.

Importantly, through PX, we have transformed the problem of multiscale shrinkage on the regression coefficients $\boldsymbol{\beta}$ across multiple scales on $\mathcal{T}_{h}$ into a structured shrinkage problem on the expanded parameters $\boldsymbol{\gamma}$, which can be conveniently addressed via the proposed fGDP prior. The hierarchical-Bayes representation of the multiscale shrinkage prior on $\boldsymbol{\beta}$ can be obtained via integrating out $\boldsymbol{\gamma}$; we have the conditional prior $\boldsymbol{\beta}|\boldsymbol{\tau},\boldsymbol{\phi}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{D}\boldsymbol{\Lambda}(\boldsymbol{\tau},\boldsymbol{\phi})^{-1}\boldsymbol{D}^{T})$ and the priors on the latent variables $\{\boldsymbol{\tau},\boldsymbol{\phi},\boldsymbol{\lambda},\boldsymbol{\nu}\}$ do not change. However, the precision matrix $\boldsymbol{D}\boldsymbol{\Lambda}(\boldsymbol{\tau},\boldsymbol{\phi})^{-1}\boldsymbol{D}^{T}$ no longer exhibits a sparse block-diagonal structure as in (4), and the resulting EM procedure of estimating $\boldsymbol{\beta}$ is less tractable than estimating $\boldsymbol{\gamma}$ as the former involves intractable expectations.

4 Parameter Estimation

We further assume each replicate is collected within a time window of length $t_{i}$ (known), $i=1,\ldots,n$. To accommodate potential overdispersion and dependencies, we incorporate random effects in the model. The Poisson log-linear mixed regression model is written as follows,

[TABLE]

where $(\beta_{1},\ldots,\beta_{m})$ is the fixed effect slope parameter for the $m$ simple Poisson mixed regression model (Hall et al., 2011b). The fixed effects measure the common association between the predictors and response, while the random effects allow replicates or leaf groups to have their own baseline rates. The total, column and row random effects are $a$, $\boldsymbol{b}$, and $\boldsymbol{c}$, respectively. Constraints are needed for the identifiability of row and column scores $b_{i}$ and $c_{j}$, so we use the corner constraint (Yee and Hastie, 2003) $b_{1}\equiv c_{1}\equiv 0$ in this article. Gaussian priors on $a$, $\boldsymbol{b}=[b_{2},\ldots,b_{n}]$ and $\boldsymbol{c}=[c_{2},\ldots,c_{m}]$ are specified as follows,

[TABLE]

with unknown variance parameters $\boldsymbol{\omega}=[\omega_{a},\omega_{b},\omega_{c}]$. Since $\beta_{j}=\boldsymbol{d}_{j}^{T}\boldsymbol{\gamma}$, we have $\eta_{i,j}=a+b_{i}+c_{j}+y_{i}\boldsymbol{d}_{j}^{T}\boldsymbol{\gamma}$. Note that the sufficiency of $R_{\boldsymbol{\gamma}}(\boldsymbol{x}_{i})$ established in Section 2.1 still holds conditional on the random effect terms (Taddy, 2013). In the next section, we introduce a penalized likelihood estimator of $\boldsymbol{\gamma}$, under the conditional likelihood ${\ell}(\boldsymbol{\gamma},\boldsymbol{\omega}):=\ln{p(\boldsymbol{X}|\boldsymbol{y},\boldsymbol{\gamma},\boldsymbol{\omega})}$ with the proposed fGDP prior $\boldsymbol{\gamma}\sim\textrm{fGDP}(\alpha_{1},\eta_{1},\alpha_{2},\eta_{2})$ guided by $\mathcal{T}_{h}$.

4.1 Variational Expectation Maximization (EM)

The hierarchical-Bayes representation of the fGDP prior facilitates an iterative EM-type algorithm for penalized estimation with (6). We adopt the type-I estimation framework (Figueiredo, 2003), which treats $\boldsymbol{\Psi}:=\{\boldsymbol{\tau},\boldsymbol{\phi},\boldsymbol{\lambda},\boldsymbol{\nu}\}$ as latent variables and $\{\boldsymbol{\gamma},\boldsymbol{\omega}\}$ as parameters to optimize. The conditional likelihood of the model in (7) involves a $(n+m-1)$-dimensional integral,

[TABLE]

which is nonanalytic. Alternatively, we take a Gaussian variational approximation (GVA) of the posteriors of random effect variables $\boldsymbol{U}:=\{a,\boldsymbol{b},\boldsymbol{c}\}$, which provides a lower bound $\underline{\ell}(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$ of ${\ell}(\boldsymbol{\gamma},\boldsymbol{\omega})$. Statistical properties of GVA for generalized linear mixed models are studied in Hall et al. (2011a, b) and Ormerod and Wand (2012), from a likelihood-based perspective. The resulting GVA estimator differs from the MLE but is asymptotically valid.

In our setting, the alternating steps are guaranteed to increase the following objective function (Neal and Hinton, 1998),

[TABLE]

where $q(\boldsymbol{\Psi},\boldsymbol{U})=q(\boldsymbol{\Psi})q_{\boldsymbol{\zeta},\boldsymbol{\kappa}}(\boldsymbol{U})$ naturally decouples into a factorized form, in which $q(\boldsymbol{\Psi})$ is left in free-form and $q_{\boldsymbol{\zeta},\boldsymbol{\kappa}}(\boldsymbol{U})$ is parameterized as Gaussian with diagonal covariance. With $t$ indexing the iterations, the overall algorithm contains the following alternating steps:

• •

E-step: Optimize $\mathcal{F}(q,\boldsymbol{\gamma},\boldsymbol{\omega})$ w.r.t. the distribution of latent variables $q(\boldsymbol{\Psi})$

[TABLE]

• •

Variational E-step: Update the Gaussian variational parameters $\{\boldsymbol{\zeta},\boldsymbol{\kappa}\}$ such that

[TABLE]

• •

M-step. Update the model parameters $\{\boldsymbol{\gamma},\boldsymbol{\omega}\}$ such that

[TABLE]

In this algorithm, the variational parameters $\{\boldsymbol{\zeta},\boldsymbol{\kappa}\}$ and the model parameters $\{\boldsymbol{\gamma},\boldsymbol{\omega}\}$ are updated with gradient updates instead of exact maximization (Lange, 1995a, b). Therefore, it is a generalized EM algorithm (Dempster et al., 1977; Neal and Hinton, 1998), as both the E-step and M-step are taken partially. Note that the latent variables $\boldsymbol{\Psi}$ only appear in the prior, and the random effect terms $\boldsymbol{U}$ only appear in the likelihood,

[TABLE]

so we can discuss them separately.

4.2 Closed-Form Expectations in the Shrinkage Prior

We compute the expected value w.r.t $\boldsymbol{\Psi}$ in the complete log-posterior, given the current parameter estimates and the observed data. Note that the entropy term does not depend on $(\boldsymbol{\gamma},\boldsymbol{\omega})$, so the relevant term in the E-step is

[TABLE]

where

[TABLE]

According to the Gaussian scale mixture (GSM) representation of the GDP prior,

[TABLE]

and denoting the pairwise differences as $\delta_{r,2\ell-1,2\ell}:=\gamma_{r,2\ell-1}-\gamma_{r,2\ell}$, $r=1,\ldots,h$, we have

[TABLE]

Given the estimates from the previous iteration ${\boldsymbol{\gamma}}^{(t)}$, the conditional posterior of latent variables $(\boldsymbol{\tau},\boldsymbol{\lambda})$ factorizes as

[TABLE]

where $({\tau}_{l}|{\lambda}_{l},-)\sim\mathrm{GIG}(0.5,\lambda_{l}^{2},\gamma_{l}^{2})$. Integrating out $\tau_{l}$, we have $(\gamma_{l}|\lambda_{l})\sim\mathrm{DE}(\gamma_{l};0,1/\lambda_{l})$ and $(\lambda_{l}|\alpha_{1},\eta_{1})\sim\mathrm{Ga}(\lambda_{l};\alpha_{1},\eta_{1})$, so $(\lambda_{l}|\gamma_{l},\alpha_{1},\eta_{1})\sim\mathrm{Ga}(\alpha_{1}+1,{|\gamma_{l}|}+\eta_{1})$, where $\mathrm{DE}(x;\mu=0,b)=\exp{(-|x|/b)}/{2b}$ refers to the Laplace distribution with scale parameter $b=1/\lambda$ and $\textrm{GIG}(x;a,b,p)$ denotes the Generalized Inverse Gaussian (GIG) distribution, $\textrm{GIG}(x;a,b,p)={0.5(a/b)^{{p}/{2}}}x^{p-1}\exp\left(-(ax+{b}/{x})/2\right)/K_{p}(\sqrt{ab})$, $(x>0)$, and $K_{p}(\theta)$ is the modified Bessel function of the second kind.

Similarly, given the estimates of the ${\boldsymbol{\delta}}^{(t)}$, the conditional posterior of latent variables $(\boldsymbol{\phi},\boldsymbol{\nu})$ also factorizes. Thus for every $(u,w)\in\mathcal{C}_{s^{\prime},\ell}$, $(s^{\prime},\ell)\in\mathcal{I}_{h-1}$, we have $({\phi}_{u,w}|{\nu}_{u,w},-)\sim\mathrm{GIG}(0.5,\nu_{u,w}^{2},\delta_{u,w}^{2})$, and integrating out $\phi_{u,w}$, we have $(\delta_{u,w}|\nu_{u,w})\sim\mathrm{DE}(\delta_{u,w};0,1/\nu_{u,w})$ and $(\nu_{u,w}|\alpha_{2},\eta_{2})\sim\mathrm{Ga}(\nu_{u,w};\alpha_{2},\eta_{2})$, so $(\nu_{u,w}|\delta_{u,w},\alpha_{2},\eta_{2})\sim\mathrm{Ga}(\alpha_{2}+1,{|\delta_{u,w}|}+\eta_{2})$. Therefore,

[TABLE]

We only need to find $\langle{\tau_{l}}^{-1}\rangle:=\langle{\rho_{l}}\rangle$. According to the change of variable formula, $f(\rho_{l})=\mathrm{GIG}(\tau_{l}^{-1};p,a,b)\rho_{l}^{-1}=\mathrm{GIG}(\rho_{l};-0.5,b,a)=\mathrm{InvGau}(\rho_{l};\sqrt{{\lambda_{l}^{2}}/{{\gamma_{l}}^{2(t)}}},\lambda_{l}^{2})$, we have

[TABLE]

and similarly denoting $\langle{\phi_{u,w}}^{-1}\rangle:=\langle{\upsilon_{u,w}}\rangle$, we obtain

[TABLE]

The GSM representation of the GDP priors determines a reweighting rule, as shown in (10) and (11), in which the weights depend only on the current estimate of the parameter $\boldsymbol{\gamma}^{(t)}$ (or its differences), and the hyperparameters $(\alpha_{1},\eta_{1},\alpha_{2},\eta_{2})$.

The penalty terms in (8) can be organized in a quadratic form, where the block-diagonal matrix $\widetilde{\boldsymbol{\Lambda}}$ can again be found by square completing,

[TABLE]

where

[TABLE]

Therefore, $\langle\ell_{1}(\boldsymbol{\Psi};\boldsymbol{\gamma})\rangle_{q(\boldsymbol{\Psi})}=\langle\ell_{1}(\boldsymbol{\Psi};\boldsymbol{\gamma})\rangle_{p(\boldsymbol{\Psi}|\boldsymbol{\gamma}^{(t)},-)}=-\boldsymbol{\gamma}^{T}\widetilde{\boldsymbol{\Lambda}}\boldsymbol{\gamma}/2$, which is the structured penalty that favors models with simpler structures conforming to $\mathcal{T}_{h}$. The quadratic form is concave and differentiable, which makes gradient-based optimization methods suitable.

4.3 Gaussian Variational Approximation of the Likelihood

In our case, the log-likelihood for the Poisson mixed model in (7) is

[TABLE]

The log-priors for the random effect terms are

[TABLE]

[TABLE]

[TABLE]

Let $q({a})={\mathcal{N}}(\zeta^{a},\kappa^{a})$, $q(b_{i})={\mathcal{N}}(\zeta^{b}_{i},\kappa^{b}_{i})$, $q(c_{j})={\mathcal{N}}(\zeta^{c}_{j},\kappa^{c}_{j})$111Note that for $b_{1}$ and $c_{1}$, we have fixed their value to [math] therefore for notational convenience, we assume $\zeta^{b}_{1}=\zeta^{c}_{1}=0$ and $\kappa^{b}_{1}=\kappa^{c}_{1}=0$ in the likelihood term. , where $\{\zeta^{a},\boldsymbol{\zeta}^{b},\boldsymbol{\zeta}^{c}\}$ are the mean parameters and $\{\kappa^{a},\boldsymbol{\kappa}^{b},\boldsymbol{\kappa}^{c}\}$ are all positive parameters for the variances. Assuming the variational proposals are independent, the lower bound of $\ln{p(\boldsymbol{X}|\boldsymbol{y},\boldsymbol{\gamma},\boldsymbol{\omega})}$ is

[TABLE]

Denoting $\underline{\ell}(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})=\langle\ell_{2}(\boldsymbol{U};\boldsymbol{\gamma},\boldsymbol{\omega})\rangle_{q(\boldsymbol{U})}$, we have

[TABLE]

So in the variational E-step and the M-step, we update the variational parameters $\{\boldsymbol{\zeta},\boldsymbol{\kappa}\}$ and the model parameters $\boldsymbol{\gamma}$ through a quasi-Newton method with objective function $Q(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa}):=\underline{\ell}(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})-\boldsymbol{\gamma}^{T}\widetilde{\boldsymbol{\Lambda}}\boldsymbol{\gamma}/2$, which only requires us to specify the first-order gradients (detailed in Appendix B).

In each M-step, we can also choose to optimize the prior parameter $\boldsymbol{\omega}$ via a fixed-point update. Setting the gradients $D_{\omega_{a}}{Q}=D_{\omega_{b}}{Q}=D_{\omega_{c}}{Q}=0$, we obtain,

[TABLE]

As suggested in Armagan et al. (2013), the hyper-parameters $\{\alpha_{1},\eta_{1},\alpha_{2},\eta_{2}\}$ can be either fixed or pre-learned from an initial Bayesian analysis based on griddy Gibbs sampling (Ritter and Tanner, 1992).

4.4 Computational Complexity

Our spinlets method is scalable to handle millions of POs. The complexity analysis and the running time of the recursive graph partitioning algorithm are detailed in Appendix A. In balancing the per-iteration cost with the convergence rate in the variational EM algorithm, we adopt a quasi-Newton method with the L-BFGS algorithm (Liu and Nocedal, 1989), which uses a predetermined $c_{0}=100$ number of previous steps to form a low-rank Hessian approximation with complexity $O(mc_{0})$. As will be illustrated in Figure 6, the variational EM algorithm converges very fast in practice. We measure the CPU time of these procedures on a standard laptop computer (Macbook Air, 1.6 GHz Intel Core i5, 8 GB 1600 MHz DDR3, Intel HD Graphics 6000 1536 MB). For the StatsBomb World Cup data with $h=9$ and $n=128$, each iteration of the variational EM algorithm takes about $14.91$ seconds to run.

5 Applications

In this section, we first compare the performance of various penalties within the GDP and fGDP families in a number of simulated examples (Section 5.1), and then use spinlets as an exploratory tool for visualizing spatial interaction networks from the FIFA World Cup 2018 dataset, while using the supervised information of interest (Section 5.2).

5.1 Simulation Study

We generate $n=\{25,50,100,200\}$ observations from the model $x_{i,j}\sim\textrm{Poisson}(\mu_{i,j})$, $\mu_{i,j}=t_{i}e^{\eta_{i,j}}$, $\eta_{i,j}=a+b_{i}+c_{j}+y_{i}\beta_{j}^{*}$, $t_{i}\sim\textrm{Ga}(2,1)$, $y_{i}\sim\textrm{Poisson}(0.5)$, $i=1,\ldots,n$, $j=1,\ldots,m$, and $a\sim{\mathcal{N}}(0,0.1)$, $b_{i^{\prime}}\sim{\mathcal{N}}(0,0.1)$, $c_{j^{\prime}}\sim{\mathcal{N}}(0,0.1)$, for $i^{\prime}=2,\ldots,n$, $j^{\prime}=2,\ldots,m$, $m=2^{h}$, $h=5$. We assume $\boldsymbol{\beta}^{*}$ to be signals with multiscale structures in the following configurations:

1. (a)

$\boldsymbol{\beta}^{*}=(1,1,0,0,1,1,0,0,1,1,-1,-1,0,0,-1,-1,1,1,0,0,1,1,0,0,-1,-1,1,1,0,0,1,1)$, 2. (b)

$\boldsymbol{\beta}^{*}=(\underbrace{1,\ldots,1}_{4},\underbrace{0,\ldots,0}_{4},\underbrace{-1,\ldots,-1}_{4},\underbrace{0,\ldots,0}_{4},\underbrace{1,\ldots,1}_{4},\underbrace{-1,\ldots,-1}_{4},\underbrace{0,\ldots,0}_{4},\underbrace{1,\ldots,1}_{4})$, 3. (c)

$\boldsymbol{\beta}^{*}=(\underbrace{1,\ldots,1}_{8},\underbrace{0,\ldots,0}_{8},\underbrace{-1,\ldots,-1}_{8},\underbrace{0,\ldots,0}_{8})$, 4. (d)

$\boldsymbol{\beta}^{*}=(1,1,0,0,\underbrace{-1,\ldots,-1}_{4},,\underbrace{0,\ldots,0}_{8},\underbrace{1,\ldots,1}_{16})$.

5.1.1 Baseline Methods

Denoting the GDP prior as $\textrm{GDP}(\alpha,\eta)$, we consider the following baseline methods (detailed derivations provided in Appendix C) operating on the original parameter $\boldsymbol{\beta}$:

1. (1)

GDP-0: no prior, $\beta_{j}\sim\textrm{GDP}(-1,1)$, and $j=1,\ldots,m$, 2. (2)

GDP prior with default parameters, $\beta_{j}\sim\textrm{GDP}(1,1)$, and $j=1,\ldots,m$, 3. (3)

Fused lasso signal approximation (FLSA) (Friedman et al., 2007) implemented with the GDP prior: $\beta_{j}\sim\textrm{GDP}(1,1)$, $\beta_{j^{\prime}}-\beta_{j^{\prime}+1}\sim\textrm{GDP}(1,1)$, and $j^{\prime}=1,\ldots,m-1$, 4. (4)

Pairwise fused lasso (PFL-S) (Petry et al., 2011) with the deletion and fusion penalties weighted by $0.8$ and $0.2$ respectively, implemented with the GDP prior: $\beta_{j}\sim\textrm{GDP}(1,1)$, and $\beta_{j}-\beta_{k}\sim\textrm{GDP}(1,1)$, $j,k=1,\ldots,m$, and $j\neq k$. 5. (5)

Pairwise fused lasso (PFL-F) (Petry et al., 2011) with the deletion and fusion penalties weighted by $0.2$ and $0.8$ respectively (implementation is similar as above).

We also consider the following structural regularizers on the expanded parameter $\boldsymbol{\gamma}$. All of them are derived from the $\textrm{fGDP}(\alpha_{1},\eta_{1},\alpha_{2},\eta_{2})$ family:

1. (6)

fGDP-S: deletion only, $\boldsymbol{\gamma}\sim\textrm{fGDP}(1,1,-1,1)$, 2. (7)

fGDP-F: fusion only, $\boldsymbol{\gamma}\sim\textrm{fGDP}(-1,1,1,1)$, 3. (8)

fGDP: with default parameters, $\boldsymbol{\gamma}\sim\textrm{fGDP}(1,1,1,1)$, 4. (9)

fGDP-NJ: with Normal-Jeffrey’s parameters, $\boldsymbol{\gamma}\sim\textrm{fGDP}(0,0,0,0)$.

To investigate the sensitivity to the hyperparameters $\{\alpha_{1},\eta_{1},\alpha_{2},\eta_{2}\}$, we additionally include: (10) fGDP1: $\boldsymbol{\gamma}\sim\textrm{fGDP}(1,0.1,1,0.1)$, (11) fGDP2: $\boldsymbol{\gamma}\sim\textrm{fGDP}(1,0.01,1,0.01)$, (12) fGDP3: $\boldsymbol{\gamma}\sim\textrm{fGDP}(1,0.001,1,0.001)$, (13) fGDP4: $\boldsymbol{\gamma}\sim\textrm{fGDP}(0.5,0.01,0.5,0.01)$, (14) fGDP5: $\boldsymbol{\gamma}\sim\textrm{fGDP}(2,0.01,2,0.01)$, and (15) fGDP6: $\boldsymbol{\gamma}\sim\textrm{fGDP}(5,0.01,5,0.01)$. All these $15$ methods share the same random initialization of parameters. The TFL penalties (Wang and Zhao, 2017) are not considered, as they only deal with tree-guided variable fusion and are limited to linear models in a completely different forward regression setting.

5.1.2 Performance Evaluation

We perform $50$ replications for every simulation scenario. To evaluate the performance, we use the F1 score = 2$\times$(precision$\times$ recall)/(precision+recall) as a metric. In particular, for selection, we examine whether the non-zero regression coefficients in the $m$ leaf groups are detected. For fusion, we make $2^{h}-1$ binary decisions on whether the regression coefficients within the descendant set $\mathcal{D}_{s^{\prime},\ell}$ are all equal, for every internal node $(s^{\prime},\ell)\in\mathcal{I}_{h-1}$. We also provide the relative mean square error (RMSE) as a reference metric for signal recovery $\textrm{RMSE}=||\widehat{\boldsymbol{\beta}}-\boldsymbol{\beta}^{*}||_{\textrm{F}}/||\boldsymbol{\beta}^{*}||_{\textrm{F}}$.

The results are illustrated in Figure 5. The results across rows demonstrate the adaptation ability of our PX-based fGDP approaches in recovering signals with multi-level and multiplex smoothness. In general, the performance improves with the replicate size. When replicate size is relatively small, regularization helps boost the performance. Methods encouraging fusion produce better results in terms of F1 score (fusion) than those only encouraging selection. Note that by encouraging the expanded parameters on the root-to-leaf path to be sparse, the fGDP-S method is also able to encourage fusion in an indirect way.

5.1.3 Hyperparameter Sensitivity

We observe that the default choice $\alpha=\eta=1$ and the Normal-Jeffrey’s choice $\alpha=\eta=0$ yield sub-optimal results. As the parameters in the PX space are redundant, a high level of shrinkage is required, demanding a small $\eta$ parameter. The PX-based fGDP approaches (fGDP1 $\sim$ fGDP6) with $\eta\in(10^{-3},10^{-1})$ compare favorably to the other methods in terms of both F1 scores (fusion and selection) and RMSE. Due to the adaptive shrinkage mechanism, the $\alpha$ parameters are less sensitive, we observed that $\alpha\in(0.5,5)$ works reasonably well. In the original parameter space, the only method standing out in terms of the F1 scores is the fusion-dominated pairwise fused lasso (PFL-F) approach; however, this approach suffers from the inaccurate prior knowledge discussed in Section 2.2. As a result, it yields the worst RMSE performance among all the baselines considered.

5.1.4 Convergence

We choose the convergence criteria to be $||\boldsymbol{\beta}^{(t)}-\boldsymbol{\beta}^{(t-1)}||_{\textrm{F}}<10^{-6}$, and for the PX-based method, this is $||\boldsymbol{D}\boldsymbol{\gamma}^{(t)}-\boldsymbol{D}\boldsymbol{\gamma}^{(t-1)}||_{\textrm{F}}<10^{-6}$. The quasi-Newton updates are performed with L-BFGS, with $1,000$ the maximum number of iterations allowed. We use the MATLAB solver minFunc (Schmidt, 2005). Setting the maximum iterations of variational EM to be $50$, the number of iterations until convergence is reported in Figure 6. As expected, the PX-based approaches converge much faster than the non-PX based approaches.

The PX-based approaches only converge to one of the many local optima. We repeated the fGDP2 approach with hyper-parameters $\alpha_{1}=\alpha_{2}=1$, $\eta_{1}=\eta_{2}=0.01$ with $500$ different random initializations. For every pair of different runs $a,b\in\{1,\ldots,500\}$ and $a\neq b$, we calculate the pairwise distances $||\widehat{\boldsymbol{\beta}}_{a}-\widehat{\boldsymbol{\beta}}_{b}||_{\textrm{F}}$ and $||\widehat{\boldsymbol{\gamma}}_{a}-\widehat{\boldsymbol{\gamma}}_{b}||_{\textrm{F}}$ as a metric for variability. The results under configuration (d) with $n=200$ are summarized in Figure 7, which show that the algorithm converges to a large number of different local optimal solutions in the auxiliary space; however, when reducing to the original space, they all map to solutions with comparable quality lying within a close neighborhood around the ground truth $\boldsymbol{\beta}^{*}$.

5.2 Supervised Dimension Reduction of Soccer Passing Networks

As a team sport, soccer is characterized by its free-flowing nature (Gudmundsson and Horton, 2017). The spatial interaction networks provide a concise abstraction of the team play on the pitch, which capture the essence of soccer as an invasion-territorial sport. Our spinlets approach provides a useful tool for further reducing complexity, taking advantage of both the passing network structure and the response relevance. On the FIFA World Cup 2018 dataset, we first consider two types of response in the two tasks below:

Task 1: Team performance measured by the goal difference (i.e., goal scored minus goal conceded) in the $128$ game plays (unique team-game pairs). 2. 2.

Task 2: The urgency and tiredness status of the game, indicated by whether the data are collected after $70$ minutes. This division of game play into two game phases results in $128\times 2=256$ replicates with binary responses.

We compare three algorithms: (1) rMETIS with $h=9$, (2) a single-scale selection-only (S3O) approach with the default GDP prior on $\boldsymbol{\beta}$, and (3) Spinlets with the fGDP2 prior on $\boldsymbol{\gamma}$. The threshold for the regression coefficients to be considered as (approximately) equal is set to be $0.005$. The rMETIS algorithm provides a common starting point of CG representation for both S3O and Spinlets. Interestingly, the same game can be seen with different eyes, when the practitioner picks a different response variable. In task 1 we consider team performance measured by the goal differences, while in task 2 we consider whether the passing network is collected in a relatively late phase of the game, when the time is running out and the players are getting physically and mentally tired. Figure 8 and 9 illustrate the reductive representations of exactly the same games with two different types of response. Spinlets merges functionally similar passes together in the same group, removes the ones irrelevant to the response, and produces visual outputs that are readily interpretable to the human eye.

The magnitude of the estimated regression coefficients $\widehat{\boldsymbol{\beta}}$ indicates the strength of association with the response. For the sake of clarity, we display them in separate plots in Figure 10. With both selection and fusion, Spinlets produces more concise representations than S3O and allows for non-uniform resolutions. Note that the number of passes in the different partition sets defines a passing strategy — this strategy can be related to the competition outcome or a situational game factor. These plots find the common pattern from all the games. For example, Figure 10 (b) suggests that the cross-passes from the two wings seem relatively inefficient in winning the game, or adopted more by the losing team. In contrast, Figure 10 (d) indicates that teams tend to control the midfield more in the first $70$ minutes, while pushing for the goal by passing more in the opposing half after $70$ minutes.

To illustrate the information preserved in the reductive representation, we plot the SDR score $R_{\widehat{\boldsymbol{\beta}}}(\boldsymbol{x}_{i})$ versus $y_{i}$ in Figure 11. For the PIR mixed model, the sufficiency of the SDR score holds conditionally on the replicate-specific random effects $[a,a+\boldsymbol{b}]$, so we cut the estimated mean of random effects $[{\zeta}_{a},{\zeta}_{a}+\boldsymbol{\zeta}_{b}]$ into three intervals, and present the results in three separate panels. Both S3O and Spinlets are regularized approaches, in which fitting the data is not the only goal. Comparing to S3O, Spinlets produces more parsimonious results with little or no sacrifice in discriminative power, as shown in Figure 11.

5.3 Predictive Powers of the Vectorial Representation

Different vectorial representation of the SPIN data can be obtained via different operations on the partitioning tree. For example, the rMETIS method cuts the tree at a single height, the S3O method prunes irrelevant leaf groups, and our Spinlets method enables both deletion and fusion across multiple scales. Although not the main focus of this article, the predictive powers of the vectorial representations induced by these methods are also evaluated. In the dimension reduction step, we first perform rMETIS with $h=9$ for all the data, and then split them into a training set and a testing set in a ratio of $90:10$. To refine the initial vectorial representation, the tree is cut at $s=7$ or $s=8$ unsupervisedly, or trimmed by S3O or Spinlets with the number of goals scored as supervision. The vectorial representations on the testing set can be determined accordingly.

Next, these vectorial representations are plugged in a forward Poisson regression model as predictors with the number of goals scored as the response. Specifically, we use the glmnet implementation (Friedman et al., 2010) of the Poisson log-linear regression model penalized by lasso ($\alpha=1$), elastic-net ($\alpha=0.5$), and ridge ($\alpha=0$). For each case, we perform $10$-fold cross-validation on the training data, and evaluate the performance on the held-out test data, both using the mean absolute error (MAE) metric. We conduct $100$ replicated simulations and summarize the mean and standard deviation of MAE in Table 1. Predicting the number of goals scored from passing strategies is, admittedly, a challenging task, and our spinlets method gives slightly improved generalization performance.

6 Discussion

In this article, we have introduced spinlets—a supervised dimension reduction method for spatial interaction networks using information on the spatial locations of each pass and also a response variable in constructing a multiscale representation. Particularly, the dimension reduction is conducted via a top-down partitioning of the similarity graph and a bottom-up pruning of the partition tree. Instead of cutting at a single height of a given tree, we select multiple tree heights for different branches of the tree adaptively, which yields representations with mixed granularities. The regularization prevents the information preserved in the lower-dimensional representation from being dominated by the supervisory signal without enough conformity to the spatial network data structures. In addition, our approach can be interpreted as an empirical Bayes approach, which estimates a hierarchical tree organization of the data in a first stage.

Besides the sports application studied in this article, our spinlets approach can accommodate massive-scale network predictors or temporally-indexed predictors with high sampling rate; both are pressing needs in neuroscience. Potentially further improvements on the flexibility and utility of spinlets are possible, such as zero-inflated variants, SDR with multiple responses and covariate adjustment.

7 Acknowledgment

This research was supported by grant W911NF-16-1-0544 from the U.S. Army Research Institute for the Behavioral and Social Sciences (ARI). The authors thank StatsBomb for making the data freely available in public for academic research. The authors would also like to thank Li Ma and Andrés Felipe Barrientos for inspirational discussions.

A Preprocessing—Recursive METIS (rMETIS)

To concisely represent the proximity information of POs (e.g., soccer passes), we define a tree structure in three steps: $(i)$ choose the Euclidean distance metric between pairs of POs and compute a $Q\times Q$ distance matrix, $Q=\sum_{i}q_{i}$; $(ii)$ construct a sparse similarity graph of POs $\mathcal{G}$ by considering $K$ nearest neighbors; $(iii)$ build a partition tree $\mathcal{T}_{h}$ via recursively applying the METIS partitioning algorithm (Karypis and Kumar, 1998) on $\mathcal{G}$. We construct a full binary tree via recursive METIS (rMETIS). In each step, a set of POs is split into two disjoint subsets. Figure 12 illustrates a sub-branch of $\mathcal{T}_{10}$ with all assigned POs in the $16$ groups plotted.

The $Q\times Q$ distance matrix is constructed approximately by the k-d tree nearest neighbor search, which has $O(Q)$ worst case complexity. Second, the rMETIS algorithm assigns the POs into $m$ primary groups. After $h$ phases, the set of POs $\mathcal{P}$ is partitioned into $m=2^{h}$ leaf groups on $\mathcal{T}_{h}$. The complexity of the bisection algorithm is ${O}(|E|\log{m})$, where $|E|$ is the number of edges in the similarity graph. For $Q=49,988$, constructing the nearest neighbor set (with $K=1,500$ nearest neighbors) requires $15.97$ seconds, computing the similarity matrix costs $68.58$ seconds, and running the rMETIS algorithm (with tree depth $h=9$) takes $52.87$ seconds.

B Gradients in the M-Step

The variational parameters $\boldsymbol{\zeta}$ and $\boldsymbol{\kappa}$ should be chosen to make $\underline{\ell}(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$ as close as possible to $\ell(\boldsymbol{\gamma},\boldsymbol{\omega})$. Given $q(\boldsymbol{\Psi})$ and denoting $\mathcal{F}(q(\boldsymbol{\Psi}),\boldsymbol{\gamma},\boldsymbol{\omega})=Q(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$, we have

• •

Gradients of $\boldsymbol{\gamma}$: Denoting the gradients of $Q(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$ over $\gamma_{l}$ as $D_{\gamma_{l}}{Q}$, we have

[TABLE]

where $\tilde{x}_{i,j}=t_{i}\exp{[\zeta^{a}+\zeta^{b}_{i}+\zeta^{c}_{j}+\frac{1}{2}(\kappa^{a}+\kappa^{b}_{i}+\kappa^{c}_{j})+y_{i}\sum_{l=1}^{L}{d}_{j,l}{\gamma}_{l}]}$.

• •

Gradients of $\boldsymbol{\zeta}$: Denoting the gradients of $Q(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$ over $\zeta^{a}$, $\zeta_{i}^{b}$, and $\zeta_{j}^{c}$ as $D_{\zeta^{a}}{Q}$, $D_{\zeta_{i}^{b}}{Q}$, and $D_{\zeta_{j}^{c}}{Q}$ respectively, we have

[TABLE]

• •

Gradients of $\boldsymbol{k}$: To ensure positivity of the variance parameters $\boldsymbol{\kappa}$, we apply the $\log$ reparameterization on $\boldsymbol{\kappa}$. Denoting the gradients of $Q(\boldsymbol{\gamma},\boldsymbol{\omega},\boldsymbol{\zeta},\boldsymbol{\kappa})$ over $k^{a}=\log{(\kappa^{a})}$, $k_{i}^{b}=\log{(\kappa_{i}^{b})}$, and $k_{j}^{c}=\log{(\kappa_{j}^{c})}$ as $D_{k^{a}}{Q}$, $D_{k_{i}^{b}}{Q}$, and $D_{k_{j}^{c}}{Q}$ respectively, we have

[TABLE]

C GDP Family Priors on $\boldsymbol{\beta}$

The GDP Prior

The GDP penalty term is $\log{p(\boldsymbol{\beta})}=\sum_{j=1}^{m}\log{p({\beta}_{j})}=\sum_{j=1}^{m}[-\log{(2\xi)}-(\alpha+1)\log{(1+\frac{|{\beta}_{j}|}{\alpha\xi})]}$. In the E-Step, we have $\langle{\rho_{j}}\rangle={(\alpha+1)}/{[|{\beta_{j}}^{(t)}|({|{\beta_{j}}^{(t)}|}+\eta)]}$. In the M-step, the gradient w.r.t $\beta_{j}$ can be written as,

[TABLE]

where $\tilde{x}_{i,j}:=t_{i}\exp{[\zeta^{a}+\zeta^{b}_{i}+\zeta^{c}_{j}+\frac{1}{2}(\kappa^{a}+\kappa^{b}_{i}+\kappa^{c}_{j})+y_{i}\beta_{j}]}$ and the gradients w.r.t other parameters do not change.

GDP-based FLSA

The GDP-based FLSA penalty term is

[TABLE]

In the E-Step, we have

[TABLE]

where $\delta_{j^{\prime}}^{(t)}={\beta_{j^{\prime}}}^{(t)}-{\beta_{j^{\prime}+1}}^{(t)}$. In the M-step, the gradient w.r.t $\beta_{j}$ can be written as $D_{\beta_{j}}{Q}=\sum_{i=1}^{n}\epsilon_{i,j}y_{i}-[\widetilde{\boldsymbol{\Lambda}}^{(t)}\boldsymbol{\beta}]_{j}$, where

[TABLE]

GDP-based PFL

The GDP-based PFL penalty term is

[TABLE]

where $\theta\in[0,1]$ is the weight parameter that balances selection with fusion. In the E-Step, we have

[TABLE]

where $\delta_{j,k}^{(t)}={\beta_{j}}^{(t)}-{\beta_{k}}^{(t)}$. Similarly, in the M-step the gradient w.r.t $\beta_{j}$ can be written as $D_{\beta_{j}}{Q}=\sum_{i=1}^{n}\epsilon_{i,j}y_{i}-[\widetilde{\boldsymbol{\Lambda}}^{(t)}\boldsymbol{\beta}]_{j}$, where

[TABLE]