Joint and Individual Component Regression

Many fields of scientific research involve the analysis of heterogeneous data. In particular, data may appear in the form of multiple matrices, with data heterogeneity arising from either variables or samples. One example is the multi-view/source data, which measure different sets of variables on the same set of samples. The sets of variables may come from different platforms/sources/modalities. For instance, in genomics studies, measurements are collected as different biomarkers, such as mRNA and miRNA (Muniategui et al.,, 2013). Another example is the multi-group data, which measure the same set of variables on disparate sets of samples, which leads to heterogenous subpopulations/subgroups in the entire population. For instance, in the Alzheimer’s Disease (AD) study, subjects can have different subtypes, such as Normal Control (NC), Mild Cognitive Impairment (MCI), and AD.

We study the classical regression problem with one continuous response for multi-group data. Although there are many well-established regression techniques for homogeneous data (Hoerl and Kennard,, 1970; Tibshirani,, 1996), they may not be suitable for multi-group data. One naive approach is to ignore data heterogeneity and fit a global model using these techniques. However, a single global model can be too restrictive because the diverse information from different subgroups may not be identified. On the other hand, one can train separate group-specific models. Despite its flexibility, the information that is jointly shared across different groups is not sufficiently captured. Therefore, it is desirable to build a flexible statistical model that can simultaneously quantify the jointly shared global information and individual group-specific information for heterogeneous data.

There are several existing methods proposed in the literature under the context of regression for multi-group data. Meinshausen and Bühlmann, (2015) took a conservative approach and proposed a maxmin effect method that is reliable for all possible subsets of the data. Zhao et al., (2016) proposed a partially linear regression framework for massive heterogeneous data, and the goal is to extract common features across all subpopulations while exploring heterogeneity of each subpopulation. Tang and Song, (2016); Ma and Huang, (2017) proposed fused penalties to estimate regression coefficients that capture subgroup structures in a linear regression framework. Wang et al., (2018) proposed a locally-weighted penalized model to perform subject-wise variable selection. Wang et al., (2023) proposed a factor regression model for heterogeneous subpopulations under the high-dimensional factor decomposition. However, these models either are not specifically designed to identify the globally-shared and group-specific structures, or impose strong theoretical assumptions on the covariates. On the other hand, there exist some works that study this manner for multi-source data. Lock et al., (2013) proposed JIVE to learn joint and individual structures from multiple data matrices by low-rank approximations. Some extensions of JIVE can be found in Feng et al., (2018); Gaynanova and Li, (2019). All of these decomposition methods are fully unsupervised. Recently, Li and Li, (2021) proposed a supervised integrative factor regression model for mult-source data and studied its statistical properties with hypothesis tests. Palzer et al., (2022) proposed sJIVE that extends JIVE with the supervision from the response. These methods are supervised, but both focused on regressions for multi-source data.

In this paper, we consider the supervised learning problem of predicting a response with multi-group data. We propose a Joint and Individual COmponent Regression (JICO), a novel latent component regression model that covers JIVE as a special case. Our proposed model decomposes the response into jointly shared and group-specific components, which are driven by low-rank approximations of joint and individual structures from the predictors respectively. The joint structure shares the same coefficients across all groups, whereas individual structures have group-specific coefficients. Moreover, by choosing different ranks of joint and individual structures, our model covers global and group-specific models as special cases. To estimate JICO, we propose an iterative algorithm to solve for joint and individual scores using latent component representation. To construct the latent scores, we utilize the Continuum Regression (CR) (Stone and Brooks,, 1990), which provides a unified framework that covers OLS, PLS, and PCR as special cases. Some follow-up studies and modern extensions of CR can be found in Björkström and Sundberg, (1996); Lee and Liu, (2013). Embracing this flexibility and generaliziblity from CR, our proposed JICO model extends to the heterogeneous data setup and is able to achieve different model configurations on the spectrum of CR under this more complicated setting. JICO attains a good balance between global and group-specific models, and further achieves its flexibility by extending CR.

The rest of this paper is organized as follows. In Section 2, we briefly review JIVE and introduce our proposed JICO model. We further provide sufficient conditions to make JICO identifiable. In Section 3, after two motivating special cases, we introduce our iterative algorithm. In Sections 4 and 5, we evaluate the performance of JICO by simulation studies and real data analysis on the Alzheimer’s disease dataset, respectively. In Section 6, we conclude this paper with some discussion and possible extensions.

Suppose we observe data pairs $(\mathbf{X}_{g},\bm{Y}_{g})_{g=1}^{G}$ from $G$ groups, where $\mathbf{X}_{g}\in\mathbb{R}^{n_{g}\times p}$ and $\bm{Y}_{g}\in\mathbb{R}^{n_{g}}$ are the data matrix and the response vector for the $g$th group respectively. In this setting, each data matrix has the same set of $p$ explanatory variables, whereas the samples vary across groups. We let $\mathbf{X}=[\mathbf{X}_{1}^{\prime},\ldots,\mathbf{X}_{G}^{\prime}]^{\prime}\in\mathbb{R}^{n\times p}$ and $\bm{Y}=[\bm{Y}_{1}^{\prime},\ldots,\bm{Y}_{G}^{\prime}]^{\prime}\in\mathbb{R}^{n}$, where $n=\sum_{g=1}^{G}n_{g}$.

Our model is closely related to JIVE, which provides a general formulation to decompose multiple data matrices into joint and individual structures. The JIVE decomposes $\mathbf{X}_{g}$ as

where $\mathbf{J}_{g}\in\mathbb{R}^{n_{g}\times p}$ represents the submatrix of the joint structure of $\mathbf{X}_{g}$, $\mathbf{A}_{g}\in\mathbb{R}^{n_{g}\times p}$ represents the individual structure of $\mathbf{X}_{g}$, and $\mathbf{E}_{g}\in\mathbb{R}^{n_{g}\times p}$ is the error matrix. We consider that $\bm{Y}_{g}$ has a similar decomposition into joint and individual signals

where $\mathbf{e}_{g}\in\mathbb{R}^{n_{g}}$ is the noise from the $g$-th group. Let $\tilde{\mathbf{X}}_{g}=\mathbf{J}_{g}+\mathbf{A}_{g}$ and $\tilde{\bm{Y}}_{g}=\bm{J}^{Y}_{g}+\bm{A}^{Y}_{g}$ be the noiseless counterparts of $\mathbf{X}_{g}$ and $\bm{Y}_{g}$. Lemma 1 gives conditions to ensure that $\mathbf{J}_{g}$, $\mathbf{A}_{g}$, $\bm{J}_{g}^{Y}$, and $\bm{A}_{g}^{Y}$ are identifiable.

Lemma 1 shows that $\tilde{\mathbf{X}}_{g}$ can be uniquely decomposed into the sum of $\mathbf{J}_{g}$ and $\mathbf{A}_{g}$ if we require them to satisfy conditions $(i)-(iii)$, following similar statements as in Feng et al., (2018). To ensure the unique decomposition of $\tilde{\mathbf{Y}}_{g}$, we need to further require $\text{col}(\mathbf{J}_{g})\perp\text{col}(\mathbf{A}_{g})$, which is different from Palzer et al., (2022), that requires $\text{row}(\mathbf{J}_{g})\perp\text{row}(\mathbf{A}_{g})$.

In practice, only $\mathbf{X}_{g}$ and $\bm{Y}_{g}$ are observable. In Lemma 2, we show in $(a)$ that the identifiable conditions in Lemma 1 can still be achieved given the observed $\{\mathbf{X}_{g},\bm{Y}_{g}\}_{g=1}^{G}$. Moreover, in Lemma 2$(b)$, we show that if $\mathbf{J}_{g}$ and $\mathbf{A}_{g}$ are assumed to have low ranks, they can be further decomposed as $\mathbf{J}_{g}=\mathbf{S}_{g}\mathbf{U}$, where $\mathbf{S}_{g}$ is a $n_{g}\times K$ score matrix, $\mathbf{U}$ is a $K\times p$ loading matrix, and $K=\text{rank}(\mathbf{J}_{g})$; and $\mathbf{A}_{g}=\mathbf{T}_{g}\mathbf{U}_{g}$, where $\mathbf{T}_{g}$ is a $n_{g}\times K_{g}$ score matrix and $\mathbf{U}_{g}$ is a $K_{g}\times p$ loading matrix, and $K_{g}=\text{rank}(\mathbf{A}_{g})$. Under this formulation, if $\mathbf{S}_{g}^{\prime}\mathbf{T}_{g}=\bm{0}_{K\times K_{g}}$, then (2.1) and (2.2) can be expressed as

where $\bm{\alpha}\in\mathbb{R}^{K}$ and $\bm{\alpha}_{g}\in\mathbb{R}^{K_{g}}$ are the coefficients of the joint and individual components respectively. Model (2.4) gives a unified framework to model multi-group data. When $K=0$, the joint term $\mathbf{S}_{g}\bm{\alpha}$ vanishes and (2.4) reduces to a group-specific model of $\bm{Y}_{g}=\mathbf{T}_{g}\bm{\alpha}_{g}+\mathbf{e}_{g}$. On the other hand, when $K_{1}=\cdots=K_{g}=0$, the individual term $\mathbf{T}_{g}\bm{\alpha}_{g}$ vanishes and (2.4) reduces to a global model of $\bm{Y}_{g}=\mathbf{S}_{g}\bm{\alpha}+\mathbf{e}_{g}$. When $K\neq 0$ and $K_{g}\neq 0$, (2.4) has both global and group-specific components, thus lies between the above two extreme cases.

Corollary 1 follows directly from the proof of Lemma 2$(b)$. As a remark, the columns of $\mathbf{W}$ and $\mathbf{W}_{g}$ form the sets of bases that span the row spaces of $\mathbf{J}_{g}$ and $\mathbf{A}_{g}$ respectively. Hence, $\mathbf{W}^{\prime}\mathbf{W}_{g}=\bm{0}$ is a sufficient and necessary condition for $\text{row}(\mathbf{J}_{g})\perp\text{row}(\mathbf{A}_{g})$. Moreover, note that in Lemma 2$(b)$, $\mathbf{S}_{g}^{\prime}\mathbf{T}_{g}=\bm{0}$ directly implies that $\mathbf{A}_{g}^{\prime}\mathbf{J}_{g}=\bm{0}$, the latter being a sufficient condition for $\text{col}(\mathbf{J}_{g})\perp\text{col}(\mathbf{A}_{g})$. Therefore, in Corollary 1, $\mathbf{W}^{\prime}{\mathbf{X}}_{g}^{\prime}{\mathbf{X}}_{g}\mathbf{W}_{g}=\bm{0}$ provides a sufficient condition for $\text{col}(\mathbf{J}_{g})\perp\text{col}(\mathbf{A}_{g})$, which satisfies one of the identifiability constraints for the unique decomposition of $\mathbf{X}_{g}$ in Lemma 1. In Section 3, we describe the algorithm to solve for $\mathbf{W}$ and $\mathbf{W}_{g}$ respectively.

The key to estimate (2.3) and (2.4) is the constructions of score matrices $\mathbf{S}_{g}$ and $\mathbf{T}_{g}$. To motivate our estimation procedure, in Sections 3.1, we discuss the joint and individual score estimation under two special cases respectively. In Section 3.2, we introduce an iterative algorithm for the general case.

One significant advantage of our proposed model is its flexibility of lying in between global and group-specific models. Moreover, the choice of parameter $\gamma$ in CR allows it to identify the model that best fits the data. In this section, we conduct multiple simulation studies to further demonstrate the advantage of our proposed model.

We consider three simulation settings in this section. In the first two settings, we generate data according to models that contain both global and group-specific components. The data are generated in a way that PCR and PLS solutions are favored respectively. In the last setting, we simulate data from two special cases: a global model and a group-specific model. The data are simulated so that the OLS is favored for both cases. For all three settings, JICO can adaptively choose the correct model parameter $\gamma$ so that it has the optimal performance. Moreover, we further illustrate how the rank selection impacts the performance of JICO by examining the results using mis-specified ranks.

We fix $G=2,p=200,n_{1}=n_{2}=50$. In each replication, we generate 100 training samples to train the models and evaluate the corresponding Mean Squared Error (MSE) in an independent test set of 100 samples. We repeat simulations for 50 times.

For $g=1,\ldots,G$, we generate $\mathbf{X}_{g}$ as i.i.d. samples from $\mathcal{N}(\bm{0},\mathbf{I}_{p\times p})$. For the sake of simplicity, we generate $\bm{Y}_{g}$ by the following model with two latent components:

where $\bm{S}_{g}=\mathbf{X}_{g}\mathbf{w}\in\mathbb{R}^{n_{g}}$ is the joint latent score vector with an coefficient $\alpha$, $\bm{T}_{g}=\mathbf{X}_{g}\mathbf{w}_{g}\in\mathbb{R}^{n_{g}}$ is the individual latent score vector with an coefficient $\alpha_{g}$, and $\mathbf{e}_{g}$ is generated as i.i.d. samples from $\mathcal{N}(0,0.04)$. Here, $\mathbf{w}$ and $\mathbf{w}_{g}$ are all $p\times 1$ vectors, and they are constructed such that $\mathbf{w}^{\prime}\mathbf{w}_{g}=0$. We vary the choices of $\mathbf{w}$, $\mathbf{w}_{g}$, $\alpha$, and $\alpha_{g}$, which will be discussed later.

We apply our proposed method to analyze data from the Alzheimer’s Disease (AD) Neuroimaging Initiative (ADNI). It is well known that AD accounts for most forms of dementia characterized by progressive cognitive and memory deficits. The increasing incidence of AD makes it a very important health issue and has attracted a lot of scientific attentions. To predict the AD progression, it is very important and useful to develop efficient methods for the prediction of disease status and clinical scores (e.g., the Mini Mental State Examination (MMSE) score and the AD Assessment Scale-Cognitive Subscale (ADAS-Cog) score). In this analysis, we are interested in predicting the ADAS-Cog score by features extracted from 93 brain regions scanned from structural magnetic resonance imaging (MRI). All subjects in our analysis are from the ADNI2 phase of the study. There are 494 subjects in total in our analysis and 3 subgroups: NC (178), eMCI (178) and AD (145), where the numbers in parentheses indicate the sample sizes for each subgroup. As a reminder, NC stands for the Normal Control, and eMCI stands for the early stage of Mild Cognitive Impairment in AD progression.

For each group, we randomly partition the data into two parts: 80% for training the model and the rest for testing the performance. We repeat the random split for 50 times. The testing MSEs and the corresponding standard errors are reported in Table 2. Both groupwise and overall performance are summarized. We compare our proposed JICO model with four methods: ridge regression (Ridge), PLS, and PCR. We perform both a global and a group-specific fit for Ridge, PLS, and PCR, where the regularization parameter in Ridge and the number of components in PCR or PLS are tuned by 10-fold cross validation (CV). For our proposed JICO model, we demonstrate the result by fitting the model with fixed $\gamma=0,0.25,1,\infty$, or tuned $\gamma$ respectively. In practice, using exhaustive search to select the optimal values for $K$ and $K_{g}$ can be computationally cumbersome, because the number of combinations grows exponentially with the number of candidates for each parameter. Based on our numerical experience, we find that choosing $K_{g}$ to be the same does not affect the performance on prediction too much. Details are discussed in Appendix D. Therefore, in all these cases, the optimal ranks for JICO are chosen by an exhaustive search in $K\in\{0,1\}$ and $K_{1}=K_{2}=K_{3}\in\{0,1\}$ to see which combination gives the best MSE. We choose $K$ and $K_{g}$ to be small to improve our model interpretations. The optimal value of $\gamma$ is chosen by 10-fold CV.

As shown in Table 2, JICO performs the best among all competitors. Fitting JICO with $\gamma=0.25$ yields the smallest overall MSE. JICO with parameters chosen by CV performs slightly worse, but is still better than the other global or group-specific methods. The results of JICO with $\gamma=0,1$ and $\infty$, which correspond to OLS, PLS, and PCR, are also provided in Table 2. Even though their prediction is not the best, an interesting observation is that they always have better performance than their global or group-specific counterparts. For example, when $\gamma=1$, JICO has much better overall prediction than the group-specific PLS. This indicates that it is beneficial to capture global and individual structures for regression when subpopulations exist in the data.

In Table 2, global models perform the worst, because they do not take group heterogeneity into consideration. The group-specific Ridge appears to be the most competitive one among group-specific methods. Note that for the AD group, our JICO model with $\gamma=0.25$ or tuned $\gamma$ outperforms the group-specific Ridge method by a great margin.

To get our results more interpretable, we further apply the JICO model to NC and AD groups. We run 50 replications of 10-fold CV to see which combination of tuning parameters gives the smallest overall MSE. The best choice is $\gamma=\infty,K=1,K_{NC}=K_{AD}=3$. Then, we apply JICO using this choice and tuning parameters and give the heatmaps of the estimated $\hat{\mathbf{J}}_{g}$ (left column) and $\hat{\mathbf{A}}_{g}$ (right column) in Figure 5.1. Rows of heatmaps represent samples and columns represent MRI features. We use the Ward’s linkage to perform hierarchical clustering on the rows of $\hat{\mathbf{J}}_{g}$, and arrange the rows of $\hat{\mathbf{J}}_{g}$ and $\hat{\mathbf{A}}_{g}$ in the same order for each group. Moreover, we apply the same clustering algorithm to the columns of $\hat{\mathbf{J}}_{g}$ to arrange the columns in the same order across the two disease groups for both joint and individual structures. Figure 5.1 shows that JICO separates joint and individual structures effectively. The joint structures across different disease groups share a very similar pattern, whereas the individual structures appear to be very distinct. We further magnify the right column of Figure 5.1 in Figure 5.2 with the brain region names listed. We find that the variation in $\hat{\mathbf{A}}_{g}$ for the AD group is much larger than the counterpart for the NC group. We highlight the brain regions that differ the most between the two groups. The highlighted regions play crucial roles in human’s cognition, thus important in AD early diagnosis (Michon et al.,, 1994; Killiany et al.,, 1993). For example, Michon et al., (1994) suggested that anosognosia in AD results in part from frontal dysfunction. Killiany et al., (1993) showed that the temporal horn of the lateral ventricles can be used as antemortem markers of AD.

In this paper, we propose JICO, a latent component regression model for multi-group heterogeneous data. Our proposed model decomposes the response into jointly shared and group-specific components, which are driven by low-rank approximations of joint and individual structures from the predictors respectively. For model estimation, we propose an iterative procedure to solve for model components, and utilize CR algorithm that covers OLS, PLS, and PCR as special cases. As a result, the proposed procedure is able to extend many regression algorithms covered by CR to the setting of heterogeneous data. Extensive simulation studies and a real data analysis on ADNI data further demonstrate the competitive performance of JICO.

JICO is designed to be very flexible for multi-group data. It is able to choose the optimal parameter to determine the regression algorithm that suits the data the best, so that the prediction power is guaranteed. At the same time, it is also able to select the optimal joint and individual ranks that best describe the degree of heterogeneity residing in each subgroup. The JICO application to ADNI data has effectively illustrated that our proposed model can provide nice visualization on identifying joint and individual components from the entire dataset without losing much of the prediction power.

The data that support the findings of this paper are available in the Alzheimer’s Disease Neuroimaging Initiative at https://adni.loni.usc.edu/.