I am trying to determine the correct amount of variance explained by each mode of an Empirical Orthogonal Function (EOF) analysis (similar to "PCA") as applied to a gappy data set. (i.e., containing NaNs). The following question builds on an earlier one that I had regarding the differing results obtained from the decomposition of the data set's covariance matrix using either eigen or svd. In essence, the problem is that I have read that both decompositions can be used interchangeably for obtaining the EOFs from a square covariance matrix. This does seem to be the case when the data set is not gappy (as illustrated below):

###Make complete and gappy data set
x <- 1:100
y <- 1:100
grd <- expand.grid(x=x, y=y)

#complete data
z <- matrix(rnorm(dim(grd)[1]), length(x), length(y))
image(x,y,z, col=rainbow(100))

#gappy data
zg <- replace(z, sample(seq(z), length(z)*0.5), NaN)
image(x,y,zg, col=rainbow(100))

###Covariance matrix decomposition
#complete data
C <- cov(scale(z), use="pair")
E <- eigen(C)
S <- svd(C)

#sum of lambda

The sum of lambda in both eigen and svd equals the sum of the diagonal of the covariance matrix. So far, so good - Both methods explain the correct amount of variance. The next example does the same routine for a gappy version of the data set (50% NaNs):

#gappy data (50%)
Cg <- cov(scale(zg), use="pair")
Eg <- eigen(Cg)
Sg <- svd(Cg)

#sum of lambda

And here we see that the lambda values calculated by svd are greater than the sum of the diagonal of the covariance matrix. Those calculated by eigen are equal. However, because the covariance matrix is no longer positive definite, there are some negative trailing lambda values. In my previous question I showed that this tendency becomes greater with increasing gappiness.
So, I can live with this if need be, but now I'm concerned about how to correctly assign how much of the data set's variance is explained by each EOF. This should be lambda/sum(lambda). When I plot the cumulative explained variance of the EOFs, you will see the problem - because the eigen decomposition contains some negative eigenvalues, the slope of cumulative explained variance is steeper and bell-shaped:

#cumulative explained variance of the EOFs
E.cumexplvar <- cumsum(E$values/sum(E$values))
S.cumexplvar <- cumsum(S$d/sum(S$d))
Eg.cumexplvar <- cumsum(Eg$values/sum(Eg$values))
Sg.cumexplvar <- cumsum(Sg$d/sum(Sg$d))

###plot the cumulative explained variance
png("cumexplvar.png", width=8, height=4, units="in", res=200)
YLIM <- range(c(E.cumexplvar, S.cumexplvar, Eg.cumexplvar, Sg.cumexplvar))

plot(E.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="non-gappy")
points(S.cumexplvar, t="o", pch=2, col=2)
abline(h=1, col=8, lty=2)
legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)

plot(Eg.cumexplvar, t="o", col=1, ylim=YLIM, xlab="EOF", ylab="cum. expl. var.", main="gappy")
points(Sg.cumexplvar, t="o", pch=2, col=2)
abline(h=1, col=8, lty=2)
legend("bottomright", legend=c("Eigen", "SVD"), col=c(1,2), pch=c(1,2), lty=1)


cumulative explained variance

The problem may be that I should be using the sum of the absolute eigenvalues to assign their explained variance, but this also leaves me to wonder how to interpret the explained variance of the negative eigenvalues. I would be very grateful for any insight, as this is not an issue that I have come across in any reference regarding EOF as applied to gappy data.

  • 1
    $\begingroup$ As I've commented on your previous question, you shouldn't believe SVD when your cov matrix is not positive simidefinite. Because negative eigenvalues (returned by SVD as positive) do matter, still, they can't be interpreted as explaining real variance; rather, they explain the extent to which the data do not converge geometrically in euclidean space. The overal variance (trace of cov. matrix) is split into real, substantial variance (sum of positive eigenvalues) and "imaginary" anti-variance (sum of negative eigenvalues) $\endgroup$ – ttnphns Aug 22 '12 at 9:28
  • $\begingroup$ So, here is no single notion of "true" variance. You have to choose either to count overall ambivalent variance as "true" or to count just positive inflated variance as "true". That depends on your tasks and situation. $\endgroup$ – ttnphns Aug 22 '12 at 9:44
  • $\begingroup$ Thanks @ttnphns. I meant to add a comment to the other question to link to this one - will do that now. I would basically like to illustrate the possible pitfalls of EOF truncation based on a predefined cumulative explained variance (i.e. 80%) rather than some other measure of EOF significance (i.e. "Rule N" or "North's Rule of Thumb"). I hoping to identify if the real part of the eigen decomposition is actually more accurate in this respect. To put it in terms of data reconstruction - decomposition's lambda values accurately characterize the amount of data encoded? $\endgroup$ – Marc in the box Aug 22 '12 at 9:47
  • 1
    $\begingroup$ To put it in terms of data reconstruction... My own experience with reconstructing original data by means of Principal Coordinates Analysis (a version of linear MDS based on PCA idea) from their twisted (noised) cov matrix tells that setting negative eigenvalues to 0 and proportional alteration of positive ones so that they sum to the trace of the matrix gives most accurate reconstruction. But your problem may be different than mine was. $\endgroup$ – ttnphns Aug 22 '12 at 9:59
  • $\begingroup$ Great advise - I will look into this in more detail. Thanks for all your help over the past couple days. $\endgroup$ – Marc in the box Aug 22 '12 at 10:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.