Comment déterminer les composants principaux significatifs en utilisant l’amorçage ou l’approche de Monte Carlo?

40

Je suis intéressé par la détermination du nombre de régularités significatives issues d'une analyse en composantes principales (ACP) ou d'une fonction empirique orthogonale (EOF). Je suis particulièrement intéressé par l'application de cette méthode aux données climatiques. Le champ de données est une matrice MxN, M étant la dimension temporelle (par exemple, les jours) et N étant la dimension spatiale (par exemple, les emplacements lon / lat). J'ai lu une méthode de bootstrap possible pour déterminer les PC significatifs, mais je n'ai pas pu trouver une description plus détaillée. Jusqu'à présent, j'ai appliqué la règle de base de North (North et al ., 1982) pour déterminer ce seuil, mais je me demandais si une méthode plus robuste était disponible.

Par exemple:

###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="")

entrez la description de l'image ici

Et voici la méthode que j’utilise pour déterminer l’importance du PC. Fondamentalement, la règle empirique est que la différence entre les Lambdas voisins doit être supérieure à l'erreur qui leur est associée.

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

entrez la description de l'image ici

J'ai trouvé utiles les sections de chapitres de Björnsson et Venegas ( 1997 ) sur les tests de signification. Elles renvoient à trois catégories de tests, dont le type de variance dominant correspond probablement à ce que j'espère utiliser. Le terme fait référence à un type d’approche de Monte Carlo consistant à mélanger la dimension temporelle et à recalculer les Lambdas sur de nombreuses permutations. von Storch et Zweiers (1999) font également référence à un test qui compare le spectre Lambda à un spectre de référence "bruit". Dans les deux cas, je suis un peu incertain de la façon dont cela pourrait être fait et de la manière dont le test de signification est effectué compte tenu des intervalles de confiance identifiés par les permutations.

Merci de votre aide.

Références: Björnsson, H. et Venegas, SA (1997). "Manuel pour l'analyse de données climatiques effectuée par EOF et SVD", Université McGill, Rapport du CCGCR n ° 97-1, Montréal, Québec, 52 p. http://andvari.vedur.is/%7Efolk/halldor/PICKUP/eof.pdf

GR Nord, TL Bell, RF Cahalan et FJ Moeng. (1982). Erreurs d'échantillonnage dans l'estimation des fonctions orthogonales empiriques. Lun. Wea. Rev., 110: 699–706.

von Storch, H. Zwiers, FW (1999). Analyse statistique dans la recherche climatique. La presse de l'Universite de Cambridge.

Marc dans la boîte
la source
Quelle est votre référence sur l'approche bootstrap?
Michael R. Chernick
4
Un bootstrap ne fonctionnera pas ici. Cela ne fonctionnera pas avec les ensembles de données dans lesquels presque toutes les observations sont corrélées avec à peu près toutes les autres observations; il a besoin d'indépendance, ou du moins d'indépendance approximative (conditions de mélange dans une série chronologique, par exemple) pour produire des répliques justifiables des données. Bien sûr, il existe des systèmes d'amorçage spéciaux, tels que l'amorçage sauvage, qui peuvent contourner ces problèmes. Mais je ne parierai pas beaucoup là-dessus. Et vous devez vraiment consulter les statistiques multivariées et les suivre afin de ne pas obtenir un autre bâton de hockey indéfendable en guise de réponse.
StasK
2
@Marc dans la zone, vous faites peut-être référence à différents bootstrap de bloc utilisés pour les séries chronologiques, MBB (bootstrap de bloc mobile), CBB (bootstrap de bloc circulaire) ou SBB (bootstrap de bloc fixe) qui utilisent des blocs de temps des données pour estimer le modèle. paramètres.
Michael R. Chernick
3
@StasK Je ne sais pas pourquoi vous pensez avoir besoin de conditions de mélange pour appliquer bootstrap aux séries chronologiques. Les méthodes basées sur un modèle exigent simplement que vous ajustiez une structure de série chronologique et que vous puissiez ensuite amorcer des résidus. Vous pouvez donc avoir une série chronologique avec les tendances et les composantes saisonnières tout en effectuant un bootstrap basé sur un modèle.
Michael R. Chernick
2
Je n'ai pas accès au texte intégral, mais vous pouvez essayer d'y jeter un coup d'œil: "Hamid Babamoradi, Frans van den Berg, Åsmund Rinnan, Limites de confiance basées sur Bootstrap dans l'analyse en composantes principales - Étude de cas, Chimiométrie et systèmes de laboratoire intelligents, Volume 120, 15 janvier 2013, pages 97-105, ISSN 0169-7439, 10.1016 / j.chemolab.2012.10.007. ( Sciencedirect.com/science/article/pii/S0169743912002171 ) Mots-clés: Bootstrap; PCA; Limites de confiance; BC < sous "a </ sub>; Incertitude"
tomasz74

Réponses:

19

Je vais essayer de faire avancer un peu le dialogue ici même si telle est ma question. Cela fait six mois que j'ai posé cette question et, malheureusement, aucune réponse complète n'a été donnée. Je vais essayer de résumer ce que j'ai rassemblé jusqu'à présent et de voir si quelqu'un peut élaborer sur les problèmes en suspens. Veuillez excuser la longue réponse, mais je ne vois pas d'autre moyen ...

Tout d’abord, je vais démontrer plusieurs approches en utilisant un ensemble de données synthétiques éventuellement meilleur. Il provient d'un article de Beckers et Rixon ( 2003 ) illustrant l'utilisation d'un algorithme pour effectuer une analyse EOF sur des données gappy. J'ai reproduit l'algorithme en R si quelqu'un est intéressé ( lien ).

Ensemble de données synthétiques:

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

entrez la description de l'image ici

Ainsi, le champ de données exact Xtest composé de 9 signaux et j’y ai ajouté du bruit pour créer le champ observé Xp, qui sera utilisé dans les exemples ci-dessous. Les EOF sont déterminés comme tels:

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

Suivant l'exemple que j'ai utilisé dans mon exemple initial, je déterminerai les EOF "significatives" selon la règle empirique de North.

La règle du Nord

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)

entrez la description de l'image ici

Etant donné que les valeurs Lambda de 2: 4 sont très proches les unes des autres en amplitude, elles sont considérées comme non significatives par la règle générale - en d'autres termes, leurs modèles EOF respectifs peuvent se chevaucher et se mélanger compte tenu de leurs amplitudes similaires. C'est malheureux étant donné que nous savons que 9 signaux existent réellement sur le terrain.

Une approche plus subjective consiste à afficher les valeurs Lambda transformées en journal ("tracé de balayage"), puis à ajuster une régression aux valeurs de fin. On peut alors déterminer visuellement à quel niveau les valeurs lambda se situent au-dessus de cette ligne:

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)

entrez la description de l'image ici

Ainsi, les 5 premiers EOF se situent au-dessus de cette ligne. J'ai essayé cet exemple lorsqu'aucun Xpbruit supplémentaire n'a été ajouté et que les résultats révèlent les 9 signaux d'origine. Ainsi, l’insignifiance des EOF 6: 9 est due au fait que leur amplitude est inférieure au bruit sur le terrain.

Une méthode plus objective est le critère "Rule N" de Overland et Preisendorfer (1982). Il y a une implémentation dans le wqpaquet, que je montre ci-dessous.

Règle 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)

entrez la description de l'image ici

La règle N identifiait 4 EOF significatifs. Personnellement, je dois mieux comprendre cette méthode. Pourquoi est-il possible d'évaluer le niveau d'erreur en fonction d'un champ aléatoire qui n'utilise pas la même distribution que celle-ci Xp? Une variante de cette méthode consisterait à rééchantillonner les données Xpafin que chaque colonne soit redistribuée de manière aléatoire. De cette manière, nous nous assurons que la variance totale du champ aléatoire est la même que Xp. En rééchantillonnant plusieurs fois, nous sommes alors en mesure de calculer une erreur de base de la décomposition.

Monte Carlo avec champ aléatoire (comparaison de modèles nuls)

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)

entrez la description de l'image ici

Encore une fois, 4 EOF sont au-dessus des distributions pour les champs aléatoires. Mon souci avec cette approche, et celle de la règle N, est que celles-ci ne traitent pas vraiment les intervalles de confiance des valeurs de Lambda; Par exemple, une première valeur lambda élevée se traduira automatiquement par une quantité de variance inférieure à expliquer par les valeurs de fin. Ainsi, la valeur Lambda calculée à partir de champs aléatoires aura toujours une pente moins raide et peut entraîner une sélection de trop peu d'EOF significatifs. [NOTE: La eofNum()fonction suppose que les EOF sont calculés à partir d'une matrice de corrélation. Ce nombre peut être différent si vous utilisez par exemple une matrice de covariance (données centrées mais non mises à l'échelle).]

Enfin, @ tomasz74 a mentionné le document de Babamoradi et al. (2013), que j'ai examiné brièvement. C'est très intéressant, mais semble être plus axé sur le calcul des CI des charges et des coefficients EOF, plutôt que sur Lambda. Néanmoins, j'estime qu'il pourrait être adopté pour évaluer l'erreur Lambda en utilisant la même méthodologie. Un rééchantillonnage bootstrap du champ de données est effectué en rééchantillonnant les lignes jusqu'à la production d'un nouveau champ. La même ligne peut être rééchantillonnée plus d'une fois, ce qui constitue une approche non paramétrique et il n'est pas nécessaire de faire d'hypothèses sur la distribution des données.

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)

entrez la description de l'image ici

Bien que cela puisse sembler plus robuste que la règle empirique de North pour calculer l'erreur des valeurs Lambda, je crois maintenant que la question de la signification de l'EOF revient à donner des opinions différentes sur ce que cela signifie. Pour les méthodes empiriques et de bootstrap du Nord, la signification semble être davantage basée sur le fait de savoir s'il existe ou non un chevauchement des valeurs Lambda. Si tel est le cas, ces EOF peuvent être mélangés dans leurs signaux et ne pas représenter de "véritables" modèles. D'autre part, ces deux EOF peuvent décrire une quantité de variance significative (par rapport à la décomposition d'un champ aléatoire - par exemple la règle N). Donc, si on est intéressé par le filtrage du bruit (c'est-à-dire via la troncature EOF), alors la règle N sera suffisante. Si l'on souhaite déterminer des modèles réels dans un ensemble de données, les critères plus rigoureux de chevauchement de Lambda peuvent être plus robustes.

Encore une fois, je ne suis pas un expert sur ces questions, j'espère donc que quelqu'un avec plus d'expérience peut ajouter à cette explication.

Les références:

Beckers, Jean-Marie et M. Rixen. "Calculs EOF et remplissage de données à partir d'ensembles de données océanographiques incomplètes." Journal de la technologie atmosphérique et océanique 20.12 (2003): 1839-1856.

Overland, J., et R. Preisendorfer, Un test de signification pour les principaux composants appliqués à une climatologie cyclonique, Mon. Wea. Rev., 110, 1-4, 1982.

Marc dans la boîte
la source