¿Cómo determinar componentes principales significativos usando bootstrapping o el enfoque de Monte Carlo?

Estoy interesado en determinar la cantidad de patrones significativos que surgen de un Análisis de Componentes Principales (PCA) o un Análisis de Función Ortogonal Empírica (EOF). Estoy particularmente interesado en aplicar este método a los datos climáticos. El campo de datos es una matriz MxN con M siendo la dimensión de tiempo (por ejemplo, días) y N siendo la dimensión espacial (por ejemplo, ubicaciones de lon / lat). He leído sobre un posible método de arranque para determinar PC importantes, pero no he podido encontrar una descripción más detallada. Hasta ahora, he estado aplicando la regla general de North (North et al ., 1982) para determinar este límite, pero me preguntaba si había un método más robusto disponible.

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 usando para determinar la importancia de la PC. Básicamente, la regla general es que la diferencia entre las 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)

Me pareció útil la sección del capítulo de Björnsson y Venegas ( 1997 ) sobre las pruebas de significación: se refieren a tres categorías de pruebas, de las cuales el tipo de varianza dominante es probablemente lo que espero usar. Se refieren a un tipo de enfoque de Monte Carlo de barajar la dimensión del tiempo y volver a calcular las Lambdas en muchas permutaciones. von Storch y Zweiers (1999) también se refieren a la prueba que compara el espectro Lambda con un espectro de "ruido" de referencia. En ambos casos, estoy un poco inseguro de cómo se puede hacer esto, y también de cómo se realiza la prueba de significación dados los intervalos de confianza identificados por las permutaciones.

Gracias por tu ayuda.

Referencias: Björnsson, H. y Venegas, SA (1997). "Un manual para análisis EOF y SVD de datos climáticos", McGill University, CCGCR Report No. 97-1, Montreal, Québec, 52pp. http://andvari.vedur.is/%7Efolk/halldor/PICKUP/eof.pdf

GR North, TL Bell, RF Cahalan y FJ Moeng. (mil novecientos ochenta y dos). Errores de muestreo en la estimación de funciones empíricas ortogonales. Lun. Wea. Rev., 110: 699–706.

von Storch, H, Zwiers, FW (1999). Análisis estadístico en investigación climática. Prensa de la Universidad de Cambridge.

Voy a tratar de avanzar un poco en el diálogo aquí, aunque esta es mi pregunta. Han pasado 6 meses desde que pregunté esto y desafortunadamente no se han dado respuestas completas. Intentaré resumir lo que he reunido hasta ahora y ver si alguien puede dar más detalles sobre los problemas restantes. Disculpe la larga respuesta, pero no veo otra manera ...

Primero, demostraré varios enfoques usando un conjunto de datos sintéticos posiblemente mejor. Proviene de un artículo de Beckers y Rixon ( 2003 ) que ilustra el uso de un algoritmo para realizar EOF en datos breves. He reproducido el algoritmo en R si alguien está interesado ( enlace ).

Conjunto de datos sintéticos:

#color palette
pal <- colorRampPalette(c("blue", "cyan", "yellow", "red"))

#Generate data
m=50
n=100
frac.gaps <- 0.5 # the fraction of data with NaNs
N.S.ratio <- 0.25 # the Noise to Signal ratio for adding noise to data

x <- (seq(m)*2*pi)/m
t <- (seq(n)*2*pi)/n


#True field
Xt <- 
 outer(sin(x), sin(t)) + 
 outer(sin(2.1*x), sin(2.1*t)) + 
 outer(sin(3.1*x), sin(3.1*t)) +
 outer(tanh(x), cos(t)) + 
 outer(tanh(2*x), cos(2.1*t)) + 
 outer(tanh(4*x), cos(0.1*t)) + 
 outer(tanh(2.4*x), cos(1.1*t)) + 
 tanh(outer(x, t, FUN="+")) + 
 tanh(outer(x, 2*t, FUN="+"))

Xt <- t(Xt)
image(Xt, col=pal(100))

#Noise field
set.seed(1)
RAND <- matrix(runif(length(Xt), min=-1, max=1), nrow=nrow(Xt), ncol=ncol(Xt))
R <- RAND * N.S.ratio * Xt

#True field + Noise field
Xp <- Xt + R
image(Xp, col=pal(100))

Entonces, el campo de datos verdadero Xtestá compuesto por 9 señales y le he agregado algo de ruido para crear el campo observado Xp, que se utilizará en los ejemplos a continuación. Los EOF se determinan como tales:

EOF

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

#Eigen decomposition
E <- svd(C)

#EOFs (U) and associated Lambda (L) 
U <- E$u
L <- E$d

#projection of data onto EOFs (U) to derive principle components (A)
A <- Xp %*% U

Siguiendo el ejemplo que utilicé en mi ejemplo original, determinaré los EOF "significativos" mediante la regla general de North.

La regla de oro del norte

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)

Dado que los valores de Lambda de 2: 4 están muy cerca uno del otro en amplitud, la regla general los considera insignificantes, es decir, sus respectivos patrones EOF pueden superponerse y mezclarse debido a sus amplitudes similares. Esto es lamentable dado que sabemos que en realidad existen 9 señales en el campo.

Un enfoque más subjetivo es ver los valores de Lambda transformados logarítmicamente ("diagrama de pantalla") y luego ajustar una regresión a los valores finales. Entonces se puede determinar visualmente a qué nivel los valores lambda se encuentran por encima de esta línea:

Gráfico de sedimentación

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)

Entonces, los 5 EOF principales se encuentran por encima de esta línea. He intentado este ejemplo cuandoXp no se ha agregado ruido adicional y los resultados revelan las 9 señales originales. Entonces, la insignificancia de EOFs 6: 9 se debe al hecho de que su amplitud es menor que el ruido en el campo.

Un método más objetivo es el criterio de la "Regla N" de Overland y Preisendorfer (1982). Hay una implementación dentro delwq paquete, que muestro a continuación.

Regla 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)

La Regla N identificó 4 EOF significativos. Personalmente, necesito entender mejor este método; ¿Por qué es posible medir el nivel de error en función de un campo aleatorio que no utiliza la misma distribución que en Xp? Una variación de este método sería volver a muestrear los datos Xppara que cada columna se reorganice aleatoriamente. De esta manera, nos aseguramos de que la varianza total del campo aleatorio sea la misma que Xp. Al volver a muestrear muchas veces, podemos calcular un error de referencia de la descomposición.

Monte Carlo con campo aleatorio (es decir, comparación de modelo nulo)

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)

Nuevamente, 4 EOF están por encima de las distribuciones para los campos aleatorios. Mi preocupación con este enfoque, y el de la Regla N, es que estos no abordan realmente los intervalos de confianza de los valores de Lambda; por ejemplo, un primer valor Lambda alto dará como resultado automáticamente una menor cantidad de variación que se explicará por los finales. Por lo tanto, el Lambda calculado a partir de campos aleatorios siempre tendrá una pendiente menos pronunciada y puede resultar en la selección de muy pocos EOF significativos. [NOTA: La eofNum()función supone que los EOF se calculan a partir de una matriz de correlación. Este número puede ser diferente si se utiliza, por ejemplo, una matriz de covarianza (datos centrados pero no escalados).]

Finalmente, @ tomasz74 mencionó el artículo de Babamoradi et al. (2013), que he visto brevemente. Es muy interesante, pero parece estar más centrado en el cálculo de CI de cargas y coeficientes EOF, en lugar de Lambda. Sin embargo, creo que podría adoptarse para evaluar el error de Lambda utilizando la misma metodología. Se realiza un remuestreo de arranque del campo de datos mediante el remuestreo de las filas hasta que se produce un nuevo campo. La misma fila se puede volver a muestrear más de una vez, lo cual es un enfoque no paramétrico y no es necesario hacer suposiciones sobre la distribución de datos.

Bootstrap de valores Lambda

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)

Si bien esto puede ser más robusto que la regla general de North para calcular el error de los valores de Lambda, creo que ahora la cuestión de la importancia de EOF se reduce a diferentes opiniones sobre lo que esto significa. Para la regla general del Norte y los métodos bootstrap, la importancia parece estar más basada en si existe o no una superposición entre los valores de Lambda. Si es así, estos EOF pueden mezclarse en sus señales y no representar patrones "verdaderos". Por otro lado, estos dos EOF pueden describir una cantidad significativa de variación (en comparación con la descomposición de un campo aleatorio, por ejemplo, la Regla N). Entonces, si uno está interesado en filtrar el ruido (es decir, a través del truncamiento EOF), entonces la Regla N sería suficiente. Si uno está interesado en determinar patrones reales en un conjunto de datos, entonces los criterios más estrictos de superposición de Lambda pueden ser más sólidos.

Nuevamente, no soy un experto en estos temas, así que todavía espero que alguien con más experiencia pueda agregar a esta explicación.

Referencias

Beckers, Jean-Marie y M. Rixen. "Cálculos EOF y llenado de datos de conjuntos de datos oceanográficos incompletos". Revista de tecnología atmosférica y oceánica 20.12 (2003): 1839-1856.

Overland, J. y R. Preisendorfer, una prueba de significación para los componentes principales aplicados a una climatología ciclónica, lunes. Wea. Rev., 110, 1-4, 1982.