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

Question

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.

Answer 1

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

Answer 2

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.