Selecting PCA components which separate groups

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.

• I'd want to make just two comments about your interesting question. 1) Describe in words your approach besides showing its code (remember that people here use various software, not necessarily R). 2) The scatterplot isn't very much convincing without spikes to its floor. Also, if you have any specific doubts about your approach, please speak them out to make the question more focused. Sep 27, 2013 at 11:00
• edited: To clarify, you're first conducting PCA and then trying to isolate the principal components that are best explained by some particular variable? Are you cross-validating these with a scree plot? It might be the case that some x you pick out of your data set happens to explain a lot of the variance in a principal component, but I'm not sure that means anything if the variance is very low along that principal component. Dec 19, 2013 at 12:55
• @ssdecontrol Well, I am doing a poor mans randomization to see whether the calculated $R^2$ is well above background noise. As for whether it means anything -- the point is that it usually does, actually. Since I am doing all sets of classical analyses anyway and / or supervised machine learning, whenever I see that PCX is explained in a significant portion by a classifier, I will (i) find several variables that differ between the groups of this classifier and (ii) that I can successfully train an SML. Dec 20, 2013 at 11:23
• to find out what share of the overall variance in the data matrix is explained by a given classification If want to know just this you need no PCA. Just compute the proportion of the between-group sum-of-squares to the total sum-of-squares: (SStotal-SSwithin)/SStotal where SSwithin is the pooled within-group sum-of-squares. Dec 22, 2013 at 3:07
• I don't see any problems with the way you exploit PCA, but I don't understand why you really need it. (Just because you like it perhaps?) Because I can't see your precise aim I can't tell anything for your Is there anything better?. Dec 22, 2013 at 3:13

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.

• I have started a new bounty to award this answer. Jan 7, 2014 at 10:58
• @January: This is a nice answer, but I would like to point out that "direct LDA" is a very strange algorithm at best, see Gao and Davis, 2005, Why direct LDA is not equivalent to LDA: "we demonstrate that... D-LDA may impose a significant performance limitation in general applications", so be careful with it. Aug 28, 2014 at 13:03
• @amoeba Thanks for that citation. I have had concerns about DLDA for a while, since there is no justification for selecting components in that particular way. I see it as a very problem-specific solution which does not necessarily generalize beyond face discrimination problems, although it can be easily adapted to any problem with some knowledge of which components are most useful for discrimination. Every solution implementing high-dimensional discrimination with an assumed factor structure suffers from problems... have you found any better approaches? I'm interested in your opinion here. Aug 28, 2014 at 13:13
• @ahfoss: I am no specialist here, but as much as I can understand, "direct LDA" is one of the attempts to apply LDA in a $n \ll k$ case. A more standard approach would be to use regularized LDA (rLDA) with cross-validation to find optimal regularization parameter; I have no idea, though, how it actually performs in face recognition. Aug 28, 2014 at 13:28

@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.

• +1 for cross-validating models. Extremely important. However, I would like to hear from OP @January, who has stated he isn't interested in discrimination, although his/her problem seems very well suited to discrimination/classification analysis. Jan 8, 2014 at 15:54
• I disagree with your assertion that k-means/PCA/etc belong to the same family. This implies that they are special cases of the same model or algorithm, which is not true. The PCA algorithm is a simple matrix calculation, whereas k-means is an iterative algorithm that has been compared to the EM algorithm (not technically correct since there is no likelihood function, but still a useful comparison in some respects IMHO). Jan 8, 2014 at 15:58
• Are you referring to the plsgenomics::pls.lda function? If not how is your package different/improved? I will also point out to interested readers that PLS-LDA in general is superior to a commonly used technique of simply running PLS with a dummy-coded outcome variable. While this latter approach is not necessarily incorrect, it is definitely kludgey, not least because you can obtain predicted probabilities less than zero or greater than one! Jan 8, 2014 at 16:08
• @ahfoss: I did not mean the algorithms, as the same underlying model could be calculated by different algorithms. E.g. for PCA you can use iterative (NIPALS, POWER) or non-iterative (EVD, SVD) algorithms. Maybe a better term instead of k-means would have been "cluster analysis minimizing within-cluster sum of squares, for which e.g. k-means is a heuristic approximation". I don't have time now, will look over the answer later or we could meet in the chat room and find a better description. Jan 8, 2014 at 16:18
• ... A technical difference is that I use pls::plsr for the pls (which allows to choose from different algorithms). And I have a bunch of post-processing functions e.g. for flipping and rotating the model, which is sometimes useful for interpretation. Jan 8, 2014 at 16:28