13
$\begingroup$

I'm a biologist, and I have a large dataset that I'm trying to analyze. Here are the variables I'm working with:

  • levels of 211 different metabolites in 16 different blood samples (predictor variables)
  • how well each of the 16 blood samples performed in a specific test (response variable)

I am trying to figure out which variable(s) are the most important in predicting the performance of a blood sample in the test.

I would like to make these data more manageable by doing a PCA, but I'm new to this sort of analysis. I understand that a PCA will create groups of principal components (it will group metabolites that covary with each other and label each of these groups a principal component), but I'm not sure how to take into account the response variable in this analysis.

Any help would be much appreciated!!! Thank you.

> summary(metabolites_princomp)
Importance of components:
                          Comp.1     Comp.2      Comp.3      Comp.4
Standard deviation     3.9608225 0.40128486 0.259868774 0.215004349
Proportion of Variance 0.9805072 0.01006435 0.004220736 0.002889179
Cumulative Proportion  0.9805072 0.99057153 0.994792267 0.997681446
[...]
$\endgroup$
16
  • $\begingroup$ What language or software are you planning to work with? (It doesn't make a difference as to what PCA means but it could help shape a response) $\endgroup$ Jul 26, 2012 at 1:20
  • $\begingroup$ I've been working in R so far, but I'm new to it. I've done a PCA on the metabolite information alone (not taking into account the response variable) and I got 16 principal components, the first of which accounts for 98% of the variance. $\endgroup$
    – Lindsay
    Jul 26, 2012 at 1:29
  • $\begingroup$ I can already tell you you forgot to center the rows of your matrix! $\endgroup$ Jul 26, 2012 at 1:43
  • 1
    $\begingroup$ @David: X - rowMeans(X) will be quite a bit more efficient if you want to center the rows. :-) $\endgroup$
    – cardinal
    Jul 26, 2012 at 2:05
  • 2
    $\begingroup$ @Lindsay: Unfortunately, it is not the rows you want to be centering in your case. If you want things standardized (centered and rescaled) you can just do princomp(X,cor=T). $\endgroup$
    – cardinal
    Jul 26, 2012 at 2:14

2 Answers 2

10
$\begingroup$

First, note that your understanding of PCA is slightly off. PCA doesn't "group" variables into principal components. Each principal component is, rather, a new variable (a "new metabolite") of the same length as each of your original variables (in your case, 16).

It's true that each principal component can represent a group of original metabolites, but it's more complicated and sophisticated than that. Each metabolite can actually be correlated with multiple principal components, each a different amount. It could be true that a group of metabolites are purely represented principal component A and another group is purely represented by principal component B, but it is far more likely that each is correlated with multiple PCs, each to a different extent.

So think of each PC as a new metabolite (in PCA parlance we'd call it an "eigenmetabolite") that is representative of a general pattern across the metabolites. Some metabolites are highly representative of this pattern, while some aren't, and some share some properties with multiple patterns.

Now we come to your question. The way we bring in a response variable is to look at the relationship- for example, a correlation or linear regression- between each principal component (each eigenmetabolite) and the response variable. Here is a neat little example in R code. (Note that I use the svd function rather than the pls package, as you might be using).

set.seed(1337)

center.rows = function(m) t(apply(m, 1, function(r) r - mean(r)))

response.var = rnorm(16, 0, 1)

relevant.metabolites = t(replicate(50, rnorm(16, response.var, 2)))
unrelated.metabolites = t(replicate(110, rnorm(16, 0, 2)))
metabolite.data = rbind(relevant.metabolites, unrelated.metabolites)

s = svd(center.rows(metabolite.data))

What I did is create a matrix of simulated metabolite data with 50 genes that are somewhat correlated with the response variable, and 110 genes that aren't. We can see that only one PC is significant, and it explains ~33% of the variance:

var.explained = s$d^2 / sum(s$d^2)
plot(var.explained)

enter image description here

Now let's compare that PC to the response variable:

plot(response.var, s$v[, 1])

enter image description here

As you can see, the first PC almost perfectly recaptures the response variable that we used to create the matrix (correlation is -.993), even though we put a lot of extra noise in. (Don't worry that the correlation is negative- the direction is arbitrary). Indeed, this correlation with the "eigenmetabolite" is considerably greater than the correlation of the response variable to any of the individual significant genes:

hist(apply(relevant.metabolites, 1, function(r) cor(r, response.var)))

enter image description here

Of course, that's exactly why we're applying PCA in the first place. While the relationship between the response variable and any individual metabolite may be too weak to measure over the rest of the noise, when you combine many genes together you might find that a principal component captures that response well!

$\endgroup$
13
  • $\begingroup$ @ David: Thank you so much for such a detailed explanation! I'm going to give it a go right now. $\endgroup$
    – Lindsay
    Jul 26, 2012 at 2:25
  • 1
    $\begingroup$ It doesn't look like your "unrelated" metabolites are unrelated. :-) One aspect of this answer that may be misleading is that even if there were some PC associated with the response, it may not be among the ones explaining the variation in the data. So, some caution is advised. $\endgroup$
    – cardinal
    Jul 26, 2012 at 2:26
  • $\begingroup$ @David: I'm getting stuck at the step where you plot PC1 against the response variable. I have a data.frame with 212 rows and 16 columns. The last row is my response variable, so I'm typing in plot(metabolites[212,], s$v[,1]). Why won't that work for me? $\endgroup$
    – Lindsay
    Jul 26, 2012 at 2:41
  • $\begingroup$ What is your error? $\endgroup$ Jul 26, 2012 at 2:42
  • 1
    $\begingroup$ @cardinal: right you are, that was a typo, since fixed (I'm not at a computer so I'll have to adjust the figures later in response to that change). Lindsay, what are the dimensions of your s$v matrix? $\endgroup$ Jul 26, 2012 at 2:47
2
$\begingroup$

After the great discussion, above, I'm hesitant to add another idea, but it seems to me that partial least squares may do what you want, more easily. I usually do this in SAS, so I don't know the details of doing it in R, but there is a PLS package that does both partial least squares regression and principal component regression.

The essential idea behind PLS is that it incorporates both relations among the independent variables and between the independent and dependent variables.

$\endgroup$
5
  • $\begingroup$ Thanks for your input. I've been doing a bit of research about this today, and I was drifting towards PCR as well. I'm not sure that I quite understand PLS and PCR yet enough to be confident about executing them, but it does sound like it might make more sense considering my goals. $\endgroup$
    – Lindsay
    Jul 26, 2012 at 19:52
  • $\begingroup$ Are 'partial least squares' and 'squared partial correlation' the same thing? I'm wondering whether this is similar to Kruskal's Drivers Analysis $\endgroup$ Apr 20, 2017 at 20:28
  • $\begingroup$ @MarkRamotowski I dont think they are teh same, but I never heard of Kruskal's drivers analysis. $\endgroup$
    – Peter Flom
    Apr 21, 2017 at 11:34
  • $\begingroup$ It's used in market research - I'm not sure how valid an approach it is if you've never heard of it! wiki.q-researchsoftware.com/wiki/… $\endgroup$ Apr 21, 2017 at 12:21
  • $\begingroup$ @MarkRamotowski There's a LOT of stuff I haven't heard of! I don't do market research, so, my ignorance of its methods doesn't mean much. $\endgroup$
    – Peter Flom
    Apr 22, 2017 at 12:52

Not the answer you're looking for? Browse other questions tagged or ask your own question.