# Choosing number of PCA components when multiple samples for each data point are available

Update: I posted an answer below describing my current attempts to approach this problem.

I am facing a perennial problem of identifying significant principal components. This question has been discussed here a couple of times (one, two and a very impressive thread three), and the usual advice is to use either so called "parallel analysis" (a very misleading name in my opinion; it is essentially Monte Carlo estimation of eigenvalue spectrum under null) or Laplace/BIC approaches from Minka, 2000, Automatic Choice of Dimensionality for PCA. Apart from that, I am aware of some ways to perform cross-validation over data points (see Bro et al., 2008, Cross-validation of component models), plus there is of course a good old approach of staring at the eigenvalue spectrum trying to locate an "elbow".

All of these methods assume that you have one set of data points, let us say $N$ points from $\mathbb{R}^D$, and nothing apart from that. In my case, however, each point is an average of $M$ points, which are repeated measurements of the same value. This makes me think that I can somehow infer which PCA components are significant, by using the data that goes into averages; maybe by some cross-validation or bootstrapping approach. My question is: how?

Update. Following discussion in the comments, I should clarify what I mean by "significant". For this I will introduce some notation. There are $N$ points $\bar{\mathbf{x}}_i \in \mathbb{R}^D$, and each point is an average over $M$ repeated measurements $\bar{\mathbf{x}}_i=\frac{1}{M}\sum_{j=1}^M \mathbf{x}_i^{(j)}$. Each measurement is a noisy observation of a true position of the corresponding point, i.e. $\mathbf{x}_i^{(j)}=\mathbf{a}_i + \mathbf{\xi}_i^{(j)}$ (note that $\xi$ is a $D$-dimensional vector as well, but I cannot make Greek letters bold here). I am interested in principal components of $\{\mathbf{a}_i\}$, but can only perform PCA on $\{\bar{\mathbf{x}}_i\}$, and the small components can be completely distorted by the noise. The larger the noise, the less principal components of $\{\mathbf{a}_i\}$ I will be able to recover. The challenge is to select only those leading PCs that definitely "come from" $\{\mathbf{a}_i\}$.

In my answer below I describe two approaches that I currently use. I am not fully happy with either.

• I upvoted because of the interesting references you cited. That's a good way to frame a question. Nice. Jan 29, 2014 at 18:56
• The meaning of "significant" is usually not in question because a probability model and null hypothesis are implied by the context. Not so here. The replication reveals effects of measurement error on the PCA results. However, if you view your $N$ points as being a single realization of some process, then "significant" might mean with respect to a null hypothesis like "all components are equal." (This issue is rarely discussed for PCA because it's usually used as an exploratory rather than a confirmatory analysis.) So: what is your probability model and what hypothesis are you testing?
– whuber
Jan 29, 2014 at 22:19
• @whuber: This is a fair question. I do use PCA as an exploratory analysis, and am not testing any hypothesis. I assume that my $N$ points come from an underlying process, but are corrupted by some noise. In the limit of infinite number of measurements (each point $\mathbf{f}(t_i)$ is observed a lot of times, $M \to \infty$) I would say that all my components are "significant" (i.e. represent the actual signal), no matter how small. In reality, small components of the signal will not be observed because of the noise. Question is: how to separate signal components from the noise components? Jan 29, 2014 at 22:38
• That is an interesting question--but it is not what you seem to be asking here. One usually performs the separation in the obvious way: by splitting each set of multiple observations into their mean plus their residuals. The latter estimate the noise while the former estimates the points.
– whuber
Jan 30, 2014 at 16:41
• I do not follow your terminology where you use the words "represent" and "signal." However, I think I get the gist of it: it will be difficult to identify components whose variance is smaller than the noise variance. But that is a criterion right there! You have an estimate of the noise variance from the residuals and you have estimates of the component variances from the squares of the eigenvalues. How you determine which components are "significantly" greater than the noise ought to depend on what you are trying to accomplish and what the consequences are of retaining too many or too few.
– whuber
Jan 30, 2014 at 17:26

I'm not sure that I completely understand the data. Do you have M replicates at N timepoints of a response in $\Re^d$? That would make $N \times M \times D$ actual numbers? Are the trajectories vector valued functions?

If I have understood this correctly, I don't think that your replicates are going to tell you anything about the number of PC's. The replicates enable you to obtain a more accurate estimate of the covariance matrix than you would otherwise have, but in themselves, they don't tell you anything about the relationship between successive time points (which is what the eigen-decomposition is about).

To determine significance, you need some sort of probability model to the effect that the data belong to something in a reduced space (the PC space) with full rank error. Factor analysis posits models of this kind and the factor analysis literature is full of pointers on how to choose the right number of factors (it's as much an art as a science, I believe). Jeremy Anglim's blog is a useful resource, especially This post

• To your questions: yes, yes, and yes, so I think you understood correctly (maybe I should add some more explicit notation to my question). However, I think I disagree. Here is some intuition: I can take 1st out of $M$ replications for each point and run PCA on the resulting $N$ points; then I can take the 2nd, 3rd and so on. I will obtain $M$ PCA decompositions. Now if a component is significant, I would expect it to appear in each of those PCA replications. But if it is noise, then it will look different every time. Jan 29, 2014 at 21:30
• @amoeba Yes, that is certainly true. You would be using your reps to validate the model. The model being so validated would be the usual sort of model, however. You would still need some sort of criterion for determining fitness. You could do exploratory FA on half the data and a confirmatory on the other half, I suppose. Jan 29, 2014 at 22:58
• This suggestion sounds interesting, I will look into it. Jan 30, 2014 at 10:48
• I updated my question with some more explicit notation and clarifications and posted an answer describing my two current approaches. It seems to me that my approach #2 is related to your suggestion about confirmatory FA. Feb 5, 2014 at 23:59

In order to advance the discussion here I will describe two approaches that I am currently using.

For convenience I will repeat the notation here. There are $N$ points $\bar{\mathbf{x}}_i \in \mathbb{R}^D$, and each point is an average over $M$ repeated measurements $\bar{\mathbf{x}}_i=\frac{1}{M}\sum_{j=1}^M \mathbf{x}_i^{(j)}$. Each measurement is a noisy observation of a true position of the corresponding point, i.e. $\mathbf{x}_i^{(j)}=\mathbf{a}_i + \mathbf{\xi}_i^{(j)}$ (note that $\xi$ is a $D$-dimensional vector as well, but I cannot make Greek letters bold here). I am interested in principal components of $\{\mathbf{a}_i\}$, but can only perform PCA on $\{\bar{\mathbf{x}}_i\}$, and the small components can be completely distorted by the noise. The larger the noise, the less principal components of $\{\mathbf{a}_i\}$ I will be able to recover. The challenge is to select only those leading components that definitely "come from" $\{\mathbf{a}_i\}$.

First method

If the measurement noise in one particular dimension has variance $\sigma^2$ (which can be different for different dimensions and data points), then the standard error of the mean is given by $\frac{\sigma^2}{M}$. Informally, this is the amount of noise that will remain after averaging and can potentially screw the principal components. I want to estimate how much variance can possibly come from this remaining noise. Note that if I subtract two measurements $\Delta \mathbf{x}_i=\mathbf{x}_i^{(j_1)}-\mathbf{x}_i^{(j_2)}$, I get rid of the signal $\mathbf{a}_i$ and get twice the noise variance. So to estimate the maximum amount of overall "noise variance" I take $\frac{1}{\sqrt{2M}}\Delta \mathbf{x}_i$ as data points, run PCA on them and take the variance of the first PC as my noise floor. All PCs from $\bar{\mathbf{x}}_i$ that have variance above this noise floor I declare significant. See a figure below.

In fact I can take different pairs of measurements to construct my noise estimates, and this will give slightly different noise floors. I run this procedure multiple times, randomly selecting pairs each time, and then take the average noise floor over repetitions.

This procedure seems to follow @whuber's suggestion in the comment above.

Second method

This is a cross-validation procedure. I split my $M$ measurements in two parts, and average them separately, obtaining training data $\bar{\mathbf{x}}_i$ and test data $\bar{\mathbf{y}}_i$. Then I compute PCA on $\bar{\mathbf{x}}_i$ and get the projection of my data onto the principal axes: $\mathbf{z}_i$. The question is now: how many dimensions of $\mathbf{z}$ should I take to minimize the reconstruction error of the test dataset? More precisely, for each number $k$ of components I project the training data onto the $k$-dimensional subspace spanned by the first $k$ principal axes $\mathbf{z}_i=WW^\top \bar{\mathbf{x}}_i$ where $W$ is the $D\times k$ matrix of the first $k$ principal axes, and compute the reconstruction error $e(k)=\sum_i||\bar{\mathbf{y}}_i-\mathbf{z}_i||^2$. This error should have a minimum at certain $k$ and this number of components I declare significant.

Here again I repeat this procedure many times for different random splits of the data, average obtained reconstruction errors $e(k)$ and then find the minimum.

This procedure is probably related to @Placidia's suggestion ("do exploratory FA on half the data and a confirmatory on the other half"), but I know too little about confirmatory factor analysis to say if it's really the same thing.

Results

The following figure shows the outcomes of both methods on four different datasets. Four columns are datasets, top row is noise floor method, second row is cross-validation method. On the top, red line shows the empirical eigenvalue spectrum, blue lines are noise spectra from different repetitions, and the dashed horizontal line shows the mean noise floor. Black dots mark components above the noise floor.

On the bottom, blue lines are reconstruction errors on various cross-validation folds, red line is the mean, black dots mark the components until the minimum of the red curve.

The two methods give similar, but not identical results. In the fourth dataset I get 4 components with each method. In the second and third datasets cross-validation results in a bit more significant components than the noise floor. However, in the first dataset there is one component well above the noise floor, but cross-validation shows that minimum reconstruction error is obtained with zero components.

Both methods seem conservative to me, however I don't like that sometimes one and sometimes another turns out to be more sensitive. At the moment I compute both estimates and take the maximum; this is reasonable if both estimates are indeed conservative, but not very elegant.