Estoy interesado en determinar el número de patrones significativos que salen de un análisis de componentes principales (PCA) o de un análisis de funciones ortogonales empíricas (EOF). Estoy especialmente interesado en aplicar este método a los datos climáticos. El campo de datos es una matriz MxN, siendo M la dimensión temporal (por ejemplo, días) y N la dimensión espacial (por ejemplo, ubicaciones lon/lat). He leído sobre un posible método bootstrap para determinar los PCs significativos, pero no he podido encontrar una descripción más detallada. Hasta ahora, he estado aplicando la regla de oro de North (North y otros ., 1982) para determinar este límite, pero me preguntaba si existía un método más sólido.
Como ejemplo:
###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$values
#projection of data onto EOFs (U) to derive principle components (A)
A <- Fsc %*% U
dim(U)
dim(A)
#plot of top 10 Lambda
plot(L[1:10], log="y")
#plot of explained variance (explvar, %) by each EOF
explvar <- L/sum(L) * 100
plot(explvar[1:20], log="y")
#plot original patterns versus those identified by EOF
layout(matrix(1:12, nrow=4, ncol=3, byrow=TRUE), widths=c(1,1,1), heights=c(1,0.5,1,0.5))
layout.show(12)
par(mar=c(4,4,3,1))
image(tmp1, main="pattern 1")
image(tmp2, main="pattern 2")
image(tmp3, main="pattern 3")
par(mar=c(4,4,0,1))
plot(T, tp1, t="l", xlab="", ylab="")
plot(T, tp2, t="l", xlab="", ylab="")
plot(T, tp3, t="l", xlab="", ylab="")
par(mar=c(4,4,3,1))
image(matrix(U[,1], length(x), length(y)), main="eof 1")
image(matrix(U[,2], length(x), length(y)), main="eof 2")
image(matrix(U[,3], length(x), length(y)), main="eof 3")
par(mar=c(4,4,0,1))
plot(T, A[,1], t="l", xlab="", ylab="")
plot(T, A[,2], t="l", xlab="", ylab="")
plot(T, A[,3], t="l", xlab="", ylab="")
Y, aquí está el método que he estado utilizando para determinar la importancia del PC. Básicamente, la regla general es que la diferencia entre lambdas vecinas debe ser mayor que su error asociado.
###Determine significant EOFs
#North's Rule of Thumb
Lambda_err <- sqrt(2/dim(F)[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
plot(L[1:10],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)
He encontrado la sección del capítulo de Björnsson y Venegas ( 1997 ) sobre las pruebas de significación para ser útil - se refieren a tres categorías de pruebas, de las cuales la varianza dominante -tipo es probablemente lo que espero utilizar. Se refieren a un tipo de enfoque de Monte Carlo que consiste en barajar la dimensión temporal y volver a calcular las lambdas en muchas permutaciones. von Storch y Zweiers (1999) también se refieren a una prueba que compara el espectro lambda con un espectro de "ruido" de referencia. En ambos casos, no estoy seguro de cómo se puede hacer esto, y también de cómo se hace la prueba de significación dados los intervalos de confianza identificados por las permutaciones.
Gracias por su ayuda.
Referencias: Björnsson, H. y Venegas, S.A. (1997). "A manual for EOF and SVD analyses of climate data", McGill University, CCGCR Report No. 97-1, Montréal, Québec, 52pp. http://andvari.vedur.is/%7Efolk/halldor/PICKUP/eof.pdf
G.R. North, T.L. Bell, R.F. Cahalan y F.J. Moeng. (1982). Sampling errors in the estimation of empirical orthogonal functions. Mon. Wea. Rev., 110:699-706.
von Storch, H, Zwiers, F.W. (1999). Statistical analysis in climate research. Cambridge University Press.
