Selecting PCA components which separate groups

Question

I frequently used to diagnose my multivariate data using PCA (omics data with hundreds of thousands of variables and dozens or hundreds of samples). The data often come from experiments with several categorical independent variables defining some groups, and I often have to go through a few components before I can find the ones that show a separation between the groups of interest. I have came up with a rather primitive way of finding such discriminating components, and I wonder

1. to what extent this is reasonable / justifiable, and
2. whether there are better ways of achieving the same.

Note that this is exploratory. Before I convince anyone else, I want to convince myself. If I see that there are components that clearly distinguish the groups of interest (e.g. control vs treatment), even if they are responsible for a minor portion of the variance of the responses, I trust it more than a result from, say, supervised machine learning.

Here is my approach. I will use the "metabo" example data set from pca3d in R.

The idea is to assess how much variance of each of the component can be explained by the independent variable. For this, I calculate a simple model for each component and use $R^2$ as a metric to order the components from "most interesting" to "least interesting".

require( pca3d )
# data on metabolic profiles of TB patients and controls
data( metabo )
# first column is the independent variable
pca <- prcomp( metabo[,-1], scale.= T ) 

# create a model for each component
lm.m <- lm( pca$x ~ metabo[,1] )
lm.s <- summary( lm.m )
lm.r2 <- sapply( lm.s, function( x ) x$r.squared )
plot( lm.r2, type= "l" )
text( 1:length( lm.r2 ), lm.r2, 1:length( lm.r2 ), pos= 3 )

Here is the result. The plot shows the percentage in variance of each component explained by the independent variable in metabo[,1].

We can sort the components by $r^2$ to find out which ones to display with order( lm.r2, decreasing= TRUE ); the first three components are 2, 1 and 7.

pca3d( pca, components= c( 1, 2, 7 ), group= metabo[,1] )

Here is the plot:

(The red and green categories are two groups of subjects who are not patients, and it is to be expected that they cannot be distinguished.)

To reformulate my questions,

1. Does this approach make sense to you? My problem is that it looks too much like data dredging. Also, intuitively I think maybe I should turn the table and ask what part of the variance in the independent variable is explained by each variable? Finally, I'm (almost) sure that I am reinventing the wheel, poorly, so my second question is
2. Is there anything better?

Note that I do not want to switch to partial least squares or anything similar at this stage; I just want to diagnose the PCA in the context of my classification.

Answer 1

The answer to your question #1 is yes, your solution amounts to data dredging. The answer to your question #2 is yes, there are superior methods in the literature.

The central problem with your approach is that you are not addressing the high-dimensional data problem, i.e. problems that arise when $n << p$. Your solution is quite arbitrary and lacks any sort of theoretical justification: I will point you to some literature that can help you find adequate methods below.

You are running an analysis that resembles principal components regression, except that you have swapped your independent and dependent variables, resulting in a large multivariate (as opposed to multiple) regression analysis. Multivariate regression requires that your sample size be larger than the number of dependent variables, a requirement which you are thoroughly violating in your example.

If you are truly committed to running PCA on your data and then using multivariate regression, you must use an appropriate method. For example, look into MRCE and related methods [1].

However, despite a few puzzling comments you have made, everything in your analysis as currently presented suggests that your ultimate goal is to identify relationships between a large set of continuous variables (metabo[,-1]) and a single categorical variable (metabo[,1]). PCA is a poor way of accomplishing this. There are two general classes of solutions to this problem in the high-dimensional case: first, solutions that assume sparsity, and solutions that assume a factor structure.

Sparsity-based solutions typically assume that only a very small proportion of variables are in fact related to the categorical variable of interest, and attempt to find this small subset; for example see DALASS [2]. Factor-structure-based methods assume that your discriminator variables are manifestations of underlying latent variables with a true relationship to the categorical variable. An example of this class of methods is DLDA [3].

Note that I am not necessarily recommending any methods I have mentioned for your data; you must carefully consider your goals and a priori knowledge of the problem in selecting an appropriate method.

[1] Rothman, Levina, Zhu (2010). Sparse Multivariate Regression With Covariance Estimation. Journal of Computational and Graphical Statistics, Volume 19, Number 4, Pages 947–962.

[2] Nickolay T. Trendafilov, Ian T. Jolliffe, DALASS: Variable selection in discriminant analysis via the LASSO, Computational Statistics & Data Analysis, Volume 51, Issue 8, 1 May 2007, Pages 3718-3736.

[3] Yu, Yang (2001). A direct LDA algorithm for high-dimensional data with application to face recognition. Pattern Recognition 34, 2067-2070.

Answer 2

@ahfoss already pointed you to LDA as the classification analogon to PCA. Actually, these two methods are related to each other and also to PLS:

nature of dependent variable (for supervised)     unsupervised    supervised
or structure of data (unsupervised)
continuous                                        PCA             PLS
factor/groups/classes                                             LDA

The relation is in the projection: PCA projects data so that the variance-covariance matrix of the scores will be $\mathbf I$. LDA does a similar projection, but instead of the variance-covariance matrix of the whole data set, the pooled within-class variance-covariance matrix becomes $\mathbf I$.

As @ahfoss also said, the usual LDA is not feasible for $n \ll p$. But PLS-LDA is a feasible way of dealing with that, see e.g. Barker, M. & Rayens, W.: Partial least squares for discrimination, J Chemom, 17, 166-173 (2003). DOI: 10.1002/cem.785 (this paper also discusses the relation between PLS and LDA).

PLS can be seen as a regularization like the LASSO, and also sparse PLS is available (though I haven't used it: my data is more suitable to normal PLS, which doesn't assume sparsity). For a nice discussion of different regularization methods, see e.g. the Elements of Statistical Learning.

One nice property of doing PLS-LDA is that you can write the final model in bilinear form if you take care to use the same centering for both steps. This allows to interpret the model just the same way you'd interpret the LDA model that you could calculate if you had more cases and $n \gg p$. (In case you're forced to do PCA, this applies for PCA-LDA as well).

Say, PLS scores $\mathbf T = $ data matrix $\mathbf X \times$ PLS weights $\mathbf W$ and
for a usual LDA, the scores $\mathbf L = $ data matrix $\mathbf X\times$ LDA coefficients $\mathbf B$.

Then doing the LDA in PLS (X-)score space, we get:
$\mathbf L'^{(n \times k - 1)} = \mathbf T^{(n \times m)} \mathbf B'^{(m \times k - 1)}$
$\phantom{\mathbf L^{(n \times k - 1)}} = \mathbf X^{(n \times p)} \mathbf W^{(p \times m)} \mathbf B'^{(m \times k - 1)}$
$\phantom{\mathbf L^{(n \times k - 1)}} = \mathbf X^{(n \times p)} \mathbf B''^{(p \times k - 1)}$
The dashes mark that these LDA scores ($\mathbf L'$) may be (slightly) different from LDA scores you'd obtain without regularization, and so are the coefficients $\mathbf B''$. $\mathbf B'$ will be very different from $\mathbf B$ as they apply to PLS score space and not to the original data space.

Practical note: in case you work in R, I have a package under development that provides PLS-LDA and PCA-LDA models. Let me know if you'd like to give it a try.

---

In order to avoid data dredging, you need to validate your final model (= measure its performance) with independent data.

Independent here means that this case (patient?) did not contribute to the model fitting in any way. In particular,

• did not enter any kind of preprocessing that involves multiple cases, such as centering or standardization
• did not enter the PCA/PLS/... calculation.
• was not used for hyperparameter estimation.

As you have only few cases, a resampling strategy would be appropriate. In this situation, it is probably best to fix any hyperparameters (like number of PCs or PLS latent variables, or the LASSO bound) by external knowledge in order to avoid a second, inner split of your training data for the optimization of the hyperparameter.