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
- to what extent this is reasonable / justifiable, and
- 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
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,
- 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
- 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.