Bi-clustering for time-varying relational count data analysis

TL;DR

This paper introduces two novel bi-clustering models, dPIRM and a zero-inflated extension, for analyzing time-varying relational count data, effectively capturing latent structures and temporal dynamics.

Contribution

The paper proposes two new models for bi-clustering relational count data over time, including a zero-inflated extension for sparse data, with detailed derivations and analysis.

Findings

Both models effectively capture temporal dynamics in relational data.

Simulation and real data demonstrate the models' effectiveness.

Zero-inflated model handles sparse, short-interval data better.

Abstract

Relational count data are often obtained from sources such as simultaneous purchase in online shops and social networking service information. Bi-clustering such relational count data reveals the latent structure of the relationship between objects such as household items or people. When relational count data observed at multiple time points are available, it is worthwhile incorporating the time structure into the bi-clustering result to understand how objects move between the cluster over time. In this paper, we propose two bi-clustering methods for analyzing time-varying relational count data. The first model, the dynamic Poisson infinite relational model (dPIRM), handles time-varying relational count data. In the second model, which we call the dynamic zero-inflated Poisson infinite relational model, we further extend the dPIRM so that it can handle zero-inflated data. Proposing both two models is important as zero-inflated data are often encountered, especially when the time intervals are short. In addition, by explicitly deriving the relevant full conditional distributions, we describe the features of the estimated parameters and, in turn, the relationship between the two models. We show the effectiveness of both models through a simulation study and a real data example.

Tables

Table 1. Table 1: RI results

$m$	$s$	PIRM	dPIRM	dZIPIRM
0.1	0.3	0.79(0.16)	0.89(0.14)	0.91(0.13)
	0.5	0.77(0.16)	0.94(0.10)	0.93(0.10)
	0.7	0.76(0.18)	0.90(0.14)	0.94(0.11)
0.2	0.3	0.77(0.16)	0.87(0.15)	0.96(0.07)
	0.5	0.76(0.16)	0.94(0.10)	0.93(0.11)
	0.7	0.76(0.16)	0.92(0.12)	0.96(0.08)
0.3	0.3	0.73(0.16)	0.89(0.14)	0.88(0.13)
	0.5	0.73(0.17)	0.89(0.14)	0.90(0.13)
	0.7	0.74(0.17)	0.87(0.15)	0.90(0.14)

Table 2. Table 2: Log-likelihood results

$m$	$s$	PIRM	dPIRM	dZIPIRM
0.1	0.3	-2.28(0.38)	-2.04(0.28)	-2.01(0.25)
	0.5	-2.23(0.40)	-2.02(0.27)	-2.02(0.27)
	0.7	-2.26(0.32)	-2.10(0.25)	-2.09(0.25)
0.2	0.3	-2.24(0.33)	-2.04(0.26)	-2.00(0.25)
	0.5	-2.20(0.33)	-2.02(0.23)	-2.04(0.25)
	0.7	-2.20(0.30)	-2.04(0.27)	-2.03(0.29)
0.3	0.3	-2.14(0.25)	-2.12(0.22)	-2.11(0.22)
	0.5	-2.22(0.25)	-2.09(0.27)	-2.08(0.28)
	0.7	-2.18(0.26)	-2.02(0.26)	-2.03(0.25)

Table 3. Table 3: Log-likelihood results for the real data

	PIRM	dPIRM	dZIPIRM
Log-likelihood	-6.88	-4.84	-3.84

Table 4. Table 4: Time-varying clustering assignments for selected objects in the dPIRM

	January	February	March	April	May	June
cluster 1	small item	small item	small item	small item	small item	small item
	tie	tie	tie	tie	tie	tie
	suit	suit	suit	suit	suit	suit
	casual shirt	casual shirt	casual shirt	casual shirt	casual shirt	casual shirt
	jeans		jeans	jeans	jeans	jeans
				slacks	slacks	slacks
				short sleeve T-shirt	short sleeve T-shirt	short sleeve T-shirt
						short pants
cluster 2			short pants	short pants
				long sleeve T-shirt
cluster 3	slacks	slacks
	short sleeve T-shirt	short sleeve T-shirt
cluster 4	short pants	short pants			long sleeve T-shirt	long sleeve T-shirt
cluster 5		jeans	slacks		short pants
			short sleeve T-shirt
cluster 6	long sleeve T-shirt	long sleeve T-shirt	long sleeve T-shirt

Table 5. Table 5: Time-varying clustering assignments for selected objects in the dZIPIRM

	January	February	March	April	May	June
cluster 1	small item	small item	small item	small item	small item	small item
	tie	tie	tie	tie	tie	tie
	suit	suit	suit	suit	suit	suit
			casual shirt	casual shirt	casual shirt	casual shirt
			jeans	jeans	jeans	jeans
				slacks	slacks	slacks
				short sleeve T-shirt	short sleeve T-shirt	short sleeve T-shirt
					short pants	short pants
cluster 2	casual shirt			short sleeve T-shirt	short sleeve T-shirt	short sleeve T-shirt
	jeans					short pants
cluster 4		casual shirt
		jeans
cluster 7			slacks
cluster 9		slacks
cluster 10	short sleeve T-shirt	short sleeve T-shirt
cluster 11			short sleeve T-shirt		short pants
cluster 12	slacks		long sleeve T-shirt	short pants
cluster 13	short pants	short pants	short pants	long sleeve T-shirt	long sleeve T-shirt	long sleeve T-shirt
	long sleeve T-shirt	long sleeve T-shirt

Equations

$$\mathbf{\beta }|\gamma$$

$${\mathbf{\pi }}_{tk}|{\alpha }_0,\kappa ,\mathbf{\beta }$$

$$z_{ti}|z_{(t-1)i},{\Pi }_t$$

$${\lambda }_{k\ell }|a,b$$

$$x_{tij}|Z_t,\{{\lambda }_{k\ell }\}$$

$$\mathbf{\beta }|\gamma$$

$${\mathbf{\pi }}_{tk}|{\alpha }_0,\kappa ,\mathbf{\beta }$$

$$z_{ti}|z_{(t-1)i},{\Pi }_t$$

$$w_{tij}|c,d$$

$$r_{tij}|w_{tij}$$

$${\lambda }_{k\ell }|a,b$$

$$x_{tij}|Z_t,\{{\lambda }_{k\ell }\},r_{tij}$$

$$P(x|\lambda ,w)=(1-w)I(x=0)+w\frac{{\lambda }^x{e}^{-\lambda }}{x!}.$$

$$r|w$$

$$x|\lambda ,r$$

$$P(x,r|\lambda ,w)$$

$$={\left(\frac{{\lambda }^x{e}^{-\lambda }}{x!}\right)}^{r}I(x=0{)}^{1-r}{w}^r(1-w{)}^{1-r}$$

$$\displaystyle=\biggl{\{}w\frac{\lambda^{x}e^{-\lambda}}{x!}\biggr{\}}^{r}\biggl{\{}(1-w)I(x=0)\biggr{\}}^{1-r}.$$

$$\sum _{r}P(x,r|\lambda ,w)=(1-w)I(x=0)+w\frac{{\lambda }^x{e}^{-\lambda }}{x!}.$$

$$h(\lambda |x,Z)$$

$$\propto {\lambda }_{k\ell }^{(a+{\sum }_{t,i,j\in C_{k\ell }}{x}_{tij})-1}{e}^{-[b+{\sum }_{t,i,j}I\{t,i,j\in C_{k\ell }\}]{\lambda }_{k\ell }}.$$

$${\lambda }_{k\ell }$$

$$h(\lambda |w,r,x,Z)$$

$$\propto {\lambda }_{k\ell }^{(a+{\sum }_{t,i,j\in C_{k\ell }}{r}_{tij}{x}_{tij})-1}{e}^{-[b+{\sum }_{t,i,j\in C_{k\ell }}{r}_{tij}]{\lambda }_{k\ell }}.$$

$${\lambda }_{k\ell }$$

$$\displaystyle\begin{split}x_{tij}=0&\Longrightarrow\left\{\begin{array}[]{ll}r_{tij}=1&(\mathrm{Poisson})\\ r_{tij}=0&(\mathrm{Degenerate})\end{array}\right.,\\ x_{tij}>0&\Longrightarrow\quad r_{tij}=1\,(\mathrm{Poisson}).\end{split}$$

$$\sum _{t,i,j\in C_{k\ell }}{r}_{tij}{x}_{tij}=\sum _{t,i,j\in C_{k\ell }}{x}_{tij},$$

$${\lambda }_{k\ell }$$

$$h(w|\lambda ,r,x)$$

$$\propto \prod _t\prod _i{w}_{tij}^{c+r_{tij}-1}(1-w_{tij}{)}^{d+1-r_{tij}-1}.$$

$$w_{tij}$$

$$h(r|x,\lambda ,w)$$

$$\times (I(x_{tij}=0){)}^{(1-r_{tij})}{w}_{tij}^{r_{tij}}(1-w_{tij}{)}^{1-r_{tij}}$$

$$\propto \prod _t\prod _i\prod _j\{w_{tij}{e}^{-{\lambda }_{z_{ti}{z}_{tj}}}({\lambda }_{z_{ti}{z}_{tj}}{)}^{x_{tij}}{\}}^{r_{tij}}$$

$$\times \{(1-w_{tij})I(x_{tij}=0){\}}^{(1-r_{tij})}.$$

$$r_{tij}$$

$$\begin{array}{rl}E\left[{\lambda }_{k\ell }^{(dP)}\right]& =\frac{a+{\sum }_{t,i,j\in C_{k\ell }}{x}_{tij}}{b+{\sum }_{t,i,j}I\{t,i,j\in C_{k\ell }\}},\\ V\left[{\lambda }_{k\ell }^{(dP)}\right]& =\frac{a+{\sum }_{t,i,j\in C_{k\ell }}{x}_{tij}}{(b+{\sum }_{t,i,j}I\{t,i,j\in C_{k\ell }\}{)}^2}.\end{array}$$

$$\begin{array}{rl}E\left[{\lambda }_{k\ell }^{(dZIP)}\right]& =\frac{a+{\sum }_{t,i,j\in C_{k\ell }}{x}_{tij}}{b+{\sum }_{t,i,j\in C_{k\ell }}{r}_{tij}},\\ V\left[{\lambda }_{k\ell }^{(dZIP)}\right]& =\frac{a+{\sum }_{t,i,j\in C_{k\ell }}{x}_{tij}}{(b+{\sum }_{t,i,j\in C_{k\ell }}{r}_{tij}{)}^2}.\end{array}$$

$$E\left[{\lambda }_{k\ell }^{(dP)}\right]$$

$$V\left[{\lambda }_{k\ell }^{(dP)}\right]$$

Full text

Keywords: bayesian model, bi-clustering, count data, time-varying relational data, zero-inflated poisson.

1 Introduction

Relational data contain information about the relationship between objects and may be obtained from, for example, point-of-sale (POS) information and the social networking service (SNS) community. These types of data are often represented as a matrix in which each element indicates the relationship between the row and column elements (e.g., the presence or absence of a relationship or its strength). Bi-clustering is one way to cluster elements along the row and column dimensions of a matrix simultaneously [14, 6]. By applying bi-clustering to relational data, we can ascertain the latent structure of the relationship. Indeed, many probabilistic models using bi-clustering of relational data have been proposed [9, 15, 8].

Another situation commonly found in practice is when relationships change over time. For example, in the case of simultaneous purchases of products, the relationships between them may change according to events such as the season. In the SNS community, on the contrary, human relationships change through events such as career moves. In this way, the relationship between objects often changes through an event. In this paper, we call relational data obtained at multiple time points “time-varying relational data,” as in [7].

A naive approach to bi-clustering time-varying relational data is to apply bi-clustering to the data at each time step separately. However, in practice, this leads to unstable clustering results [22], and it becomes difficult to interpret clusters across time points since cluster structures are inconsistent over time. On the contrary, clustering using a time-varying structure [20, 5] can help detect changes in relationships due to specific events.

Several bi-clustering methods exist for time-varying relational data [22, 20, 5]. First, the dynamic infinite relational model (dIRM) [7], which extends HMM-based method [18, 3] to relationship data, is especially useful since it can automatically determine the number of clusters from the data. The dIRM method can handle relational data containing either a zero or a one (essentially whether a relationship exists). However, in practice, we often obtain relational count data. For example, if each element of a data matrix represents how many times two products (corresponding to the row and column of the matrix) are simultaneously purchased, this data can be considered to be relational count data between the products. Since count data contain more information than binary data, it is preferable to use a model that can handle relational count data when we have count information instead of only zeros and ones.

Based on the foregoing, in this paper, we propose two new methods for the bi-clustering of count data under time changes. First, we extend the Poisson infinite relational model (PIRM) [21], a bi-clustering method for relational count data that does not consider bi-clustering relationships that change over time, to be able to handle time changes to a model. We call this first model the dynamic Poisson infinite relational model (dPIRM). Second, we propose a dPIRM that can handle data containing many zeros, which we call the dynamic zero-inflated Poisson infinite relational model (dZIPIRM). The latter extension is especially important for count data, as zero-inflated count data are often obtained in practice, especially when the interval between observation times is short [1, 11, 4].

The remainder of the paper is organized as follows. In Section 2, we begin by considering a bi-clustering model of time series count data without zero-inflated characteristics. Next, we propose a bi-clustering model for zero-inflated time-varying count data. In Section 3, we derive the posterior distributions of the relevant parameters and clarify their properties. In Section 4, we illustrate how our proposed models can handle numerically simulated data. In Section 5, we show how the two proposed bi-clustering methods may be used with real data. Finally, Section 6 concludes.

2 Bi-clustering for time-varying count data

In this section, we propose two bi-clustering methods. We first introduce one for time-varying count data and then the zero-inflated Poisson model (ZIP) [10] to extend the first method to be able to handle zero-inflated count data.

2.1 The dPIRM

We denote the data by X=\bigl{\{}x_{tij}\in\{0,1,2,\cdots\};1\leq i\leq N_{1},1\leq j\leq N_{2},1\leq t\leq T\bigr{\}}, where $N_{1}$ represents the number of row objects, $N_{2}$ represents the number of column objects, and $T$ represents the time horizon. We can assume that $N_{1}=N_{2}$ without loss of generality. In addition, since relational count data are not necessarily limited to an undirected graph, we allow asymmetric relations (i.e., $x_{tij}\neq x_{tji}$ for all $i,j\,(i\neq j)$).

We propose a bi-clustering model of time-varying count data under the following three assumptions:

(A1)

Cluster structures remain consistent over time.

(A2)

Objects can move between clusters at each time step.

(A3)

The strength of the relationship between clusters can be expressed from zero to infinity.

We describe this dPIRM as follows:

[TABLE]

Here, $\Pi_{t}=\{\bm{\pi}_{tk}:k=1,2,\cdots\}$, $z_{ti}=k,k\in\{1,2,\cdots\}$.

In Equation (1), the term “Stick” denotes a division formula known as the stick-breaking process [17]. This process can be described as $\beta_{k}=v_{k}\prod_{l=1}^{k-1}(1-v_{l})$, where $v_{k}\sim\mathrm{Beta}(1,\gamma)$. $\gamma>0$ is a shape parameter. $\bm{\beta}=(\beta_{1},\beta_{2},\cdots)$ represents time-averaged membership of the clusters (A1), and the $\bm{\beta}$ are guaranteed to satisfy $\sum_{k}\beta_{k}=1$.

DP in Equation (2) represents the Dirichlet process. $\alpha_{0}$ plays the role of a concentration parameter in the Dirichlet distribution. Thanks to the DP, the number of clusters is automatically determined. $\pi_{tk\ell}$ in $\bm{\pi}_{tk}=(\pi_{tk1}$, $\pi_{tk2}$, $\cdots$, $\pi_{tk\ell},\cdots)$ represents the transition probability that an object belonging to cluster $k$ at time $t-1$ moves to cluster $\ell$ at time $t$. $\pi_{tk\ell}$ is guaranteed to satisfy $\pi_{tk\ell}>0$ and $\sum_{\ell}\pi_{tk\ell}=1$. $\delta_{k}$ is an indicator function taking the value of one when it belongs to cluster $k$ and zero otherwise. $\kappa>0$ adjusts how likely it is for objects to stay in the same cluster at each time step. In Equation (3), cluster assignment $Z_{t}=\{z_{ti}\}_{i=1}^{N}$ is generated with a multinomial distribution using $\bm{\pi}_{tk}$ to incorporate the time structure. This enables objects to move between clusters over time, following (A2).

In Equation (4), $\lambda_{k\ell}$ represents the strength of the relationship between the cluster indicated by the $k$th row and $\ell$th column. $a$ and $b$ represent the shape and scale parameters of the Gamma distribution, respectively. By generating $\lambda_{k\ell}$ from a Gamma distribution, it is possible to represent relationships from zero to infinity, following (A3). In Equation (5), count data $x_{tij}$ are generated with the parameter $\lambda_{z_{ti}z_{tj}}$.

Next, we consider the relationship between existing models (i.e., the dIRM and PIRM) and the dPIRM. For the dIRM, the Gamma distribution for the prior of $\lambda_{k\ell}$ in Equation (4) of the dPIRM is replaced with a Beta distribution, while the Poisson distribution in Equation (5) is replaced with the Bernoulli distribution. On the contrary, for the PIRM, the parameter $\bm{\beta}$ in Equation (1) in the dPIRM is only used as a parameter for a multinomial distribution instead of $\bm{\pi}_{tz_{(t-1)i}}$ in Equation (3). Moreover, the PIRM framework does not consider the time structure as in Equation (2).

Using the dPIRM, we can handle time-varying relational count data. However, in empirical applications, zero-inflated count data are often encountered, especially when the interval between time steps is short. In such cases, the approximation accuracy of the model is significantly low, as the dPIRM cannot handle zero-inflated count data. For this reason, we propose a model to address this issue in the next subsection.

2.2 The dZIPIRM

In this section, we extend the dPIRM to handle zero-inflated count data by proposing the dZIPIRM. In addition to Assumptions (A1)–(A3), this extension further assumes the following:

(A4)

Zero-inflated data are expressed by assuming mixture distributions of count data and zero data.

Given time-varying count data $x_{tij}\,(1\leq i,j\leq N,1\leq t\leq T)$, the dZIPIRM model is defined as follows:

[TABLE]

We assume that Equations (6), (7), (8), and (11) follow the same generator processes as the dPIRM. Equations (9), (10), and (14) represent the extension of the dPIRM to zero-inflated count data. Before explaining the model in detail, we first describe how we handle zero-inflated data with this model.

First, we introduce the ZIP, which is often used to model zero-inflated data. In the Poisson distribution, we assume that zeros and nonzeros are obtained from the same generator process. On the contrary, in the ZIP, count data are generated from, mixture of a Poisson distribution with weight $w$ and a Degenerate distribution for $x=0$ with weight $1-w$. The Degenerate distribution for $x=0$ denotes a distribution of $P(x=0)=1$. Letting $x=x_{tij}$, $w=w_{tij}$, and $\lambda=\lambda_{z_{ti}z_{tj}}$, the probability function can be written as follows:

[TABLE]

$I(\cdot)$ is an indicator function.

A naive approach to extending the dPIRM for zero-inflated data is to simply assume that data $x$ follows the ZIP. However, this assumption causes a problem when estimating the parameters. Specifically, to estimate the parameters using Gibbs sampling, a technique often used in Bayesian models, we need to derive the full conditional distribution. In general, however, it tends to be difficult to derive full conditional distributions explicitly when the probability is represented as a sum, as in Equation (15) [2]. However, if we can derive the full conditional distribution, we can easily interpret the parameter properties. Therefore, to ensure tractability and interpretability, we introduce the latent variable $r$ [16, 12]. Letting $r=r_{tij}$, we consider a model by setting the distribution as follows:

[TABLE]

In Equation (19), the value of $r$ generated by Equation (16) determines whether the data follow the Poisson distribution or the Degenerate distribution for $x=0$. Thus, given $\lambda$ and $w$, the joint probability of $x$ and $r$ is

[TABLE]

In this way, the probability can be represented in the form of a product, and the joint probability of the parameters and data (i.e., $P(x,r,w,\lambda)$) is also expressed as a product. This solves the problem of the joint distribution in Equation (15) being expressed as a sum.

Furthermore, when we sum out all the possible realized values of $r$ in Equation (20), we obtain the formula on the left-hand side of the ZIP as follows:

[TABLE]

The resulting formula is equivalent to the ZIP in Equation (15). In other words, by assuming Equations (16) and (19), we obtain a model with the property of the ZIP while deriving an explicit posterior distribution. For this reason, the distributions in Equations (16) and (19) are assumed to hold in Equations (10) and (14).

Since Equations (6), (7), (8), and (11) are equivalent to the dPIRM, we explain Equations (9), (10), and (14) in detail. In Equation (9), we use the Beta distribution for the prior distribution of $w_{tij}$. Since the Beta distribution is the conjugate prior distribution of the Bernoulli distribution, we can analytically derive the posterior distribution, with $c>0$ and $d>0$ as the shape parameters. In Equations (10) and (14), $r_{tij}$ determines whether $x_{tij}$ follows the Degenerate distribution for $x_{tij}=0$ or the Poisson distribution, following (A4). This indicates that in this model, zero data, which is generated from the Degenerate distribution ($r_{tij}=0$), is not taken into account for estimation. In other words, when estimating the parameters in the dZIPIRM, we can remove the effect of zeros that follow the Degenerate distribution. This intuitive interpretation will be justified by Theorem 1 in Section 3.2.

Figure 1 illustrates the parameter relations in the PIRM, dPIRM, and dZIPIRM. The circle nodes denote the variables, square nodes are the hyperparameters, and shaded nodes indicate the observations. The figure illustrates the differences between the data generation processes considered in each model, showing that the dZIPIRM has a structure in which $x_{tij}$ depends on $r_{tij}$.

3 Sampling parameters

In this section, we derive the full conditional distribution of the parameter $\lambda_{k\ell}$ of the dPIRM as well as the parameters $\lambda_{k\ell}$, $w_{tij}$, and $r_{tij}$ of the dZIPIRM. We then explain their properties. Deriving the full conditional distributions allows us to understand the properties of the parameters when the other parameters and data are fixed. Moreover, by comparing $\lambda_{k\ell}$ for the dZIPIRM and dPIRM, we can appreciate how the models differ.

We use beam sampling [19, 7] to sample the parameters, as it is faster and more efficient than Gibbs sampling when there is strong series correlation between $z_{ti}$ and $z_{(t-1)i}$. Indeed, convergence is significantly slowed in this case for Gibbs sampling.

3.1 Derivation of the full conditional distribution

In this section, we derive the full conditional distribution of the parameters $\lambda_{k\ell}$ in the dPIRM and $\lambda_{k\ell}$, $w_{tij}$, and $r_{tij}$ in the dZIPIRM. Since the full conditional distribution of the other parameters in the dPIRM and dZIPIRM are the same as in the dIRM, we refer the reader to [7] for the derivation of those distributions. In the derivation of the full conditional distribution, it is sufficient to consider the proportionality relation between the parameters for the calculation. To simplify the notation, we define the following set, $C_{k\ell}=\{t,i,j;z_{ti}=k,z_{tj}=\ell\}$.

First, we explain the parameters of the dPIRM. $\lambda_{k\ell}$ in the dPIRM is represented by the following proportionality relation:

[TABLE]

Therefore, the full conditional distribution of $\lambda_{k\ell}$ in the dPIRM becomes

[TABLE]

Next, we derive the full conditional distributions of $\lambda_{k\ell}$, $w_{tij}$, and $r_{tij}$ in the dZIPIRM. Unlike the dPIRM, the derivation of these full conditional distributions is not straightforward. However, thanks to the formulation of the dZIPIRM in Section 2, setting conjugate prior distributions for the prior distribution of the parameters allows us to analytically derive the full conditional distribution of the parameters in the dZIPIRM. First, $\lambda_{k\ell}$ in the dZIPIRM is represented by the following proportionality relation:

[TABLE]

Following a similar reasoning as in Equation (21), the full conditional distribution of $\lambda_{k\ell}$ in the dZIPIRM may be expressed as follows:

[TABLE]

This expression may be further simplified as follows. From Equation (14):

[TABLE]

Using this, we obtain

[TABLE]

which we prove in Appendix. Using Equation (23), Equation (22) may be rewritten as follows:

[TABLE]

This expression is useful to compare the properties of the dPIRM and dZIPIRM, as explained in Section 3.2.

Next, for the parameter $w_{tij}$ of the dZIPIRM, we use the proportionality relation:

[TABLE]

Therefore, the full conditional distribution of $w_{tij}$ becomes the Beta distribution as follows:

[TABLE]

Finally, we consider the full conditional distribution of the latent variable $r_{tij}$, for which we consider the proportionality relation:

[TABLE]

The full conditional distribution of $r_{tij}$ becomes

[TABLE]

3.2 Interpreting the properties of the model using the full conditional distribution

In this section, we clarify and compare the properties of the two models by studying the full conditional distributions of the parameters derived in Section 3.1. First, we compare $\lambda_{k\ell}$ in the dPIRM and dZIPIRM to explain the properties of each model. The expected value and variance of $\lambda_{k\ell}$ in the dPIRM from Equation (21), denoted by $\lambda_{k\ell}^{(dP)}$, is expressed as follows:

[TABLE]

$\sum_{t,i,j\in C_{k\ell}}x_{tij}$ and $\sum_{t,i,j}I\{t,i,j\in C_{k\ell}\}$ represent the total value and number of $x_{tij}$ belonging to the $k$th row and $\ell$th column cluster, respectively. Ignoring $a$ and $b$, we can interpret Equation (26) as indicating that the expected value of $\lambda_{k\ell}$ is the average value of $x_{tij}$ belonging to the $k$th row and $\ell$th column cluster through time. Note that $a$ and $b$ are hyperparameters that adjust the average value and must be determined in advance. Further, from Equation (26), we can see how $a$ and $b$ affect the estimated $\lambda_{k\ell}$. Specifically, $\lambda_{k\ell}$ becomes larger when $a>b$, while $\lambda_{k\ell}$ becomes smaller when $a<b$. In addition, the larger the values of $a$ and $b$, the larger the amount of prior information. Thus, by deriving explicit relationship between parameters and hyperparameters, users can have intuition for how hyperparameters affect the estimation, and then it helps them select hyperparameters.

Next, we consider $\lambda_{k\ell}$ in Equation (24) of the dZIPIRM, denoted by $\lambda_{k\ell}^{(dZIP)}$. As before, the expected value and variance of $\lambda_{k\ell}^{(dZIP)}$ are derived to be

[TABLE]

The numerator in Equation (27) is the same as in Equation (26). $\sum_{t,i,j\in C_{k\ell}}r_{tij}$ represents the number of $x_{tij}$ generated from the Poisson distribution ($r_{tij}=1$), which belongs to the $k$th row and $\ell$th column cluster. Ignoring $a$ and $b$, Equation (27) indicates the average value of $x_{tij}$ considered to be generated from the Poisson distribution ($r_{tij}=1$), which belongs to the $k$th row and $\ell$th column cluster over time.

From Equations (26) and (27), when we compare $\lambda_{k,\ell}$ in the dPIRM and dZIPIRM, the following theorem is derived.

Theorem 1**.**

[TABLE]

where

[TABLE]

and then $\alpha_{k\ell}\in[0,1]$.

The proof of Theorem 1 is presented in Appendix.

From Equation (28) we obtain an intuitive interpretation of $\lambda_{k\ell}$. Ignoring $b$, $\alpha_{k\ell}$ represents the proportion of data that belongs to the $k$th row and $\ell$th column cluster and that are generated from a Poisson distribution ($r_{tij}=1$). That is, in the dZIPIRM, the more data in a cluster generated from the Degenerate distribution, the larger $\lambda_{k\ell}^{(dZIP)}$ tends to be compared with $\lambda_{k\ell}^{(dP)}$. Conversely, the more data generated from the Poisson distribution, the closer $\lambda_{k\ell}^{(dZIP)}$ tends to be to $\lambda_{k\ell}^{(dP)}$. In short, under the estimation framework using the full conditional distribution, the estimated $\lambda_{k\ell}^{(dZIP)}$ tends to become larger than $\lambda_{k\ell}^{(dP)}$ because in the dZIPIRM, the zero data, which are considered to be generated from the Degenerate distribution, are not taken into account to estimate $\lambda_{k\ell}^{(dZIP)}$. This theoretical result is consistent with the model interpretation in Section 2.

From Equation (29), we may also say that the dispersion of $\lambda_{k\ell}^{(dZIP)}$ tends to be larger than that of $\lambda_{k\ell}^{(dP)}$. However, as shown in the simulation study in Section 4, this result does not necessarily affect the clustering result.

Next, we consider $w_{tij}$ in Equation (25). Again, the expected value and variance of $w_{tij}$ is expressed as

[TABLE]

From Equation (30), we can see how the hyperparameters $c$ and $d$ affect the estimation. Especially, if $c$ and $d$ are set to be equal and small, the estimation of $w_{tij}$ is less affected by the hyperparameters $c$ and $d$.

4 Simulation Study

In this section, we conduct a simulation to evaluate the fit to the data and clustering accuracy of the dPIRM and dZIPIRM.

4.1 Data generation and setting

We prepared several synthetic datasets. To create these datasets, we first stipulated that the number of clusters is $K=4$, the number of time steps is $T=5$, and the number of objects is $N=16$. Next, we generated $\lambda_{k\ell}$ from the discrete uniform distribution $\mathrm{Uniform}\{0,9\}$ ($t=1,\cdots,T$, $k,\ell=1,\cdots,K$), and then $x_{tij}$ from $\mathrm{Poisson}(\lambda_{z_{ti}z_{tj}})$ ($i,j=1,\cdots,N$, $z_{ti}=1,\cdots,K$), as shown in Figure 2. When the time changes from $t$ to $t+1$, we randomly move objects to other clusters in three patterns of movement ratio $m$, where $m=0.1,0.2,0.3$. Finally, we randomly convert the numerical value of the data into zero with a zero ratio $s$, where $s=0.3,0.5,0.7$. We thus generate data for all nine ($=3\times 3$) settings.

To estimate the parameters, we need to prespecify the hyperparameters. In Section 3.2, we showed the effect of the hyperparameters on the expected value. In this simulation, since there is no prior information, we reduce the influence of the hyperparameters as much as possible. Therefore, to reduce the influence of hyperparameters, we set $a=b=1$. Similarly, we set $c=d=1$ as well as $\kappa=\gamma=\alpha_{0}=1$.

We compare the PIRM, dPIRM, and dZIPIRM in this scenario by sampling the parameters 300 times for each method. We use two measurements for the evaluation. One is the Rand index (RI), which computes the similarity between the true and estimated clustering results [13, 7]. The RI is coded one (zero) when two clustering results match (do not match). As with the computation of the RI in the dIRM, we compute the RI between the assignment of correct answers $Z_{t}$ and the estimated $\hat{Z}_{t}$ for each time step, and then average the RI for $T$ time points. For the evaluation of cluster allocation in PIRM, where one cluster assignment is obtained in each time step (and thus T cluster assignments in total), we calculate $T$ RIs between the cluster assignment obtained by PIRM and true cluster assignment $Z_{t}$ at each time point, and pick up the highest RI result of all $T$ RI results. The second measurement is log-likelihood, for which a large value indicates that the model fits the data well [7].

We calculate the mean and standard deviation for the two measurements above for the data generated for 50 data points in each of the nine settings described earlier.

4.2 Results

Tables 1 and 2 show the results for the RI and log-likelihood, respectively. $m$ represents the rate at which the object moved between clusters at each time step and $s$ represents the ratio of objects converted into zeros.

As shown in Table 1, the RIs of the dPIRM and dZIPIRM are higher than that of the PIRM in every setting. This finding means that models incorporating time structures such as the dPIRM and dZIPIRM can cluster data better than those that cluster data at each time step, the PIRM. When comparing the dPIRM and dZIPIRM, the dZIPIRM has a higher RI than the dPIRM in several settings. Furthermore, the dZIPIRM often yields a lower standard deviation than the dPIRM. This finding means that the dZIPIRM can perform stable clustering with accurate estimates even when many zero values exist in the data.

Table 2 shows that the dPIRM and dZIPIRM yield higher log-likelihood values than the PIRM in any setting, just as with the RI. This means that we can better describe the assumed scenario by employing models incorporating time structures such as the dPIRM and dZIPIRM. When $m=0.1$, the dZIPIRM yields high log-likelihood values for every zero ratio. This finding shows that the dZIPIRM fits the data well when considering zero-inflated data for objects unlikely to move at each time step. Further, the dZIPIRM has higher log-likelihood values than the dPIRM for many patterns, like the RI, while the standard deviation is similar between the dPIRM and dZIPIRM, unlike the RI.

These results suggest that both the dPIRM and the dZIPIRM are superior to the PIRM in terms of RI and log-likelihood for this simulation setting. In addition, the dZIPIRM provides an accurate and stable clustering allocation as well as a good data fit compared with the dPIRM when time structures and zero-inflated data are assumed.

5 Application

In this section, we illustrate how the dPIRM and dZIPIRM may be applied to POS data to identify a group of products likely to be purchased simultaneously.

5.1 Data and setting

We analyze the POS data of a fashion brand in Japan, provided by Sothink Co., Ltd. These data provide information for two store types: shop-in-shop and franchise. Since fashion brand products are sensitive to the influence of the season in Japan, it is also important to consider marketing strategy changes depending on the season.

The data cover 78,807 customers and 12,937 products. Product categories include short pants, small items, and short-sleeved T-shirts. We view product categories as objects in the model and count simultaneous purchases by product category.

We formulate a matrix showing the number of simultaneous purchases of 40 types of objects every month from January to June 2016 ($T=6$, $N=40$). That is, we create an $N\times N$ relational count data matrix for $T$ time steps. The zero-element ratio for the entire matrix is 0.52. We analyze the resulting dataset using the dPIRM and dZIPIRM. Since there is no prior hyperparameter information, we reuse the simulation approach.

5.2 Results

In this section, we first compare the log-likelihood values of the PIRM, dPIRM, and dZIPIRM to assess the fit of each model.

Table 3 shows that the dZIPIRM has the highest log-likelihood, meaning that it provides the better fit to the data. Hence, assuming a ZIP and a time-varying structure improves the overall fit.

Next, we show the results for the dPIRM. The number of clusters in the dPIRM is estimated to be six. Figure 3 shows the estimated $\lambda_{k\ell}$ corresponding to each cluster combination, while Figure 4 shows the total number of objects in each cluster for the six-month sample period. Table 4 shows the objects belonging to the corresponding clusters. Since the parameter $\lambda_{k,\ell}$ represents the strength of the relationship, if the value of the estimated $\lambda_{k\ell}$ is large, the objects in Cluster $k$ and Cluster $\ell$ are more likely to be purchased simultaneously.

Figure 3 shows that the objects in Cluster 1 are more likely to be purchased with those in Cluster 1 simultaneously (i.e., $\hat{\lambda}_{11}$ is high). Further, the objects in Cluster 3 are more likely to be purchased simultaneously with those in Cluster 1 and Cluster 5 (i.e., $\hat{\lambda}_{31}$ and $\hat{\lambda}_{35}$ are high), while they are unlikely to be purchased with other items in Cluster 3 (i.e., $\hat{\lambda}_{33}$ is low). The items in Cluster 5 are more likely to be bought with items in Cluster 1, followed in order of those in Clusters 5 and 3 (i.e., $\hat{\lambda}_{51}>\hat{\lambda}_{55}>\hat{\lambda}_{53}$).

From the number of objects belonging to each cluster in Figure 4, we see that the number in Cluster 1 is much larger than that for the other clusters. Therefore, the objects in Cluster 1 are more likely to be purchased simultaneously with several other types of objects. As shown in Table 4, there are two types of objects in Cluster 1: objects that remain in Cluster 1 for a long period and items that move from or to Cluster 1 depending on the season.

Examples of objects that remain in Cluster 1 for a long period include small items such as handkerchiefs and necklaces, which may be interpreted as objects likely to be purchased along with other items throughout the year. We may also interpret items such as ties as being more likely to be bought with a suit throughout the year.

Meanwhile, items that move from or to Cluster 1 depending on the season, such as slacks and short-sleeved T-shirts, belong to Clusters 3 and 5 at the beginning of the year, but move to Cluster 1 from April onwards. This date is the turning point of the season in Japan, when spring starts and temperatures rise. This means that slacks and short-sleeved T-shirts are not purchased much in the winter, whereas in the spring they are more likely to be bought along with other objects.

We now show the results for the dZIPIRM. The estimated number of clusters is 13. Figures 5, 6 and Table 5 for the dZIPIRM follow Figures 3, 4 and Table 4 for the dPIRM. Figure 5 shows that the number of clusters in the dZIPIRM is larger than that in the dPIRM because of the increase in cluster expressiveness, which in turn is a consequence of the additional flexibility provided by the assumption that there are many zeros.

Beyond the estimated number of clusters, the results for the dZIPIRM are similar as those for the dPIRM. For example, Cluster 1 of the dZIPIRM is similar to Cluster 1 of the dPIRM (e.g., always including small items, ties, and suits), while some objects move from cluster to cluster depending on the season. However, using the dZIPIRM, we may interpret the result in more detail. For example, the dPIRM could not distinguish between the behavior of slacks and short-sleeved T-shirts, whereas the dZIPIRM can interpret the difference between these objects. Specifically, in January, the category “slacks” belongs to Cluster 12, while the category “short-sleeved T-shirts” belongs to Cluster 10. Products in Cluster 12 are unlikely to be purchased with objects in any cluster. Meanwhile, objects in Cluster 10 are likely to be purchased with objects in Cluster 2 (e.g., casual shirts and jeans). Thus, by assuming a zero-inflated setting using the dZIPIRM, we obtain a result that is easier to interpret.

From the above analysis, we may conclude that since the dPIRM yields fewer clusters, only rough trends may be modeled. However, as the dZIPIRM yields more clusters, a more detailed interpretation of the clusters may be obtained. As Table 3 shows, the dZIPIRM fits the data better than the dPIRM. Despite this, as is typically the case with empirical data, it is difficult to objectively validate which method should be preferred. Nevertheless, this empirical application shows that both the dPIRM and the dZIPIRM allow us to see how customer buying behavior changes with the seasons. In addition, by using the dZIPIRM, we obtain a clustering result that provides a more detailed interpretation of the clusters when the data contain many zeros.

6 Conclusion and Discussion

We proposed two bi-clustering methods for count data considering changes over time. The first method, the dPIRM, enabled the bi-clustering of time-varying count data by using the time structure and distribution information for these data. The second method, the dZIPIRM, extended the dPIRM to express zero-inflated real data. To derive the posterior distribution and clarify the properties of the parameters, zero-inflated data were formulated by using latent variables.

We derived the explicit posterior distributions of the parameters and outlined their properties under the estimation framework using a full conditional distribution. Specifically, by comparing the expected value and variance of $\lambda$ for both the dPIRM and the dZIPIRM, we showed that how the specific model assumptions affect the estimation result. We also explained the theoretical relationship between the two models. Although this theorem is specifically derived to compare our two proposed model, dPIRM and dZIPIRM, this could be also applied for other ZIP-related model.

By incorporating a time structure, we further showed that both the dPIRM and the dZIPIRM outperform the PIRM in terms of the RI and log-likelihood values in a simulation study. We also observed that the dZIPIRM yields more stable and better clustering allocations as well as a better data fit than the dPIRM.

For the empirical data, the dZIPIRM, which assumes the properties of the ZIP and a time-varying structure, resulted in the best data fit. Moreover, the clustering result of the dZIPIRM is similar to that of the dPIRM. However, using the dZIPIRM allows us to interpret the cluster structure in more detail, as this model relies on more flexible assumptions than the dPIRM.

We consider two future works. At first, although the row and column objects in the dPIRM and dZIPIRM may differ, for comprehensibility, we analyzed the data using the same object for the rows and columns. Therefore, future work will focus on real data in situations where the objects of the rows and columns are different.

Future research will also study how the choice of hyperparameters affects the estimation results in a simulation setting (i.e., how the estimation result changes depending on the choice of hyperparameters). In Section 3.2, we referred to the relationship between the hyperparameters and parameters, and we believe further work on this area is needed to provide users with better guidance about hyperparameter choice.

Appendix

In this Appendix, we show proofs for Equation (23) and Theorem 1. First, Equation (23) can be derived as follows.

Proof.

[TABLE]

∎

Next, proof of Theorem 1 is derived as follows.

Proof.

[TABLE]

Here, we derive Eq. (28)

[TABLE]

Here, we derive Eq. (29)

[TABLE]

[TABLE]

∎