Calibrating covariate informed product partition models

Abstract

Covariate informed product partition models incorporate the intuitively appealing notion that individuals or units with similar covariate values a priori have a higher probability of co-clustering than those with dissimilar covariate values. These methods have been shown to perform well if the number of covariates is relatively small. However, as the number of covariates increase, their influence on partition probabilities overwhelm any information the response may provide in clustering and often encourage partitions with either a large number of singleton clusters or one large cluster resulting in poor model fit and poor out-of-sample prediction. This same phenomenon is observed in Bayesian nonparametric regression methods that induce a conditional distribution for the response given covariates through a joint model. In light of this, we propose two methods that calibrate the covariate-dependent partition model by capping the influence that covariates have on partition probabilities. We demonstrate the new methods’ utility using simulation and two publicly available datasets.

Appendices

Appendix: MCMC algorithm for the calibrated similarity and tempered mixture of experts

Here we provide pertinent computation details for the MCMC algorithm used to fit the TME model and PPMx with calibrated similarity. We focus primarily on the updating of cluster labels, as conditional on these updating the remaining model parameters is straightforward employing a Gibbs sampler or Metropolis–Hastings steps.

1.1 Calibrated similarity

To update the cluster membership of subject i for the calibrated similarity, cluster weights are created by comparing the unnormalized posterior for the jth cluster when subject i is excluded from that when subject i is included. In addition to weights for existing clusters, algorithm 8 of Neal (2000) requires calculating weights for p empty clusters whose cluster-specific parameters are auxiliary variables generated from the prior. To make this more concrete, let $S_j^{-i}$ denote the jth cluster and $k^{-i}$ the number of clusters when subject i is not considered. Similarly $\varvec{x}_j^{\star -i}$ will denote the vector of covariates corresponding to cluster h when subject i has been removed. Then the multinomial weights associated with the $k^{-i}$ existing clusters and one empty cluster are

$$\begin{aligned}&Pr(s_i = j | - ) \propto \end{aligned}$$

(16)

$$\begin{aligned}&{\left\{ \begin{array}{ll} N(y_i ; \mu ^{\star }_{j}, \sigma ^{2\star }_{j}) \displaystyle \frac{c(S_{j}^{-i}\cup \{i\})\tilde{g}(\varvec{x}^{\star -i}_{j}\cup \{\varvec{x}_i\})}{c(S_{j}^{-i})\tilde{g}(\varvec{x}^{\star -i}_{j})}\ \quad \text{ for } \ j = 1, \ldots , k^{-i}\\ N(y_i; \mu ^{\star }_{\mathrm{new}, j}, \sigma ^{2\star }_{\mathrm{new},j}) c(\{i\}) \tilde{g}(\{\varvec{x}_i\}) p^{-1}\ \quad \text{ for } \ j = k^{-i}+1. \end{array}\right. } \end{aligned}$$

(17)

where as mentioned $\mu ^{\star }_{\mathrm{new}, j}$ and $\sigma ^{2\star }_{\mathrm{new},j}$ are auxiliary variables that are drawn from their respective prior distributions. Since $\tilde{g}(\{\varvec{x}_i\})$ needs to be accounted for when standardizing the multinomial weights, we employ the following ratios in the MCMC algorithm

$$\begin{aligned} \tilde{g}(\varvec{x}^{\star -i}_{j} \cup \varvec{x}_i)&= \frac{g(\varvec{x}^{\star -i}_{j} \cup \varvec{x}_i)}{\sum _{\ell }g(\varvec{x}^{\star -i}_{\ell } \cup \varvec{x}_i) } \\ \tilde{g}(\varvec{x}^{\star -i}_{j})&= \frac{g(\varvec{x}^{\star -i}_{j})}{\sum _{\ell }g(\varvec{x}^{\star -i}_{\ell })+ g(\{\varvec{x}_i\})}\\ \tilde{g}(\{\varvec{x}_i\})&= \frac{g(\{\varvec{x}_i\})}{\sum _{\ell }g(\varvec{x}^{\star -i}_{\ell })+ g(\{\varvec{x}_i\})}.\\ \end{aligned}$$

When $\varvec{x}_i$ is included in the j cluster then it is not able to form its own singleton. However, when it is excluded from the jth cluster, then it is completely plausible that it forms its own singleton cluster. For these reasons the similarity value $g(\{\varvec{x}_i\})$ is only included in $\tilde{g}(\varvec{x}^{\star -i}_{j})$ and $\tilde{g}(\{\varvec{x}_i\})$.

1.2 TME with unknown $\varvec{\xi }^{\star }_j$ and fixed J

Upon introducing latent component labels $s_i$ such that $Pr(s_i = j) = w(\varvec{x}_i; \varvec{\xi }^{\star }_j)$, the data model (14) can be written hierarchically as

$$\begin{aligned} p(\varvec{y} | \varvec{x}, \varvec{\mu }^{\star }, \varvec{\sigma }^{2\star }, \varvec{\xi }^{\star }, \varvec{c})&= \prod _{i=1}^m\prod _{\ell =1}^J N(y_i|\mu ^{\star }_{\ell }, \sigma ^{2\star }_{\ell })^{I[s_i = \ell ]}\end{aligned}$$

(18)

$$\begin{aligned} s_i&\sim \sum _{\ell =1}^J \delta _{\ell }w(\varvec{x}_i; \varvec{\xi }^{\star }_{\ell }) \end{aligned}$$

(19)

where $\delta _{\ell }$ is the dirac measure. With this hierarchical representation, a MCMC algorithm can be constructed by cycling through the following

• Update component labels using

  $$\begin{aligned} Pr(s_i = h | - )&\propto N(y_i | \mu ^{\star }_h, \sigma ^{2\star }_h) w(\varvec{x}_i ; \varvec{\xi }^{\star }_h) \end{aligned}$$

• If $\varvec{x}_i$ is comprised of continuous and categorical variables, then without loss of generality let $\varvec{x}_i = (x_{1i}, x_{2i})$ where $x_{1i}$ is continuous and $x_{2i}$ is categorical. Further, $\varvec{\xi }_j^{\star } = (\eta ^{\star }_j,v^{2\star }_j, \varvec{\pi }^{\star }_j)$ with $\varvec{\xi }^{\star } = (\varvec{\xi }_1^{\star },\ldots , \varvec{\xi }_J^{\star })$. Then $\varvec{\xi }_j^{\star }=(\eta ^{\star }_j,v^{2\star }_j, \varvec{\pi }^{\star }_j)$ can be updated within the MCMC algorithm by way of a Metroplis–Hastings step employing

  $$\begin{aligned} {[}\varvec{\xi }^{\star }_j | - ]&\propto \prod _{i=1}^m Pr(s_i | \varvec{\xi }^{\star })\prod _{j=1}^J p(\varvec{\xi }^{\star }_j) \\&\propto \prod _{i=1}^m w(\varvec{x}_{i}; \varvec{\xi }^{\star }_1)^{I[s_i=1]} \times \ldots \times w(\varvec{x}_{i}; \varvec{\xi }^{\star }_1)^{I[s_i=J]} p(\varvec{\xi }^{\star }_j) \\&\propto \prod _{i:s_i=j} w(\varvec{x}_{i}; \varvec{\xi }^{\star }_j) p(\varvec{\xi }^{\star }_j)\\&= \prod _{i:s_i=j} \frac{q(x_{i1}|\eta ^{\star }_j, v^{2\star }_j) q(x_{i2}|\varvec{\pi }^{\star }_j)}{\sum _{\ell }^Jq (x_{i1}|\eta ^{\star }_{\ell }, v^{2\star }_{\ell })q(x_{i2}|\varvec{\pi }^{\star }_{\ell })}p(\eta ^{\star }_j,v^{2\star }_j, \varvec{\pi }^{\star }_j) \end{aligned}$$

  where $q(x_{i1}|\eta ^{\star }_j, v^{2\star }_j)$ is normal density and and $q(x_{i2}|\varvec{\pi }^{\star }_j)$ a multinomial density. For $\varvec{\pi }^{\star }_j$ an independent Metropolis–Hastings sampler with uniform (over the simplex) candidate density may be considered. This candidate density will cancel out in the Metropolis–Hastings ratio (though this may become more inefficient as the number of categories in $x_{i2}$ increases). Updating $\eta ^{\star }_j$ and $v^{2\star }_j$ can be accomplished using a random walk Metropolis step with normal candidate density for both.

• Updating the likelihood parameters $\mu _j^{\star }$ and $\sigma _j^{2\star }$ can be carried out using Gibbs steps as their full conditionals have well known closed forms.

Computing Gower’s dissimilarity

The daisy function found in the cluster package (Maechler et al. 2016) of the statistical software R was employed to calculate the Gower dissimilarity. The calculated dissimilarity is an “average” of the individual p dissimilarities

$$\begin{aligned} d(\varvec{x}_i, \varvec{x}_j) = \frac{1}{p}\sum _{\ell =1}^p d(x_{i\ell },x_{j\ell }). \end{aligned}$$

For numeric or continuous x’s, $d(x_{i\ell },x_{j\ell }) = |x_{i\ell } - x_{j\ell }|/R_{\ell }$ where $R_{\ell } = \max _h(x_{h\ell }) - \min _h(x_{h\ell })$. For nominal variables

$$\begin{aligned} d(x_{i\ell },x_{j\ell }) =\left\{ \begin{array}{cl} 0 &{} \quad \text{ if } \,\,x_{i\ell }=x_{j\ell } \\ 1 &{} \quad \text{ otherwise }. \end{array} \right. \end{aligned}$$