Precision Therapeutic Biomarker Identification with Application to the Cancer Genome Project

TL;DR

This paper introduces a novel multivariate regression model with overlapping clusters and an EM algorithm to identify cancer-specific therapeutic biomarkers and drug resistance patterns from the Cancer Genome Project data.

Contribution

It proposes a generalized finite mixture of multivariate regression model with a new EM algorithm to analyze complex genomic and drug response data in cancer research.

Findings

Revealed new therapeutic inter-relationships between cancers

Identified existing and novel drug behaviors

Enhanced understanding of drug resistance patterns

Abstract

Cancer cell lines have frequently been used to link drug sensitivity and resistance with genomic profiles. To capture genomic complexity in cancer, the Cancer Genome Project (CGP) (Garnett et al., 2012) screened 639 human tumor cell lines with 130 drugs ranging from known chemotherapeutic agents to experimental compounds. Questions of interest include: i) can cancer-specific therapeutic biomarkers be detected, ii) can drug resistance patterns be identified along with predictive strategies to circumvent resistance using alternate drugs, iii) can biomarkers of drug synergies be predicted ? To tackle these questions, following statistical challenges still exist: i)biomarkers cluster among the cell lines; ii) clusters can overlap (e.g. a cell line may belong to multiple clusters); iii) drugs should be modeled jointly. We introduce a multivariate regression model with a latent overlapping cluster indicator variable to address above issues. A generalized finite mixture of multivariate regression (FMMR) model in connection with the new model and a new EM algorithm for fitting are proposed. Re-analysis of the dataset sheds new light on the therapeutic inter-relationships between cancers as well existing and novel drug behaviors for the treatment and management of cancer.

Tables

Table 1. Web Table 1: Accuracy of BIC in selection of K 𝐾 K .

sample size	EM	gEMseplasso
n=150	0.32	0.56
n=450	0.22	0.88

Equations

$${\mathbf{y}}_i=\sum _{k=1}^K{\mathbf{B}}_k^T{\mathbf{x}}_i{P}_{ik}+{\epsilon }_i,i=1,\dots ,n,$$

$$({\mathbf{y}}_i\mid z_i=k,{\mathbf{x}}_i,{\mathbf{B}}_k,\mathbf{\Sigma })\sim N_q({\mathbf{B}}_k^T{\mathbf{x}}_i,\mathbf{\Sigma }).$$

$$f({\mathbf{y}}_i,z_i\mid {\mathbf{x}}_i,\mathbf{\Theta })=\prod _{k=1}^K{\left\{{\pi }_{k}f({\mathbf{y}}_i\mid z_i=k,{\mathbf{x}}_i,{\mathbf{B}}_k,\mathbf{\Sigma })\right\}}^{I(z_i=k)}.$$

$$f({\mathbf{y}}_i\mid {\mathbf{x}}_i,\mathbf{\Theta })=\sum _{k=1}^K{\pi }_{k}f({\mathbf{y}}_i\mid z_i=k,{\mathbf{x}}_i,{\mathbf{B}}_k,\mathbf{\Sigma }).$$

$$l_n(\mathbf{\Theta })=\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\log {\pi }_k+\sum _{i=1}^n\sum _{k=1}^K{z}_{ik}\log f({\mathbf{y}}_i\mid z_i=k,{\mathbf{x}}_i,{\mathbf{B}}_k,\mathbf{\Sigma }),$$

$${\tilde{l}}_n(\mathbf{\Theta })=l_n(\mathbf{\Theta })-{\rho }_n(\mathbf{\Theta }),$$

$${\rho }_n(\mathbf{\Theta })=\sum _{k=1}^K{\pi }_k{\rho }_{nk}({\mathbf{B}}_k).$$

$${\rho }_{nk}({\mathbf{B}}_k)={\lambda }_k\Vert {\mathbf{B}}_k{\Vert }_1$$

$${\rho }_{nk}({\mathbf{B}}_k)={\lambda }_{k1}\Vert {\mathbf{B}}_k{\Vert }_1+{\lambda }_{k2}\Vert {\mathbf{B}}_k{\Vert }_2^2,$$

$$T={\cup }_{s=1}^K\left\{(l_1\dots l_s):\{l_1,\dots ,l_s\}\subseteq \{1,\dots ,K\}\right\}.$$

$$P\left(z_i=(l_1\dots l_s)\right)={\pi }_{(l_1\dots l_s)}.$$

$$\sum _{(l_1\dots l_s)\in T}{\pi }_{(l_1\dots l_s)}=1,{\pi }_{(l_1\dots l_s)}\ge 0.$$

$$\left({\mathbf{y}}_i\mid z_i=(l_1\dots l_s),{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,\mathbf{\Sigma }\right)\sim N_q\left(\sum _{k=1}^s{\mathbf{B}}_{l_k}^T{\mathbf{x}}_i,\mathbf{\Sigma }\right).$$

$$f({\mathbf{y}}_i,z_i$$

$$\prod _{(l_1\dots l_s)\in T}{\left\{{\pi }_{(l_1\dots l_s)}f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,\mathbf{\Sigma }\right)\right\}}^{I(z_i=(l_1\dots l_s))}.$$

$$f({\mathbf{y}}_i\mid {\mathbf{x}}_i,\mathbf{\Theta })=\sum _{(l_1\dots l_s)\in T}{\pi }_{(l_1\dots l_s)}f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,\mathbf{\Sigma }\right).$$

$$l_n^0(\mathbf{\Theta })=\sum _{i=1}^n\log \left(\sum _{(l_1\dots l_s)\in T}{\pi }_{(l_1\dots l_s)}f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,\mathbf{\Sigma }\right)\right).$$

$${\tilde{l}}_n^0(\mathbf{\Theta })=l_n^0(\mathbf{\Theta })-\sum _{k=1}^K\sum _{\begin{array}{c}(l_1\dots l_s):\\ k\in \{l_1,\dots ,l_s\}\end{array}}{\pi }_{(l_1\dots l_s)}{\rho }_{nk}({\mathbf{B}}_k),$$

$$l_n(\mathbf{\Theta })=$$

$$+\sum _{i=1}^n\sum _{(l_1\dots l_s)\in T}{z}_{i,(l_1\dots l_s)}\log f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,\mathbf{\Sigma }\right),$$

$${\tilde{l}}_n(\mathbf{\Theta })=l_n(\mathbf{\Theta })-\sum _{k=1}^K\sum _{\begin{array}{c}(l_1\dots l_s):\\ k\in \{l_1,\dots ,l_s\}\end{array}}{\pi }_{(l_1\dots l_s)}{\rho }_{nk}({\mathbf{B}}_k).$$

$$\hat{P}(z_{i,(l_1\dots l_s)}^m=1\mid {\mathbf{y}}_i,{\mathbf{x}}_i,{\mathbf{\Theta }}^m)=\frac{{\pi }_{(l_1\dots l_s)}^{m}f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k}{\}}_{k=1}^s,{\mathbf{\Sigma }}^m\right)}{{\sum }_{(l_1^{\prime }\dots l_s^{\prime })\in T}{\pi }_{(l_1^{\prime }\dots l_s^{\prime })}^{m}f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_k^{\prime }}{\}}_{k=1}^s,{\mathbf{\Sigma }}^m\right)}.$$

$${\pi }_{(l_1\dots l_s)}^{m+1}=\frac{1}{n}\sum _{i=1}^n{z}_{i,(l_1\dots l_s)}^m,(l_1\dots l_s)\in T.$$

$${\mathbf{B}}_k^{m+1}=\arg \underset{{\mathbf{B}}_k}{\max }\sum _{i=1}^n\sum _{\begin{array}{c}(l_1\dots l_s):\\ k\in \{l_1,\dots ,l_s\}\end{array}}{z}_{i,(l_1\dots l_s)}^m\log f\left({\mathbf{y}}_i\mid z_i,{\mathbf{x}}_i,\{{\mathbf{B}}_{l_j}^m{\}}_{j=1,l_j\ne k}^s,{\mathbf{\Sigma }}^m\right)$$

$$-\sum _{\begin{array}{c}(l_1\dots l_s):\\ k\in \{l_1,\dots ,l_s\}\end{array}}{\pi }_{(l_1\dots l_s)}^{m+1}{\lambda }_k\Vert {\mathbf{B}}_k{\Vert }_1.$$

$$\displaystyle\mathbf{B}_{k}^{m+1}=\arg\min_{\mathbf{B}_{k}}\sum_{\begin{subarray}{c}(l_{1}\dots l_{s}):\\ k\in\{l_{1},\dots,l_{s}\}\end{subarray}}tr\big{(}(\mathbf{Y}^{*}-\mathbf{X}\mathbf{B}_{k})^{T}\mathbf{z}_{(l_{1}\dots l_{s})}^{m}(\mathbf{Y}^{*}-\mathbf{X}\mathbf{B}_{k})(\mathbf{\Sigma}^{m})^{-1}\big{)}$$

$$+2\sum _{\begin{array}{c}(l_1\dots l_s):\\ k\in \{l_1,\dots ,l_s\}\end{array}}{\pi }_{(l_1\dots l_s)}^{m+1}{\lambda }_k\Vert {\mathbf{B}}_k{\Vert }_1,$$

$${\mathbf{y}}_i^{\ast }={\mathbf{y}}_i-\sum _{l\in \{l_1,\dots ,l_s\}\setminus k}({\mathbf{B}}_l^m{)}^T{\mathbf{x}}_i,{\mathbf{Y}}^{\ast }=({\mathbf{y}}_1^{\ast },\dots ,{\mathbf{y}}_q^{\ast }),$$

$${\mathbf{z}}_{(l_1\dots l_s)}^m=diag\left(z_{1,(l_1\dots l_s)}^m,\dots ,z_{n,(l_1\dots l_s)}^m\right).$$

$${\mathbf{\Sigma }}^{m+1}=$$

$$=$$

$$+$$

$${\mathbf{\Sigma }}^{m+1}=\frac{{\sum }_{(l_1\dots l_s)\in T}(\mathbf{Y}-\mathbf{X}{\sum }_{k=1}^s{\mathbf{B}}_{l_k}^{m+1}{)}^T{\mathbf{z}}_{(l_1\dots l_s)}^m(\mathbf{Y}-\mathbf{X}{\sum }_{k=1}^s{\mathbf{B}}_{l_k}^{m+1})}{{\sum }_{(l_1\dots l_s)\in T}{\sum }_{i=1}^n{z}_{i,(l_1\dots l_s)}^m},$$

Full text

**Precision Therapeutic Biomarker Identification with Application to the Cancer Genome Project

**Hongmei Liu and J. Sunil Rao *Division of Biostatistics, University of Miami, Miami, FL 33136

h.liu7@med.miami.edu

JRao@biostat.med.miami.edu*

Cancer cell lines have frequently been used to link drug sensitivity and resistance with genomic profiles. To capture genomic complexity in cancer, the Cancer Genome Project (CGP) (Garnett et al., 2012) screened 639 human tumor cell lines with 130 drugs ranging from known chemotherapeutic agents to experimental compounds. Questions of interest include: i) can cancer-specific therapeutic biomarkers be detected, ii) can drug resistance patterns be identified along with predictive strategies to circumvent resistance using alternate drugs, iii) can biomarkers of drug synergies be predicted ? To tackle these questions, following statistical challenges still exist: i)biomarkers cluster among the cell lines; ii) clusters can overlap (e.g. a cell line may belong to multiple clusters); iii) drugs should be modeled jointly. We introduce a multivariate regression model with a latent overlapping cluster indicator variable to address above issues. A generalized finite mixture of multivariate regression (FMMR) model in connection with the new model and a new EM algorithm for fitting are proposed. Re-analysis of the dataset sheds new light on the therapeutic inter-relationships between cancers as well existing and novel drug behaviors for the treatment and management of cancer. Key Words: Cancer biomarkers; EM algorithm; finite mixture of multivariate regression model; LASSO; overlapping clustering.

1 Introduction

1.1 Data Description

The use of drugs to selectively target specific genetic alterations in defined patient subpopulations has seen significant successes. One example can be found in the treatment of chronic myeloid leukaemia (CML) where the first consistent chromosomal abnormality associated with a human cancer was identified back in the 1960s. Fast forward to the 1980s where the consequence of this abnormality was discovered to be the production of an abnormal gene called BCR-ABL. Intense drug discovery programs were initiated to shut down the activity of BCR-ABL, and in 1992, imatinib (Gleevec) was developed. In 1998, the drug was tested in CML patients who had exhausted standard treatment options and whose life expectancy was limited, with remarkable results in their blood counts returning to normal. In 2001, the FDA approved imatinib. Today, a once commonly fatal cancer now has a five-year survival rate of 95% (Druker et al., 2006).

Achievements like this largely inspire today’s high throughput screening studies of linking cancer drugs (known or in development) to specific genomic changes which could be used as therapeutic biomarkers. The hope is that such analyses will shed light on biological mechanisms underlying drug sensitivity, tumor resistance and potential drug combination synergies.

Cancer cell lines have frequently been used as a convenient way of conducting such studies. For a systematic search of therapeutic biomarkers to a variety of cancer drugs, the Cancer Genome Project (CGP) (Garnett et al., 2012) screened 639 human tumor cell lines, which represent much of the tissue-type and genetic diversity of human cancers, with 130 drugs. These drugs, including approved drugs, drugs in development as well as experimental tool compounds, cover a wide range of targets and processes involved in cancer biology. A range of 275–507 cell lines were screened for each drug. The effect of a 72h drug treatment on cell viability was examined to derive such measures of drug sensitivity as the half-maximal inhibitory concentration ($IC_{50}$). The cell lines underwent sequencing of 64 known cancer genes, genome-wide analysis of copy number gains and losses, and expression profiling of 14,500 genes.

Given the degree of complexity of this dataset, the multivariate analysis of variance (MANOVA) and the Elastic-Net regression applied in Garnett et al. (2012) are insufficient for precise knowledge discovery. First, the marginal drug-feature associations discovered in MANOVA rarely reflect true relationships, as it is more likely that sensitivity of cancer cells to drugs depends on a multiplicity of genomic and epigenomic features with potential interactions. Second, the Elastic-Net regression fails to concern following issues: 1) since the 639 cell lines come from a variety of cancer tissue types, there is likely additional heterogeneity manifested as subpopulations with overlaps in data; 2) note that the 130 drugs (response variables) are hardly independent, one can improve the prediction accuracy by modeling with multiple drugs (Breiman and Friedman, 1997).

Moreover, there are some direct questions of interest from a subject matter perspective that we want to address. These include, i) can cancer-specific therapeutic biomarkers be detected, ii) can drug resistance patterns be identified along with predictive strategies to circumvent resistance using alternate drugs, iii) can biomarkers of combination therapies be identified to help predict synergies in drug activities ? To tackle these questions and previously discussed statistical challenges, we propose a multivariate regression model with a latent overlapping cluster indicator variable. Fitting procedures inducing concurrent variable selection are introduced.

The rest of the paper is organized as follows. In Section 1.2, we give a selective overview of existing clustering and overlapping clustering methods for general (without response variables) and regression data. In Section 2, a new statistical model is introduced, a generalized FMMR model in connection with the new model and a new EM algorithm for fitting are provided. We also establish a type of consistency optimality for estimation of the generalized FMMR model and perform some small simulation studies to empirically demonstrate this. In Section 3, we put forward another fitting solution to the new model for comparison with the generalized FMMR model. Section 4 contains a comprehensive re-analyses of the CGP data using the proposed method. Discussion is included in Section 5.

1.2 Relevant Statistical Literature Review

Clustering is a well established technique to group data elements based on a measure of similarity. Traditional clustering techniques generate partitions so that each data point belongs to one and only one cluster. It has long been recognized that such ideal partition seldom exists in real data (Needham, 1965). It is more likely that clusters overlap in some parts. To handle the overlapping issue, Lazzeroni and Owen (2002) put forward the well-known plaid model for two-sided overlapping clustering for gene expression data, see also Turner et al. (2005) for improved plaid model and Zhang (2010) for Bayesian plaid model formulation. Its numerical solution, however, produces unsatisfactory cluster retrievals (see Section 3). Other overlapping clustering methods include the “naive” finite mixture (FM) model with a hard threshold on posterior membership probabilities, the probabilistic model (Banerjee et al., 2005) and the multiplicative mixture model based approach (Fu and Banerjee, 2008).

More often interest centers on investigating the relationship between response variables $\mathbf{Y}\in\mathbb{R}^{q}$ and covariates $\mathbf{X}\in\mathbb{R}^{p}$ by fitting a regression model rather than to explore the $\mathbf{X}$ or $\mathbf{Y}$ on its own. In regression analyses with $q=1$, finite mixture of regression (FMR) models are commonly used to capture unobserved cross-sectional heterogeneity in the data (Jedidi et al., 1996). The FMR model postulates that a sample of observations come from a finite mixture of latent partitioning sub-populations with each sub-population represented by a regression model. It was first introduced to statistical literatures by Quandt (1972). DeSarbo and Cron (1988) proposed an EM algorithm (Dempster et al., 1977) based maximum likelihood estimation for the FMR model. Khalili and Chen (2012) put forward a penalized likelihood approach for variable selection in FMR models. As a natural generalization to multivariate responses case, Jones and McLachlan (1992) introduced the FMMR model. Grun and Leisch (2008) proposed a $R$ package flexmix for fitting FMMR models, however it assumes that the response variables are independent.

2 The Proposed Method

2.1 Statistical modeling

For a sample (of cell lines) of $n$ observations, denote $\mathbf{y}_{i}=(y_{i1},\dots,y_{iq})^{T}\in\mathbb{R}^{q}$ a vector of responses ($IC_{50}$ values), $\mathbf{x}_{i}=(x_{i1},\dots,x_{ip_{n}})^{T}\in\mathbb{R}^{p_{n}}$ a vector of predictors (genomic markers) and $\mathbf{\varepsilon}_{i}=(\varepsilon_{i1},\dots,\varepsilon_{iq})^{T}$ a vector of random errors for the $i$th observation. Assume $\mathbf{\varepsilon}_{i}\overset{iid}{\sim}N_{q}(\mathbf{0},\mathbf{\Sigma})$ for $i=1,\dots,n$. The following proposed model allows overlapping clustering for multivariate regression data,

[TABLE]

where $K$ is the total number of clusters, $\mathbf{B}_{k}$ is an unknown $p_{n}\times q$ coefficient matrix for the $k$th cluster, and $P_{ik}\in\{0,1\}$ is 1 if observation $i$ belongs to the $k$th cluster, otherwise 0. In traditional clustering problem, it assumes that each observation belongs to exactly one cluster, namely $\sum_{k=1}^{K}P_{ik}=1$ for all $i$. Therefore we allow $\sum_{k=1}^{K}P_{ik}\geq 1$ so that each observation can belong to multiple clusters.

We provide some interpretation for the clusters in model (1). A cluster $k$ contains a subset of observations for whom $P_{ik}=1$. We postulate that not all genomic features are relevant in describing the cluster $k$, thus assume a sparse coefficient matrix $\mathbf{B}_{k}$. Since the sparse patterns can vary with $k$, each cluster is represented by a unique set of biomarkers. For observations belong to multiple clusters, their response variables are explained by multiple sets of biomarkers, indicative of involving in several biological processes simultaneously.

2.2 The FMMR model

When $\sum_{k=1}^{K}P_{ik}=1$, model (1) can be characterized by a hierarchical structure which ends up with an FMMR model. Consider a latent cluster membership random variable $z_{i}$ for observation $i$ from (1). Given $\sum_{k=1}^{K}P_{ik}=1$, the range of $z_{i}$ equals to $\{1,\dots,K\}$. Assume that $P(z_{i}=k)=\pi_{k}$ for each $k$, then $\sum_{k=1}^{K}\pi_{k}=1$ and $\pi_{k}\geq 0$. Response $\mathbf{y}_{i}$ from model (1) satisfies

[TABLE]

Denote $\mathbf{\Theta}=(\mathbf{B}_{1},\dots,\mathbf{B}_{K},\mathbf{\Sigma},\mathbf{\pi})$ with $\mathbf{\pi}=(\pi_{1},\dots,\pi_{K-1})^{T}$. We can derive the joint density of $\mathbf{y}_{i}$ and $z_{i}$ as

[TABLE]

Summarizing (3) over $z_{i}=1,\dots,K$ yields the FMMR model

[TABLE]

When using the EM algorithm to estimate model (4), $(\mathbf{x}_{i},\mathbf{y}_{i})$ is regarded as incomplete data for missing $z_{i}$, so one works on the complete joint density in (3). For a sample of $n$ observations from (1), the complete log-likelihood function of $\mathbf{\Theta}$ becomes

[TABLE]

where $z_{ik}=I(z_{i}=k)$.

Since we do not expect all genomic markers to be informative, variable selection, to obtain a parsimonious model is necessary. Khalili and Chen (2012) and Khalili and Lin (2013) (for diverging model size) introduced a penalized likelihood approach for variable selection in FMR models, which was shown consistent in variable selection and highly efficient in computation. We define a penalized complete log-likelihood function as

[TABLE]

where

[TABLE]

The following two penalty functions are currently used to meet different demands in variable selection:

1. the $L_{1}$-penalty in LASSO (Tibshirani, 1996) for simultaneous estimation and variable selection

[TABLE]

2. a linear combination of the $L_{1}$- and $L_{2}$-penalty in Elastic-Net (Zou and Hastie, 2005) for simultaneous estimation and selection of grouped features

[TABLE]

where $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$ are respectively the $L_{1}$- and $L_{2}$-norms and the tuning parameters $\lambda_{k},\lambda_{k1},\lambda_{k2}\geq 0$. The $L_{1}$-penalty is singular at the origin, thus can shrink some coefficients to exact 0 for sufficiently large $\lambda_{k}$ or $\lambda_{k1}$ (Fan and Li, 2001).

Their revised EM algorithm can easily be generalized to maximize $\tilde{l}_{n}(\mathbf{\Theta})$ in (6), hence we bypass here. In next section, generalized FMMR model and EM algorithm are devised to fit overlapping multivariate regression data.

2.3 The Generalized FMMR Model

Finite mixture models (including FMR and FMMR) can only fit partitions, as can be seen from Section 2.2. To enable overlapping clustering, one often uses a ”naive” approach by applying a hard threshold $\alpha$ to finite mixture models. This approach assigns a sample to cluster $k$ if its posterior probability of belonging to cluster $k$ is larger than a pre-specified $\alpha$, hence enabling a sample to belong to multiple clusters. As pointed out by Banerjee et al. (2005), this method is problematic because it is not a natural generative model for overlapping clustering, since one underlying assumption of the finite mixture model is that each observation comes from one and only one mixture component. In this section however, we introduce a generalized FMMR model to retrieve overlapping clusters when $\sum_{k=1}^{K}P_{ik}\geq 1$. This generalized model retains its ability in fitting partitions.

Suppose we have $K$ objective clusters indexed by 1 to $K$. These objective clusters refer to clusters defined in (1) and can overlap with each other, resulting in $2^{K}-1$ types of overlapping patterns. Let $K$ equal to 3 for an example, overall overlapping patterns are then composed of $S=\{1,2,3,12,13,23,123\}$. Figure 1 shows the particular part each overlapping pattern in $S$ represents. By definition, overlapping patterns in $S$ are mutually exclusive. And each observation $i$ from (1) belongs to one and only one overlapping pattern.

In the next, we define $2^{K}-1$ hypothetical clusters. Each hypothetical cluster (simply called “cluster” since after) represents an overlapping pattern defined in advance and is indexed by an element in $T$,

[TABLE]

Cluster $(l_{1}\dots l_{s})$ implies its members belonging to objective clusters $l_{1},\dots,l_{s}$.

We now introduce a latent cluster membership random variable $z_{i}$ for observation $i$ from model (1) and characterize the model via a hierarchical structure. Given $\sum_{k=1}^{K}P_{ik}\geq 1$, the range of $z_{i}$ becomes $T$. Further define

[TABLE]

Then vector $\mathbf{\pi}=(\pi_{(l_{1}\dots l_{s})\in T})^{T}$ has

[TABLE]

Response $\mathbf{y}_{i}$ from (1) satisfies

[TABLE]

Let $\mathbf{\Theta}=(\mathbf{B}_{1},\dots,\mathbf{B}_{K},\mathbf{\Sigma},\mathbf{\pi})$. By (7) and (8), the joint density of $\mathbf{y}_{i}$ and $z_{i}$ equals

[TABLE]

Summarizing (9) over $z_{i}$ leads to the generalized FMMR model

[TABLE]

Note that if $\pi_{(l_{1}\dots l_{s})}=0$ for $s>1$, (10) reduces to the traditional FMMR model in (4).

Then the (conditional) log-likelihood function of $\mathbf{\Theta}$ for a sample of $n$ observations from (1) is

[TABLE]

Maximizing above ordinary likelihood function yields non-zero estimates for all regression coefficients. To induce variable selection and remove noise predictors from the regression model, we propose a penalized log-likelihood function

[TABLE]

where $\rho_{nk}(\mathbf{B}_{k})$ is the LASSO or Elastic-Net penalty function in Section 2.2. Notation $\sum_{\begin{subarray}{c}(l_{1}\dots l_{s}):\\ k\in\{l_{1},\dots,l_{s}\}\end{subarray}}$ indicates that it summarizes over $(l_{1}\dots l_{s})\in T$ for which $k\in\{l_{1},\dots,l_{s}\}$. The penalty imposed on $\mathbf{B}_{k}$ is proportional to $\sum_{\begin{subarray}{c}(l_{1}\dots l_{s}):\\ k\in\{l_{1},\dots,l_{s}\}\end{subarray}}\pi_{(l_{1}\dots l_{s})}$, which by definition is proportional to the number of observations involved in the $k$th objective cluster. This is a strategy, similar to Khalili and Chen (2012), of relating the penalty to sample sizes for enhanced power of the method.

2.4 Numerical Solution to the Generalized FMMR Model

We use the renowned EM algorithm for optimization of the generalized FMMR model. The complete log-likelihood function of $\mathbf{\Theta}$ is

[TABLE]

where $z_{i,(l_{1}\dots l_{s})}=I(z_{i}=(l_{1}\dots l_{s}))$.

The penalized complete log-likelihood function for concurrent variable selection is then

[TABLE]

Due to the overlap setting, each $\mathbf{B}_{k}$ involves in multiple clusters. Therefore coefficient matrix $\mathbf{B}_{k}$ for $k=1,\dots,K$ can not be optimized independently like in usual EM algorithm. Instead we sequentially update $\mathbf{B}_{k}$ given the rest are known. Specifically take $\rho_{nk}(\mathbf{B}_{k})$ as a LASSO penalty function for an example, our revised EM algorithm iteratively maximizes $\tilde{l}_{n}(\mathbf{\Theta})$ in two steps:

E-step: Given $\mathbf{\Theta}^{m}$, estimate $z_{i,(l_{1}\dots l_{s})}^{m}$ with its posterior probability by applying the Bayes’ rule to (9) and (10),

[TABLE]

M-step: Given $\mathbf{Z}^{m}=(z_{i,(l_{1}\dots l_{s})}^{m})$, update the mixing proportions $\mathbf{\pi}^{m+1}$ by

[TABLE]

Note that this is obtained by maximizing the leading term of (12) with respect to $\mathbf{\pi}$. This simplified updating scheme however works well in simulation studies.

Given $\mathbf{Z}^{m}$, $\mathbf{\pi}^{m+1}$ and $\mathbf{\Sigma}^{m}$, sequentially update $\mathbf{B}_{k}^{m+1}\mid\{\mathbf{B}_{s}^{m}\}_{s=1,s\neq k}^{s=K}$ by

[TABLE]

By (8),

[TABLE]

where

[TABLE]

[TABLE]

Above (14) is a multivariate regression problem with a LASSO penalty for estimation, which can be solved by the MRCE algorithm (Rothman et al., 2010). If we ignore the covariance structure of $\mathbf{\Sigma}^{m}$ merely in (14), estimation of $\mathbf{B}_{k}^{m+1}$ reduces to $q$ independent LASSO regression problems. In this case, more complex penalty functions such as the Elastic-Net and fused LASSO penalties can easily be applied. Therefore we also have investigated the performance of this simplified strategy in simulation studies. It works surprisingly comparable to original method utilizing the MRCE algorithm for estimation of $\mathbf{B}_{k}^{m+1}$.

Given $\mathbf{Z}^{m}$ and $\{\mathbf{B}_{k}^{m+1}\}_{k=1}^{K}$, update $\mathbf{\Sigma}^{m+1}$ by

[TABLE]

By taking the derivative of (15) according to $\mathbf{\Sigma}^{-1}$, we get

[TABLE]

where $\sum_{(l_{1}\dots l_{s})\in T}\sum_{i=1}^{n}z_{i,(l_{1}\dots l_{s})}^{m}=n$.

Although we can also penalize the inverse covariance matrix $\mathbf{\Sigma}^{-1}$ of (15) like in Rothman et al. (2010), simulations (not presented in the paper) show that penalizing $\mathbf{\Sigma}^{-1}$ results in a lower degree of clustering accuracy than not penalizing $\mathbf{\Sigma}^{-1}$. This may be because by penalizing $\mathbf{\Sigma}^{-1}$ one introduces bias into its estimation, a biased estimate of $\mathbf{\Sigma}^{-1}$ in return deteriorates the estimate of $z_{i,(l_{1}\dots l_{s})}$ in E-step. Thus we choose to not penalize $\mathbf{\Sigma}^{-1}$.

Commencing with an initial value $\mathbf{\Theta}^{0}$, the algorithm iterates between E- and M-steps until the relative change in log-likelihood, |\big{(}l_{n}^{m+1}(\mathbf{\hat{\Theta}})-l_{n}^{m}(\mathbf{\hat{\Theta}})\big{)}/l_{n}^{m}(\mathbf{\hat{\Theta}})|, is smaller than some threshold value, taken as $10^{-5}$ in simulation studies and $10^{-3}$ in real data analysis. Additionally, a cluster, whose mixing proportion is smaller than some threshold value taken as 0.01 in the paper, will be removed during iterations to avoid over estimations.

2.5 Selection of Tuning Parameters and the $K$

In preceding penalized likelihood approach, one needs to choose the values of component-wise tuning parameters $\lambda_{k}$ for $k=1,\dots,K$, which controls the complexity of an estimated model. The data-driven method cross-validation (CV) (Stone, 1977) is frequently adopted in literatures. Here we use a component-wise 10-fold CV method for tuning parameters selection in (14).

Moreover, selection of the number of components $K$ is essential in finite mixture models (including FMR and FMMR). In applications, the choice can be based on prior knowledge of data analysts. With respect to formatted methodologies, information criteria (IC) remains by far the most popular strategy for selection of $K$. See Claeskens et al. (2008) for general treatments on this topic. Here we choose the $K$ minimizing below $IC_{n}$,

[TABLE]

Above $N_{k}$ is the effective number of parameters in the model,

[TABLE]

where $|A|$ calculates the number of nonzero elements in $A$. In (16), $a_{n}$ is a positive sequence depending on $n$. The well known AIC (Akaike, 1974) and BIC (Schwarz et al., 1978) correspond to $a_{n}=2$ and $a_{n}=\log(n)$ respectively. It was shown in Keribin (2000) that under general regularity conditions, BIC can identify the true order of a finite mixture model asymptotically. Feasibility of his results to FMR or FMMR models is unknown yet. We examined the performance of BIC for selection of $K$ in a generalized FMMR model via simulation studies. It obtained a high degree of accuracy.

2.6 The Asymptotic Properties and Simulation Studies

Khalili and Lin (2013) showed that their approach by maximizing a penalized log-likelihood function for estimation of the FMR model is consistent in both estimation and variable selection under certain regularity conditions. Denote $f(\mathbf{w};\mathbf{\Theta})$ the joint density function of data $\mathbf{w}=(\mathbf{x},\mathbf{y})$ with $\mathbf{\Theta}\in\mathbf{\Omega}$. The conditional density function of $\mathbf{y}$ given $\mathbf{x}$ follows an FMR model in Khalili and Lin (2013) and a generalized FMMR model in (10) in our paper. As defined in (11), the penalized log-likelihhood function is a summation of the $\log(f(\mathbf{w}_{i};\mathbf{\Theta}))$’s minus a penalty function. Because the regularity conditions for asymptotic establishments in Khalili and Lin (2013) are made on $f(\mathbf{w};\mathbf{\Theta})$ and the penalty function directly, their theoretical achievements can be extended to our problem.

Denote $\mathbf{B}_{k}^{j}$ the $j$th column of $\mathbf{B}_{k}$ for $j=1,\dots,q$. Let $\mathbf{B}_{k}^{j}=(\mathbf{B}_{k1}^{j},\mathbf{B}_{k2}^{j})$ to divide the coefficient vector into non-zero and zero subsets. Denote $\mathbf{\Theta}^{0}$ the true value of $\mathbf{\Theta}$, $\mathbf{\Theta}^{0}=(\mathbf{\Theta}_{1}^{0},\mathbf{\Theta}_{2}^{0})$ is the corresponding decomposition such that $\mathbf{\Theta}_{2}^{0}$ contains all zero coefficients, and $\hat{\mathbf{\Theta}}_{n}=(\hat{\mathbf{\Theta}}_{n1},\hat{\mathbf{\Theta}}_{n2})$ is its estimate. Since the dimension of $\hat{\mathbf{\Theta}}_{n1}$ increases with $n$, we investigate the asymptotic distribution of its finite linear transformation, $\mathbf{D}_{n}\hat{\mathbf{\Theta}}_{n1}$, where $\mathbf{D}_{n}$ is an $l\times d_{n1}$ constant matrix with a finite $l$ and $d_{n1}$ is the dimension of $\mathbf{\Theta}_{1}^{0}$. Moreover $\mathbf{D}_{n}\mathbf{D}_{n}^{T}\to\mathbf{D}$ and $\mathbf{D}$ is a positive definite and symmetric matrix. We assume that $K$ is independent of the sample size $n$ and known beforehand. Its selection is discussed in Section 2.5.

Lemma 1**.**

Suppose the penalty function $\rho_{nk}(\mathbf{B}_{k})$ satisfies conditions $\mathcal{P}_{0}-\mathcal{P}_{2}$ and the joint density function $f(\mathbf{w};\mathbf{\Theta})$ satisfies conditions $R_{1}-R_{5}$ in Khalili and Lin (2013).

1. If ${p_{n}^{4}\over n}\to 0$, then there exists a local maximizer $\hat{\mathbf{\Theta}}_{n}$ for the penalized log-likelihood $\tilde{l}_{n}(\mathbf{\Theta})$ with

[TABLE]

where $q_{2n}=\max_{kij}\left\{{|\rho_{nk}^{\prime}(b_{kij}^{0})|\over\sqrt{n}}:b_{kij}^{0}\neq 0\right\}$, $b_{kij}^{0}$ is an entry of $\mathbf{B}_{k}^{0}$.

*2. If $\rho_{nk}(\mathbf{B}_{k})$ also satisfies condition $\mathcal{P}_{3}$ in Khalili and Lin (2013) and ${p_{n}^{5}\over n}\to 0$, for any $\sqrt{n/p_{n}}$-consistent maximum likelihood estimate $\hat{\mathbf{\Theta}}_{n}$, as $n\to\infty$ it has:

i. Variable selection consistency:*

[TABLE]

ii. Asymptotic normality:

[TABLE]

where $\mathcal{I}_{1}(\mathbf{\Theta}_{1}^{0})$ is the Fisher information matrix under the true subset model.

Proof.

Write $\mathbf{\Theta}=(\theta_{1},\theta_{2},\dots,\theta_{t_{n}})^{T}$, where $t_{n}$ is the total number of parameters in the model. Then Lemma 1 is a direct extension from Khalili and Lin (2013). ∎

By Lemma 1, $\hat{\mathbf{\Theta}}_{n}$ under the LASSO or Elastic Net penalty has a convergence rate $\sqrt{p_{n}/n}$ via appropriate choice of the tuning parameters. However consistent estimation does not necessarily guarantee consistent variable selection. By the Lemma, the LASSO or Elastic Net penalty does not lead to consistent variable selection, because on one hand $\lambda_{k}$ must be large enough to achieve sparsity, on the other hand the bias term $q_{2n}$ is proportional to the $\lambda_{k}/\sqrt{n}$. This problem can be solved by using adaptive LASSO or adaptive Elastic Net penalty instead, which leads to concurrent estimation and variable selection consistency.

We conducted a sequence of small simulation studies to numerically demonstrate optimality properties of the new method. They are presented in Web Appendix A.

3 Generalization of the Plaid Model

Lazzeroni and Owen (2002) proposed the plaid model to decompose a gene expression data as the sum of overlapping layers (clusters). Each layer contains a subset of genes and samples, while each gene and sample can participate in multiple layers. Denote $x_{ij}$ a gene expression data entry and $K$ is the total number of layers. The plaid model is

[TABLE]

where $\omega_{0},\omega_{k},\alpha_{ik},\beta_{jk}\in\mathbb{R}$, $P_{ik}$ is 1 if gene $i$ is in the $k$th gene-bolck and otherwise 0, $Q_{jk}$ is 1 if sample $j$ is in the $k$th sample-block and otherwise 0. An iterative algorithm was put forward to consecutively search for the layers.

In fact, model (1) is a generalized plaid model by introducing covariates into the model so that it can cluster overlapping multivariate regression data. Therefore, estimation procedures of the plaid model can also be used to estimate (1). In parallel to the plaid model, the parameters are estimated by minimizing the $Q$,

[TABLE]

where $\rho_{nk}(\mathbf{B}_{k})$ is the penalty function defined in Section 2.2.

Detailed procedures for optimizing $Q$ are shown in Algorithm 1 in Web Appendix B, which is a generalization of the improved iterative algorithm by Turner et al. (2005) to multivariate regression. At each time, the algorithm searches for a layer that explains as much of the rest data as possible; once a layer has been found, it is subtracted from the data and remains unchanged; the algorithm stops when no further layers can be found. However this algorithm introduces bias to parameter estimations due to its discrete updating scheme, where previous recruited layers can not be refined as new layers are recruited into the model. To address this issue, we further revised Algorithm 1 to enable joint optimization of all layers. Detailed procedures are presented in Algorithm 2 in Web Appendix B.

4 Therapeutic Biomarker Identification for the CGP Data

We now turn our attention back to analyses of the CGP high throughput drug sensitivity dataset for cancer described in Section 1.1. An updated version of this data

(http://www.cancerrxgene.org/downloads/) contains 707 human tumor cell lines as samples, 140 drugs as response variables, and 13831 genomic features as covariates (including the tissue type, rearrangements, mutation status of 71 cancer genes, continuous copy number data of 426 genes causally implicated in cancer, as well as genome-wide transcriptional profiles). The response variables, made up with $IC_{50}$ values from pairwise drug-cell-line screening, have some missing values. Therefore the data was first filtered by removing cell lines for which less than 50% of the drugs were tested, resulting in 591 cell lines remaining. In the remaining data, about 37.4% of the cell lines are with missing values. These missing values were then imputed via the a random forest imputation algorithm (Ishwaran et al., 2008).

Direct Q1: Can cancer-specific therapeutic biomarkers be detected?

We applied our new method to identify patterns of cancer-specific therapeutic biomarkers. Note that by “cancer-specific” we do not assume separate patterns for each cancer type but rather clusters that may be driven by one or a small number of cancer types. Throughout the analysis, we fixed $K$ to a value of 3, although as previously shown, a BIC model selection approach could also be used. However, fixing the value of $K$ does not limit the interesting findings that we can still find and saves on computational time.

Due to the scale and complexity of the data, the analysis was conducted in three steps. Although simulation studies show that modeling with multiple response variables yields much higher clustering accuracy than modeling with a single response variable, it is unreasonable to simply fit all 140 drugs (responses) in one generalized FMMR model. Because by doing so, one assumes that the cell line assignments to clusters are the same for all 140 drugs, which can hardly be true. Thus in our analysis, we first divided the 140 drugs into several groups and then fitted each drug group with a generalized FMMR model, in which to capture active grouped genomic features, the Elastic-Net penalty along with a simplified updating scheme of $\mathbf{B}_{k}$ discussed in Section 2.4 was used.

In the first step, each drug $c$ was fitted by a generalized FMMR model, from which we got drug-specific cell line assignments, $\hat{\mathbf{Z}}_{c}$, and cluster-wise coefficient estimates, $\hat{\mathbf{B}}_{k}^{c}$ for $k=1,\dots,3$. In the second step, we used the affinity-propagation clustering (APC) algorithm (Frey and Dueck, 2007) to group the 140 drugs based on results from step 1. The grouping was conducted in a two-level nested manner. In the first level, the APC algorithm was applied to all 140 drugs. The pair-wise similarity matrix required by the algorithm as input data was calculated from the Euclidean distance between $\hat{\mathbf{Z}}_{a}$ and $\hat{\mathbf{Z}}_{b}$ where $a$ and $b$ refer to two unique drugs. In the second level, the APC algorithm was re-applied to each first-level drug group. Coefficient estimates $\hat{\mathbf{B}}_{k}^{c}$ for $k=1,\dots,3$ were used to calculate the pair-wise similarities. Consequently, drugs having similar cell line assignments and cluster-wise coefficient estimates were grouped together. Resulting drug groups are shown in Web Figure 6, where the coloring and thickness of connecting lines indicate degree of closeness. In the third step, each second-level drug group $C$ from step 2 was fitted by a generalized FMMR model, from which we got drug-group-specific cell line assignments, $\hat{\mathbf{Z}}_{C}$, and cluster-wise coefficient estimates $\hat{\mathbf{B}}_{k}^{C}$ for $k=1,\dots,3$.

As an illustration, cell line members and coefficient estimates of clusters 1 and 12 for each second-level drug group are shown in Web Figures 7–8. Cell line members of a cluster are depicted via frequencies of the cancer types cell lines belong to. The frequencies are represented by the size of colored bubbles in top panel of the figures. Estimates of other clusters are not presented in the paper due to limited space. It’s very clear that clusters being driven by different cancer types have very different therapeutic genomic profiles and yet no clusters are homogeneous in a particular cancer type. This speaks to the great heterogeneity seen in cancer overall where in fact, cancer is not thought of as a single disease but rather many diseases characterized by different underlying biological changes.

**Direct Q2: Can drug resistance patterns be identified along with predictive strategies to circumvent resistance using alternative drugs?

**

We discussed in the Introduction to the paper the success story of Gleevec which was used for the treatment of CML based on the known specificity of targeting the BCR-ABL gene. It is now becoming more and more common in the treatment of cancer to first sequence known cancer genes in tumor genomes and then design therapies accordingly (Bailey et al., 2014). It is also true that tumors can develop resistance to first line therapies by accumulating mutations which confer resistance. We now show how our methodology can be used to identify predictive strategies for circumventing drug resistance based on mutation data.

Figures 2–3 are similar to Web Figures 7–8, except that they only show a subset of $\hat{\mathbf{B}}_{1}^{C}$ and $\hat{\mathbf{B}}_{12}^{C}$ with regard to the mutations of 71 cancer genes, where $\hat{\mathbf{B}}_{12}^{C}=\hat{\mathbf{B}}_{1}^{C}+\hat{\mathbf{B}}_{2}^{C}$ are coefficient estimates of cluster 12 for drug group $C$. Bottom panel of figure 2–3 are in fact regression based drug-genetic association map with red indicating drug-sensitivity biomarkers and blue indicating drug-resistance biomarkers, conditioning on a cluster of cell lines (upper panel) identified by the generalized FMMR model. We propose that above estimation results can be used for drug repurposing (Martins et al., 2015) for resistant tumors. Take soft tissue cancers as an example of this drug repurposing. We extracted those clusters which contain a high percentage of soft tissue cancers (large purple bubbles), along with the subset cluster-wise coefficient estimates. The extracted information is presented in Web Figure 5. Its upper panel shows the $IC_{50}$ values of soft tissue cell lines in extracted clusters. The drugs in Web Figure 5 were further filtered by only keeping those that have low $IC_{50}$ values without resistance biomarkers and that have high $IC_{50}$ values with resistance biomarkers, leading to Figure 4. In this figure, we conclude that the former set of drugs can be used as alternatives for the later set of drugs which progress resistance to soft tissue cancers. As a positive control of this strategy, drug Gemcitabine was found to be effective in soft tissue cancers resistant to standard chemotherapy (doxyrubicin) in Merimsky et al. (2000).

**Direct Q3: Can drug synergies be predicted?

**

Our desire to allow overlapping cell line clusters becomes apparent in answering this question. Overlapping clusters can be used to guide drug combinations and identify potential drug synergies. From $\hat{\mathbf{Z}}_{C}$, we can get the cluster-wise $IC_{50}$ values for each drug. Web Figure 9 are boxplots of the cluster-wise $IC_{50}$ values for 15 out of 140 drugs. Due to limited space, the rest are not shown. We then identify the drugs for which the magnitudes of $IC_{50}$s of overlapping cluster, say 12, are between those of clusters 1 and 2. For instance, drugs MS.275 and GW843682X show such interesting pattern. We focus on the cluster with highest $IC_{50}$s on average, say cluster 2, and use the coefficient estimates of overlapping cluster to guide the search of another drug to decrease the $IC_{50}$s of cluster 2 towards overlapping cluster.

Take GW843682X for an example, cluster 2 has the highest $IC_{50}$s among clusters 1, 2, and 12. We predicted the $IC_{50}$s in cluster 2 using $\hat{\mathbf{B}}_{1}^{gw843682x}$, $\hat{\mathbf{B}}_{2}^{gw843682x}$ and $\hat{\mathbf{B}}_{12}^{gw843682x}$ respectively. Predicted results are shown in top panel (middle) of Figure 5. By doing so, we want to see if coefficients $\hat{\mathbf{B}}_{1}^{gw843682x}$ or $\hat{\mathbf{B}}_{12}^{gw843682x}$ applied to cell lines in cluster 2 can generate low $IC_{50}$s similar to those in clusters 1 or 12. Since $\hat{\mathbf{B}}_{1}^{gw843682x}$ corresponds to the coefficient estimates of the target drug GW843682X, we need to find another drug producing similar coefficient estimates to $\hat{\mathbf{B}}_{1}^{gw843682x}$, this leads to $\hat{\mathbf{B}}_{3}^{epothilone.b}$. The predicted $IC_{50}$s of cluster 2 using $\hat{\mathbf{B}}_{3}^{epothilone.b}$, $\hat{\mathbf{B}}_{2}^{gw843682x}$ and $\hat{\mathbf{B}}_{3}^{epothilone.b}+\hat{\mathbf{B}}_{2}^{gw843682x}$ are presented in bottom panel (middle) of Figure 5. It shows that the combination of GW843682X and Epothilone.B produces much lower $IC_{50}$s for cluster 2 than GW843682X alone. Similar studies were conducted for MS.275 with respect to its clusters 1, 2, 12 and GW843682X with respect to its clusters 2, 3, 23. Results are shown in left and right columns of Figure 5 respectively. The combination of GW843682X and NVP.BEZ235 decreases the $IC_{50}$s dramatically, while the combination of MS.275 and Bicalutamide does not show such effect. As a positive control, Wildey et al. (2014) pointed out that a group of drugs (including GW843682X and Epothilone B) may provide a practical starting point to investigate combinatorial drug therapies for synergistic effect in small-cell lung cancer (SCLC). Among the group of drugs, they found synergism between GW843682X and CGP60474 via a preliminary study.

5 Discussion

In this paper, we proposed a new model for identifying therapeutic biomarkers for cancer which can answer specific questions regarding sensitivity, resistance and synergy. We used a penalized likelihood approach for the FMMR model which enforced sparsity in the genomic features. To enable overlapping clustering, the FMMR model was then generalized and a new EM algorithm derived for estimation of model parameters. While the noteworthy plaid model can also be generalized for overlapping clustering multivariate regression data, the generalized FMMR model markedly outperforms this method.

Some improvements can be made for future developments on this work. First, other useful penalty functions for multivariate regression estimations can be adopted, such as the $MAP$ (MAster Predictor) penalty in Peng et al. (2010), the $L_{2}SVS$ method by Similä and Tikka (2007) and the $L_{\infty}SVS$ method by Turlach et al. (2005). Second, the multivariate normal distribution assumption on $\mathbf{Y}_{i}$ can be extended to other more flexible parametric families of multivariate distributions, such as the skew-normal distribution (Azzalini, 2005) and the multivariate skew-$t$ distribution (Chen et al., 2014), for frequent presence of skewness and kurtosis in real data.

As described in the CGP data analyses, we had about 37% of the final analyses observations that had missing data and which we dealt with by using random forest imputation. It would be of interest to explore the impact of this step more thoroughly. An alternative would be to generalize to this problem, the recently developed E-MS algorithm Jiang et al. (2015) for model selection with incomplete data.

Acknowledgements

Hongmei Liu was supported by a doctoral fellowship from the Sylvester Comprehensive Cancer Center at the University of Miami.

Supplementary Materials

Web Appendices, Tables, and Figures referenced in Sections 2.6, 3 and 4 are available with this paper at the Biometrics website on Wiley Online Library.

Appendix A Simulation Studies

A.1 Design of the Simulations

Two scenarios were designed to investigate the non-overlapping clustering and overlapping clustering of the new method respectively. In each scenario, data were generated from model (2.1) with $K=3$, $p_{n}=15$, $q=3$ and $n=150,450$. Predictors $\mathbf{x}_{i}$, $i=1,\dots,n$, were drawn independently from $N_{p_{n}}(\mathbf{0},\mathbf{\Sigma}_{1})$ with $\mathbf{\Sigma}_{1}(i,j)=0.5^{|i-j|}$. Random errors $\mathbf{\varepsilon}_{i}$, $i=1,\dots,n$, were drawn independently from $N_{q}(\mathbf{0},\mathbf{\Sigma}_{2})$ with $\mathbf{\Sigma}_{2}(i,j)=0.75^{|i-j|}$.

The sparse coefficient matrices $\mathbf{B}_{k},k=1,2,3$, were generated as

[TABLE]

where $\otimes$ indicates the element-wise product. Each entry of $\mathbf{W}$ was drawn independently from $N(0,1)$, each entry of $\mathbf{S}$ was drawn independently from the Bernoulli distribution $B(1,p_{1})$, and each row of $\mathbf{T}$ (either all 1 or all 0) was determined by an independent draw from $B(1,p_{2})$. As a result, we expect $p_{1}p_{2}p_{n}$ relevant predictors for each response variable, and $(1-p_{2})p_{n}$ predictors are expected to be irrelevant to all $q$ response variables. We let $p_{1}=0.5$ and $p_{2}=0.9$. Each simulation was repeated for 50 times. A new sequence of $\mathbf{B}_{k},k=1,2,3$ and a new set of data $(\mathbf{X},\mathbf{Y})$ were generated for each repetition. We assumed that $K$ is known, so that performance evaluation of the new method is not interfered by selection of $K$. Its selection is separately examined in Section A.4.

Scenario 1: The 3 clusters are non-overlapping, each cluster contains $n/3$ observations.

Scenario 2: 70% of the observations involve in single cluster, 22% of the observations involve in two clusters and 8% of the observations involve in three clusters.

A.2 Quality Measures of Cluster Recovery

Quality measures in Baeza-Yates et al. [1999] have been used to evaluate the clustering performance of the proposed method. Denote $N_{S}$ the number of elements in $S$. Below quality measures were proposed to compare a retrieved cluster $\hat{A}$ with a target cluster $A$,

[TABLE]

The $F_{1}$ measure, taken as the harmonic mean of specificity and sensitivity, gives an overall measure of the clustering.

Moreover, we employed the one-to-one correspondence match approach [Turner et al., 2005] to evaluate a sequence of retrieved clusters. It includes about three steps. First make the number of retrieved clusters to be the same as the number of target clusters by adding null clusters to retrieved clusters or dropping addition poorly retrieved clusters. Retrieved clusters are then matched to target clusters in pair. Finally, calculate the mean quality measures for the sequence of paired clusters.

A.3 Results of the Simulations

A sequence of methods were implemented for comprehensive comparison with our new method, including the generalized plaid model (plaid) and its revised counterparts (aplaid); the R package flexmix (EM) where response variables are treated as independent and overlap issues are not considered; the generalized FMMR model with response variables treated as independent (gEMseplasso0); the generalized FMMR model with separate LASSO estimation of $\mathbf{B}_{k}$ (gEMseplasso); and finally the original generalized FMMR model with MRCE estimation of $\mathbf{B}_{k}$ (gEMmrce).

Their performances were evaluated via the three quality measures introduced in Section A.2 and the sum of squared errors (SSE) of coefficients estimates. Results are summarized in Web Figure 6–7. In scenario 1 when there is no overlap, EM and gEMseplasso0 do equally well showing that the generalized FMMR model retains the ability to recover non-overlapping clusters. By taking into account of the covariance structure among response variables, gEMseplasso and gEMmrce outperform EM and gEMseplasso0, but there is little difference between gEMseplasso and gEMmrce. In scenario 2 when there are 30% overlaps, specificity of EM drops down dramatically, resulting in a poor $F_{1}$ measure. The gEMseplasso0 does much better than EM by estimating the overlaps, while it is defeated by gEMseplasso and gEMmrce. In this scenario, gEMmrce outperforms gEMseplasso when $n=150$. However, they converge very fast with $n$, both achieving a median 95% clustering accuracy when $n=450$. The plaid and aplaid have a poor performance in all scenarios.

A.4 Additional Simulations

We conducted another simulation for scenario 2 to evaluate the accuracy of BIC in selection of $K$. In the simulation, $K$ is unknown in candidate procedures EM and gEMseplasso, and is chosen by minimizing the BIC. Results are summarized in Web Figure 8 and Web Table 1. Given 30% overlaps in simulated data, BIC has a poor performance in identifying the true $K$ with the EM procedure. Whereas by fitting the overlapping clusters with gEMseplasso, BIC obtains a much higher accuracy in selection of $K$ (see Web Table 1). Moreover, although $K$ is misspecified for 12% of the times when $n=450$, gEMseplasso still has achieved a median 93% clustering accuracy (see Web Figure 8).

An additional simulation was conducted for scenario 2 to investigate the advantage of modeling with multivariate response variables in the new method. We compared the performance of modeling with all 3 response variables versus modeling with each of the 3 response variables in the generalized FMMR model. Results are summarized in Web Figure 9. Using multivariate response variables leads to a much higher clustering accuracy.

Appendix B Two algorithms of the Generalized Plaid Model

Web Tables

Web Figures