How to determine significant principal components using bootstrapping or Monte Carlo approach?

I am interested in determining the number of significant patterns coming out of a Principal Component Analysis (PCA) or Empirical Orthogonal Function (EOF) Analysis. I am particularly interested in applying this method to climate data. The data field is a MxN matrix with M being the time dimension (e.g. days) and N being the spatial dimension (e.g. lon/lat locations). I have read of a possible bootstrap method to determine significant PCs, but have been unable to find a more detailed description. Until now, I have been applying North's Rule of Thumb (North et al., 1982) to determine this cutoff, but I was wondering if a more robust method was available.

As an example:

###Generate data
x <- -10:10
y <- -10:10
grd <- expand.grid(x=x, y=y)

#3 spatial patterns
sp1 <- grd$x^3+grd$y^2
tmp1 <- matrix(sp1, length(x), length(y))
image(x,y,tmp1)

sp2 <- grd$x^2+grd$y^2
tmp2 <- matrix(sp2, length(x), length(y))
image(x,y,tmp2)

sp3 <- 10*grd$y tmp3 <- matrix(sp3, length(x), length(y)) image(x,y,tmp3) #3 respective temporal patterns T <- 1:1000 tp1 <- scale(sin(seq(0,5*pi,,length(T)))) plot(tp1, t="l") tp2 <- scale(sin(seq(0,3*pi,,length(T))) + cos(seq(1,6*pi,,length(T)))) plot(tp2, t="l") tp3 <- scale(sin(seq(0,pi,,length(T))) - 0.2*cos(seq(1,10*pi,,length(T)))) plot(tp3, t="l") #make data field - time series for each spatial grid (spatial pattern multiplied by temporal pattern plus error) set.seed(1) F <- as.matrix(tp1) %*% t(as.matrix(sp1)) + as.matrix(tp2) %*% t(as.matrix(sp2)) + as.matrix(tp3) %*% t(as.matrix(sp3)) + matrix(rnorm(length(T)*dim(grd)[1], mean=0, sd=200), nrow=length(T), ncol=dim(grd)[1]) # error term dim(F) image(F) ###Empirical Orthogonal Function (EOF) Analysis #scale field Fsc <- scale(F, center=TRUE, scale=FALSE) #make covariance matrix C <- cov(Fsc) image(C) #Eigen decomposition E <- eigen(C) #EOFs (U) and associated Lambda (L) U <- E$vectors
L <- E$d #projection of data onto EOFs (U) to derive principle components (A) A <- Xp %*% U  Following the example that I used in my original example, I will determine "significant" EOFs via North's rule of thumb. North's Rule of Thumb Lambda_err <- sqrt(2/dim(Xp)[2])*L upper.lim <- L+Lambda_err lower.lim <- L-Lambda_err NORTHok=0*L for(i in seq(L)){ Lambdas <- L Lambdas[i] <- NaN nearest <- which.min(abs(L[i]-Lambdas)) if(nearest > i){ if(lower.lim[i] > upper.lim[nearest]) NORTHok[i] <- 1 } if(nearest < i){ if(upper.lim[i] < lower.lim[nearest]) NORTHok[i] <- 1 } } n_sig <- min(which(NORTHok==0))-1 n_sig plot(L[1:20],log="y", ylab="Lambda (dots) and error (vertical lines)", xlab="EOF") segments(x0=seq(L), y0=L-Lambda_err, x1=seq(L), y1=L+Lambda_err) abline(v=n_sig+0.5, col=2, lty=2) text(x=n_sig, y=mean(L[1:10]), labels="North's Rule of Thumb", srt=90, col=2)  Since the Lambda values of 2:4 are very close to each other in amplitude, these are deemed insignificant by the rule of thumb - i.e. their respective EOF patterns might overlap and mix given their similar amplitudes. This is unfortunate given that we know that 9 signals actually exist in the field. A more subjective approach is to view the log-transformed Lambda values ("Scree plot") and to then fit a regression to the trailing values. One then can determine visually at what level the lambda values lie above this line: Scree plot ntrail <- 35 tail(L, ntrail) fit <- lm(log(tail(L, ntrail)) ~ seq(length(L)-ntrail+1, length(L))) plot(log(L)) abline(fit, col=2)  So, the 5 leading EOFs lie above this line. I have tried this example when Xp has no additional noise added and the results reveal all 9 original signals. So, the insignificance of EOFs 6:9 is due to the fact that their amplitude is lower than the noise in the field. A more objective method is the "Rule N" criteria by Overland and Preisendorfer (1982). There is an implementation within the wq package, which I show below. Rule N library(wq) eofNum(Xp, distr = "normal", reps = 99) RN <- ruleN(nrow(Xp), ncol(Xp), type = "normal", reps = 99) RN eigs <- svd(cov(Xp))$d
plot(eigs, log="y")
lines(RN, col=2, lty=2)


The Rule N identified 4 significant EOFs. Personally, I need to better understand this method; Why is it possible to gauge the level of error based on a random field that does not use the same distribution as that in Xp? One variation on this method would be to resample the data in Xp so that each column is reshuffled randomly. In this way, we ensure that the total variance of the random field is the same as Xp. By resampling many times, we are then able to calculate a baseline error of the decomposition.

Monte Carlo with random field (i.e. Null model comparison)

iter <- 499
LAMBDA <- matrix(NaN, ncol=iter, nrow=dim(Xp)[2])

set.seed(1)
for(i in seq(iter)){
#i=1

#random reorganize dimensions of scaled field
Xp.tmp <- NaN*Xp
for(j in seq(dim(Xp.tmp)[2])){
#j=1
Xp.tmp[,j] <- Xp[,j][sample(nrow(Xp))]
}

#make covariance matrix
C.tmp <- t(Xp.tmp) %*% Xp.tmp #cov(Xp.tmp)

#SVD decomposition
E.tmp <- svd(C.tmp)

#record Lambda (L)
LAMBDA[,i] <- E.tmp$d print(paste(round(i/iter*100), "%", " completed", sep="")) } boxplot(t(LAMBDA), log="y", col=8, border=2, outpch="") points(L)  Again, 4 EOFs are above the distributions for the random fields. My worry with this approach, and that of Rule N, is that these are not truly addressing the confidence intervals of the Lambda values; e.g. a high first Lambda value will automatically result in a lower amount of variance to be explained by trailing ones. Thus the Lambda calculated from random fields will always have a less steep slope and may result in selecting too few significant EOFs. [NOTE: The eofNum() function assumes that EOFs are calculated from a correlation matrix. This number might be different if using a e.g. a covariance matrix (centered but not scaled data).] Finally, @tomasz74 mentioned the paper by Babamoradi et al. (2013), which I have had a brief look at. Its very interesting, but seems to be more focused on calculating CI's of EOF loadings and coefficients, rather than Lambda. Nevertheless, I believe that it might be adopted to assess Lambda error using the same methodology. A bootstrap resampling is done of the data field by resampling the rows until a new field is produced. The same row can be resampled more than once, which is a non-parametric approach and one does not need to make assumptions about the distribution of data. Bootstrap of Lambda values B <- 40 * nrow(Xp) LAMBDA <- matrix(NaN, nrow=length(L), ncol=B) for(b in seq(B)){ samp.b <- NaN*seq(nrow(Xp)) for(i in seq(nrow(Xp))){ samp.b[i] <- sample(nrow(Xp), 1) } Xp.b <- Xp[samp.b,] C.b <- t(Xp.b) %*% Xp.b E.b <- svd(C.b) LAMBDA[,b] <- E.b$d
print(paste(round(b/B*100), "%", " completed", sep=""))
}
boxplot(t(LAMBDA), log="y", col=8, outpch="", ylab="Lambda [log-scale]")
points(L, col=4)
legend("topright", legend=c("Original"), pch=1, col=4)


While this may be a more robust than North's rule of thumb for calculating the error of Lambda values, I believe now that the question of EOF significance comes down to different opinions on what this means. For the North's rule of thumb and bootstrap methods, significance appears to be more based on whether or not teere is overlap between Lambda values. If there is, then these EOFs may be mixed in their signals and not represent "true" patterns. On the other hand, these two EOFs may describe a significant amount of variance (as compared to the decomposition of a random field - e.g. Rule N). So if one is interested in filtering out noise (i.e via EOF truncation) then Rule N would be sufficient. If one is interested in determining real patterns in a data set, then the more stringent criteria of Lambda overlap may be more robust.

Again, I am not an expert in these issues, so I am still hoping that someone with more experience can add to this explanation.

References:

Beckers, Jean-Marie, and M. Rixen. "EOF Calculations and Data Filling from Incomplete Oceanographic Datasets." Journal of Atmospheric and Oceanic Technology 20.12 (2003): 1839-1856.

Overland, J., and R. Preisendorfer, A significance test for principal components applied to a cyclone climatology, Mon. Wea. Rev., 110, 1-4, 1982.

• As a remark to > My worry with this approach, and that of Rule N, is that these are not truly addressing the confidence intervals of the Lambda values; e.g. a high first Lambda value will automatically result in a lower amount of variance to be explained by trailing ones. the publication by Daniel Wilks (2016) could be helpful, where he tried to overcome this limitation: doi.org/10.1175/JCLI-D-15-0812.1 Aug 4, 2021 at 12:25