Robust Finite Mixture Regression for Heterogeneous Targets

Let $\mathbf{Y}\in\mathbb{R}^{n\times m}$ be the output/target data matrix and $\mathbf{X}\in\mathbb{R}^{n\times d}$ the input/feature data matrix, consisting of $n$ independent samples $(\mathbf{y}_{i},\mathbf{x}_{i})$, $i=1,\ldots,n$. As such, there are $m$ different targets with a common set of $d$ features. We allow $\mathbf{Y}$ to contain missing values at random; define $\Omega_{i}=\{j\in\{1,\ldots,m\}:y_{ij}\mbox{ is observed.}\}$ be the collection of indices of observed outcome in the $i$th sample $\mathbf{y}_{i}$ ($\Omega_{i}\neq\emptyset$), for $i=1,\ldots,n$.

To model multiple types of targets, such as continuous, binary, count, etc., we allow $y_{ij}$ to potentially follow different distributions in the exponential-dispersion family, for each $j=1,\ldots,m$. Specifically, we assume that given $\mathbf{x}_{i}$, the joint probability density function of $\tilde{\mathbf{y}}_{i}=\{y_{ij};j\in\Omega_{i}\}$ is

where

$\varphi_{ijr}$ is the natural parameter for the $i$th sample of the $j$th target in the $r$th mixture component, $\phi_{jr}$ is the dispersion parameter of the $j$th target in the $r$th mixture component, and the functions $a_{j},b_{j},c_{j}\ (j=1,\ldots,m)$ are determined by the specific distribution of the $j$th target. Here, the key assumption is that the $m$ tasks all correspond to the same cluster structure (e.g., the $m$ tasks all have $k$ clusters) determined by the underlying population heterogeneity; given the shared cluster label (e.g., $r$), the tasks within each mixture component then become independent of each other (depicted by the product of their probability density functions). As such, by allowing cluster label sharing, the model provides an effective way to genuinely integrate the learning of the multiple tasks.

Following the setup of GLMs, we assume a linear structure in the natural parameters, i.e.,

where $\mbox{\boldmath$\beta$}_{jr}$ is the regression coefficient vector of the $j$th response in the $r$th mixture component. Since $\mathbf{x}_{i}$ is possibly of high dimensionality, the $\mbox{\boldmath$\beta$}_{jr}$s are potentially sparse vectors. For example, when the $\mbox{\boldmath$\beta$}_{jr}$s for $j=1,\ldots,m$ share the same sparsity pattern, the tasks share the same set of relevant features within each mixture component. For $r=1,\ldots,k$, write $\mbox{\boldmath$\beta$}_{r}\in\mathbb{R}^{d\times m}=[\mbox{\boldmath$\beta$}_{1r},\mbox{\boldmath$\beta$}_{2r},\ldots,\mbox{\boldmath$\beta$}_{mr}]$ and $\phi_{r}=[\phi_{1r},\ldots,\phi_{mr}]^{T}$. Also write $\mbox{\boldmath$\beta$}\in\mathbb{R}^{(d\times m)\times k}=[\mbox{\boldmath$\beta$}_{1},\ldots,\mbox{\boldmath$\beta$}_{k}]$. Let $\theta=\{\mbox{\boldmath$\beta$},\phi_{1},\ldots,\phi_{k},\pi_{1},\ldots,\pi_{k}\}$ collecting all the unknown parameters, with the parameter space given by $\Theta=\mathbb{R}^{(d\times m)\times k}\times\mathbb{R}^{m\times k}_{>0}\times\Pi,$ where $\Pi=\{\pi;\pi_{r}>0\mbox{ for }r=1,\ldots,k\mbox{ and }\sum_{r=1}^{k}\pi_{r}=1\}$.

The data log-likelihood of the proposed model is

The missing values in $\mathbf{Y}$ simply do not contribute to the likelihood, which follows the same spirit as in matrix completion (Candès and Recht, 2009). The proposed model indeed possesses a genuine multivariate flavor, as the different outcomes share the same underlying latent cluster pattern of the heterogeneous population. We then propose to estimate $\theta$ by the following penalized likelihood criterion:

where $\mathcal{R}(\mbox{\boldmath$\beta$};\lambda)$ is some certain penalty term on the regression coefficients with $\lambda$ being a tuning parameter.

We thus name our proposed method the HEterogeneous-target Robust MIxTure regression (Hermit). The $\mathcal{R}(\mbox{\boldmath$\beta$};\lambda)$ can be flexibly chosen based on specific needs of feature selection. The first sparse penalties adopted by our model is the $\ell_{1}$ norm (lasso-type) penalty,

where $\lambda$ is the tuning parameter, $\|\cdot\|_{1}$ is the entry-wise $\ell_{1}$ norm, and $\pi_{r}^{\gamma}$s $(r=1,\ldots,k)$ are the penalty weights with $\gamma\in\{0,1/2,1\}$ being a pre-specified constant. Here the penalty also depends on the unknown mixture proportions $\pi$; when the cluster sizes are expected to be imbalanced, using this weighted penalization with some $\gamma>0$ is preferred (Städler et al, 2010). This entry-wise regularization approach allows the tasks to have independent set of relevant features. Alternatively, in order to enhance the integrative learning and potentially boost the performance of clustering, it could be beneficial to encourage the internal similarity within each sub-population. Then certain group-wise regularization of the features could be considered, which are widely adopted in multi-task learning. In particular, we consider a component-specific group sparsity pattern to achieve shared feature selection among different tasks, in which the group $\ell_{1}$ norm penalty is used (Gong et al, 2012a; Jalali et al, 2010),

where $\|\cdot\|_{1,2}$ denotes the sum of the row $\ell_{2}$ norms of the enclosed matrix, and the weights are constructed as in (5). The shared feature set in each sub-population can be used to characterize the sub-population and render the whole model more interpretable.

We propose a generalized EM (GEM) algorithm (Dempster et al, 1977) to solve the minimization problem in (4). For each $i=1,\ldots,n$, define $(\delta_{i,1},\ldots,\delta_{i,k})$ be a set of latent indicator variables, where $\delta_{i,r}=1$ if the $i$th sample $(\mathbf{y}_{i},\mathbf{x}_{i})$ belongs to the $r$th component of the mixture model (1) and $\delta_{i,r}=0$ otherwise. So $\sum_{r=1}^{k}\delta_{i,r}=1$, $\forall i$. These indicators are not observed since the cluster labels of the samples are unknown. Let $\delta$ denote the collection of all the indicator variables. By treating $\delta$ as missing, the EM algorithm proceeds by iteratively optimizing the conditional expectation of the complete log-likelihood criterion.

The complete log-likelihood is given by

where $\varphi_{ijr}=\mathbf{x}_{i}\mbox{\boldmath$\beta$}_{jr}$, for $i=1,\ldots,n,$ $j=1,\ldots,m$, and $r=1,\ldots,k$. The conditional expectation of the penalized complete negative log-likelihood is then given by

where $\mathcal{R}(\mbox{\boldmath$\beta$};\lambda,\pi)$ can be any of the penalties in (5) or (6). It is easy to show that deriving $Q_{pen}(\theta\mid\theta^{\prime})$ boils down to the computation of $\mathbb{E}[\delta_{i,r}\mid\mathbf{Y},\mathbf{X},\theta^{\prime}]$, which admits an explicit form.

The EM algorithm proceeds as follows. Let $\theta^{(0)}$ be some given initial values. We repeat the following steps for $t=0,1,2,\ldots$, until convergence of the parameters or the pre-specified maximum number of iteration $T_{out}$ is reached.

a) Update $\pi=(\pi_{1},\ldots,\pi_{k})$ by solving

b) Update $\mbox{\boldmath$\beta$},\mbox{\boldmath$\Phi$}$.

For the problem in a), Städler et al (2010) proposed a procedure to lower the objective function by a feasible point, and we find that simply setting $\pi_{r}^{(t+1)}=\sum_{i=1}^{n}\hat{\rho}_{i,r}/n$ is good enough. For the problem in b), we use an accelerated proximal gradient (APG) method introduced in Nesterov et al (2007) with the maximum number of iteration of $T_{in}$. The update steps by proximal operators correspond to the chosen penalty form. For the entry-wise $\ell_{1}$ norm penalty in (5),

where $\circ$ denotes entry-wise product, $\widetilde{\mbox{\boldmath$\beta$}}_{r}^{(t)}=\mbox{\boldmath$\beta$}_{r}^{(t)}+\tau\triangle\mbox{\boldmath$\beta$}_{r}^{(t)}$, $\tau$ denotes the step size, and $\triangle\mbox{\boldmath$\beta$}_{r}^{(t)}$ denotes the update direction of $\mbox{\boldmath$\beta$}_{r}^{(t)}$ determined by APG. For the group $\ell_{1}$ norm penalty in (6),

where $\mbox{\boldmath$\beta$}_{r,j}$ denotes the $j$th column of $\mbox{\boldmath$\beta$}_{r}$. We adopt the active set algorithm in Städler et al (2010) to speed up the computation.

The time complexity of our algorithm using the speed up technique is $O(T_{out}kmnsT_{in})$ with $s$ being the number of non-zero parameters. The algorithm performs well in practice, and we have not observed any convergence issues in our extensive numerical studies.

From the model estimation, we can get estimates of both the conditional probabilities $p(\delta_{i,r}=1\mid y_{ij^{\prime}},j^{\prime}\in\Omega_{i},\mathbf{x}_{i},\theta)$ and the conditional means $\mathbb{E}[Y_{ij}\mid\mathbf{x}_{i},\theta,\delta_{i,r}=1]$, where $\mathbb{E}[Y_{ij}\mid\mathbf{x}_{i},\theta,\delta_{i,r}=1]=\mu_{j}(\varphi_{ijr})=b^{\prime}_{j}(\mathbf{x}_{i}\mbox{\boldmath$\beta$}_{jr})$. Specifically, the conditional probabilities can be estimated by $p(\delta_{i,r}=1\mid y_{ij^{\prime}},j^{\prime}\in\Omega_{i},\mathbf{x}_{i},\hat{\theta})$ which corresponds to (7), taking $t=T_{out}$. The conditional expectations can be estimated as $\mathbb{E}[Y_{ij}\mid\mathbf{x}_{i},\hat{\theta},\delta_{i,r}=1]=b^{\prime}_{j}(\mathbf{x}_{i}\hat{\mbox{\boldmath$\beta$}}_{jr})$.

For clustering the samples, we adopt the Bayes rule, i.e., for $i=1,\ldots,n$,

Following the idea of Jacobs et al (1991), we propose to make imputation for the missing outcomes by

Unless otherwise specified, all the hyper-parameters, including regularization coefficients $\lambda$s and the number of clusters $k$, are tuned to maximize the data log-likelihood in (3) on the held-out validation data set. In other words, we fit models on training data with different specific hyper-parameter settings, and then the optimal model is chosen as the one that gives the highest log-likelihood in (3) of the held-out validation data set. This approach is fairly standard and has been widely used in existing works (Städler et al, 2010). Moreover, cross validation and various information criteria (Bhat and Kumar, 2010; Aho et al, 2014) can also be applied to determine hyper-parameters.

To perform robust estimation for parameters in the presence of outlier samples, we propose to adopt the mean shift penalization approach (She and Owen, 2011). Specifically, we extend the natural parameter model to the following additive form,

where $\zeta_{ijr}$ is a case-specific mean shift parameter to capture the potential deviation from the linear model. Apparently, when $\zeta_{ijr}$ is allowed to vary without any constraint, it can make the model fit as perfect as possible for every $y_{ijr}$. The merit of this approach is realized by assuming certain sparsity structure of the $\zeta_{ijr}$s, so that only a few of them have nonzero values corresponding to anomalies. Write $\mbox{\boldmath$\zeta$}_{r}\in\mathbb{R}^{n\times m}=[\mbox{\boldmath$\zeta$}_{1r},\mbox{\boldmath$\zeta$}_{2r},\ldots,\mbox{\boldmath$\zeta$}_{mr}]$ for $r=1,\ldots,k$, and $\mbox{\boldmath$\zeta$}\in\mathbb{R}^{(n\times m)\times k}=[\mbox{\boldmath$\zeta$}_{1},\ldots,\mbox{\boldmath$\zeta$}_{r}]$. We can then conduct joint model estimation and outlier detection by extending (4) to

where, for example, the penalty on $\zeta$ can be chosen as the group $\ell_{1}$ penalty,

so that entries of $\zeta$ are nonzero for only a few data samples.

The proposed GEM algorithm can be readily extended to handle the inclusion of $\zeta$, for which we omit the details.

Besides outlier samples, certain tasks, referred to as anomaly tasks, may not follow the assumed shared structure and thus can ruin the overall model performance. To handle anomaly tasks, though it is also intuitive to adopt the approach above, our numerical study suggests that its performance is sensitive to the tuning parameters. Here, we adopt the idea of Koller and Sahami (1996) , by utilizing the estimated conditional probabilities to measure how well a task is concordant with the estimated mixture structure. Consider the $h$th task. The main idea is to measure the discrepancy between $p(\delta_{ir}=1\mid y_{ij},j\in\Omega_{i},\mathbf{x}_{i},\theta)$, the conditional probability based on data from all observed targets, and $p(\delta_{ir}=1\mid y_{ih},\mathbf{x}_{i},\theta)$, the conditional probability based on only the $h$th task. If $h$th task is an anomaly task, it is expected that the two conditional probabilities would differ more from each other (Koller and Sahami, 1996; Law et al, 2002).

For $r=1,\ldots,k$, $i=1,\ldots,n$, let

Define $\mathbf{P}_{\Omega}=[P_{\Omega,ir}]_{n\times k}$ and $\mathbf{P}_{h}=[P_{h,ir}]_{n\times k}$. Then we define the concordant score of the $h$th task as

where $D_{KL}$ is the widely used Kullback-Leibler divergence (Cover and Thomas, 2012).

The tasks can then be ranked based on their concordant scores. As such, the detection of anomaly tasks boils down to a one-dimensional outlier detection problem. After anomaly tasks are detected, their FMR models can be built.

In practice, tasks may be clustered into groups such that each task group has its own model structure. Here we assume that each cluster of tasks shares a FMR structure defined in (1), and propose to construct a similarity matrix to discover the potential cluster pattern among tasks.

We consider a two-stage strategy. First, each task learns a FMR model on the training data independently with the same pre-fixed $k$. Then we get $\mathbf{P}_{h}=[P_{h,ir}]_{n\times k}$ for all $h=1,\ldots,m$, where

and $\delta_{h,ir}(i=1,\dots,n,r=1,\dots,k)$ and $\hat{\pi}_{hr}(r=1,\dots,k)$ are the latent variables and the estimated prior probabilities of the $h$th task, respectively.

Second, we adopt Normalized Mutual Information (NMI) (Strehl and Ghosh, 2002a) to measure the similarity between each pair of tasks. We choose NMI instead of Kullback-Leibler divergence because NMI can handle the case that the orders of clusters of two $k$-cluster structures are different. Specifically, given two methods to estimate latent variables, which are denoted by method $u$ and method $v$, let $\mathbf{P}_{u}=[P_{u,ir}]_{n\times k}$ and $\mathbf{P}_{v}=[P_{v,ir}]_{n\times k}$ denote the estimated probability of latent variables of method $u$ and method $v$, respectively, where $P_{u,ir}=p(\delta_{u,ir}=1),P_{v,ir}=p(\delta_{v,ir}=1)$ for $i=1,\ldots,n,r=1,\ldots,k$. NMI is defined as

where $I(\mathbf{P}_{u},\mathbf{P}_{v})$ denotes the mutual information between $\mathbf{P}_{u},\mathbf{P}_{v}$ such that

Following Strehl and Ghosh (2002a), we approximate $p(\delta_{u,a}=1,\delta_{v,b}=1)$, $p(\delta_{u,a}=1)$ and $p(\delta_{v,b}=1)$ by

As such, given the estimated models $\hat{\theta}_{u},\hat{\theta}_{v}$ for the $u$th and the $v$th task, respectively, we treat $p(\delta_{u,ir}=1\mid y_{iu},\mathbf{x}_{i},\hat{\theta}_{u})$ and $p(\delta_{v,ir}=1\mid y_{iv},\mathbf{x}_{i},\hat{\theta}_{v})$ as $p(\delta_{u,ir}=1)$ and $p(\delta_{v,ir}=1)$, respectively, for $i=1,\ldots,n,r=1,\ldots,k$. Then NMI between each pair of tasks are computed by (16). We note that for simplicity we set the pre-fixed $k$ to be the same, but in general $k$ can be different for different tasks by the definition of Mutual Information.

Given the similarity between each pair of tasks, any similarity-based clustering method can be applied to cluster $m$ tasks into groups. Empirically, the performance of task clustering is not sensitive to the pre-fixed $k$. As such, we set the pre-fixed $k$ to be $20$ in this paper. We then apply the proposed Hermit  approach separately for each task group.

In real-applications, one may require to use $\mathbf{x}_{i}$ only to infer the latent variables $\delta_{i,r}$ and then to predict $\mathbf{y}_{i}$, for $i=1,\ldots,n,r=1,\ldots,k$. Here we further extend our method following the idea of Mixture-Of-Experts (MOE) (Yuksel et al, 2012) model; the only modification is that $\pi_{r}$ in (1) is assumed to be function of $\mathbf{x}_{i}$, for $i=1,\ldots,n$.

To be specific, let $\mbox{\boldmath$\alpha$}=[\mbox{\boldmath$\alpha$}_{1},\ldots,\mbox{\boldmath$\alpha$}_{k}]\in\mathbb{R}^{d\times k}$ collect regression coefficient vectors for a multinomial linear model. We assume that given $\mathbf{x}_{i}$, the joint probability density function of $\tilde{\mathbf{y}}_{i}=\{y_{ij};j\in\Omega_{i}\}$ in (1) is replaced by

where

is referred to as the gating probability. All the other terms are defined the same as in (1).

Let $\theta_{2}=\{\mbox{\boldmath$\beta$},\phi_{1},\ldots,\phi_{k}\}$, with the parameter space $\Theta_{2}=\mathbb{R}^{(d\times m)\times k}\times\mathbb{R}^{m\times k}_{>0}.$ The data log-likelihood of the MOE model is

The model estimation is conducted by extending (4) to

where, for example, the penalty on $\alpha$ can be chosen as the lasso type penalty,

The minimization problem in (19) is also solved by GEM, for which the optimization procedure is similar to that in Section 3.2. The differences occur at the E-step:

where $f(y_{ij}\mid\varphi_{ijr}^{(t)},\phi_{jr}^{(t)})=\exp\{(y_{ij}\varphi_{ijr}^{(t)}-b_{j}(\varphi_{ijr}^{(t)}))/a_{j}(\phi_{jr}^{(t)})+c_{j}(y_{ij},\phi_{jr}^{(t)})\}$, at the optimization for $\alpha$:

and at the computation of $\pi$: $\pi_{r}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}p(\delta_{i,r}=1\mid\mathbf{x}_{i},\alpha_{r}^{(t+1)})$ for $r=1,\ldots,k$.

We firstly introduce some notations. Denote the regression parameters that are subject to regularization by $\beta=\mbox{vec}(\mbox{\boldmath$\beta$}_{1},\ldots,\mbox{\boldmath$\beta$}_{k}),\phi=\mbox{vec}(\Phi_{1},\ldots,\Phi_{k})$, where $\mbox{vec}(\cdot)$ is the vectorization operator. The other parameters in the mixture model are denoted by $\eta=\mbox{vec}(\log(\phi),\log(\pi))$, where $\log(\cdot)$ is entry-wisely applied. Denote the true parameter by $\theta_{0}=(\mbox{\boldmath$\beta$}_{0},\Phi_{0,1},\ldots,\Phi_{0,k},\pi_{0,1},\pi_{0,k-1})$ to be estimated under the FMR model defined in (1) and (2). In the sequel, we always use subscripts “$0$” to represent parameters or structures under the true model. To study sparsity recovery, denote the set of indices of non-zero entries of the true parameter by $S$. We use $\lesssim$ to indicate that the inequality holds up to some multiplicative numerical constants. To focus on the main idea, we consider the case of $\gamma=0$ in the following analysis.

We define average excess risk for fixed design points $\mathbf{x}_{1},\ldots,\mathbf{x}_{n}$ based on Kullback-Leibler divergence as