Flèches des variables sous-jacentes dans le biplot PCA dans R

11

Au risque de rendre la question spécifique au logiciel, et avec l'excuse de son ubiquité et de ses idiosyncrasies, je veux poser des questions sur la fonction biplot()dans R, et, plus précisément, sur le calcul et le tracé de sa valeur par défaut, des flèches rouges superposées, correspondant aux variables sous-jacentes.


[Pour donner un sens à certains des commentaires, les complots initialement affichés avaient un problème de plus en plus rare et sont maintenant effacés.]

Antoni Parellada
la source
Je ne peux pas comprendre comment vous avez obtenu vos flèches vertes. Ils ne sont pas corrects. Le fait que la longueur verte soit d'env. deux fois plus longtemps que la largeur verte permet de suspecter que vous traçiez les vecteurs relatifs aux variables non standardisées. Cela ne peut pas se produire sur un biplot de PCA basé sur des corrélations.
ttnphns
Les flèches rouges semblent correctes. Voir: ils sont de même longueur et sont symétriques par rapport à PC2. C'est la seule position possible lorsque vous effectuez l'ACP avec seulement 2 variables et basé sur des corrélations (c'est-à-dire des variables standardisées). Dans l'ACP basée sur les corrélations, les chargements (coordonnés des flèches) sont les corrélations entre les PC et les variables. Dans votre exemple, les chargements sont (vars par PC): .74752, .66424; -.74752, .66424.
ttnphns
@ttnphns Oui, les flèches rouges sont ce que j'essaie de reproduire (elles sont correctes), et elles sont tracées en R avec l' biplot(name_of_the_PCA)appel, ce qui est le cas dans ce cas biplot(PCA). J'avais centré et mis à l'échelle les données.
Antoni Parellada
Alors, quelle est votre question? Comment calculer les coordonnées des flèches rouges? Ils devraient être les chargements PCA . Parfois, les vecteurs propres sont tracés (votre commande R l'a probablement fait ??), cependant, la manière consensuelle et significative est de tracer les chargements .
ttnphns
@ttnphns Tracer les vecteurs propres (je suppose que c'est la même chose que les chargements) me donne la bonne orientation (merci), mais pas la même ampleur que les flèches rouges (je colle l'image dans l'OP).
Antoni Parellada

Réponses:

19

Pensez à voter pour @ amoeba et @ttnphns . Merci à vous deux pour votre aide et vos idées.


Ce qui suit repose sur l' ensemble de données d' iris dans R , et en particulier les trois premières variables (colonnes) Sepal.Length, Sepal.Width, Petal.Length.

Un biplot combine un tracé de chargement (vecteurs propres non normalisés) - dans le béton, les deux premiers chargements , et un tracé de score (points de données tournés et dilatés tracés par rapport aux principaux composants). En utilisant le même ensemble de données, @amoeba décrit 9 combinaisons possibles de biplot PCA basées sur 3 normalisations possibles du graphique de score des première et deuxième composantes principales, et 3 normalisations du graphique de chargement (flèches) des variables initiales. Pour voir comment R gère ces combinaisons possibles, il est intéressant de regarder la biplot()méthode:


D'abord l'algèbre linéaire prête à copier et coller:

X = as.matrix(iris[,1:3])             # Three first variables of Iris dataset
CEN = scale(X, center = T, scale = T) # Centering and scaling the data
PCA = prcomp(CEN)

# EIGENVECTORS:
(evecs.ei = eigen(cor(CEN))$vectors)       # Using eigen() method
(evecs.svd = svd(CEN)$v)                   # PCA with SVD...
(evecs = prcomp(CEN)$rotation)             # Confirming with prcomp()

# EIGENVALUES:
(evals.ei = eigen(cor(CEN))$values)        # Using the eigen() method
(evals.svd = svd(CEN)$d^2/(nrow(X) - 1))   # and SVD: sing.values^2/n - 1
(evals = prcomp(CEN)$sdev^2)               # with prcomp() (needs squaring)

# SCORES:
scr.svd = svd(CEN)$u %*% diag(svd(CEN)$d)  # with SVD
scr = prcomp(CEN)$x                        # with prcomp()
scr.mm = CEN %*% prcomp(CEN)$rotation      # "Manually" [data] [eigvecs]

# LOADINGS:

loaded = evecs %*% diag(prcomp(CEN)$sdev)  # [E-vectors] [sqrt(E-values)]

1. Reproduction du tracé de chargement (flèches):

Ici, l'interprétation géométrique de ce message de @ttnphns aide beaucoup. La notation du diagramme dans le message a été conservée: représente la variable dans l' espace sujet . est la flèche correspondante finalement tracée; et les coordonnées et sont le composant charge une variable par rapport à et :h a 1 a 2 V PC 1 PC 2VSepal L.ha1a2VPC1PC2


entrez la description de l'image ici


La composante de la variable Sepal L.par rapport à sera alors:PC1

a1=hcos(ϕ)

qui, si les scores par rapport à - appelons-les - sont standardisés afin que leurS 1PC1S1

S1=1nscores12=1 , l'équation ci-dessus est l'équivalent du produit scalaire :VS1

a1=VS1=VS1cos(ϕ)(1)=h×1×cos(ϕ)

Depuis ,V=x2

Var(V)=x2n1=Vn1V=h=var(V)n1.

Également,

S1=1=var(S1)n1.

Revenons à l'Eq. ,(1)

a1=h×1×cos(ϕ)=var(V)var(S1)cos(θ)(n1)

cos(ϕ) peut donc être considéré comme un coefficient de corrélation de Pearson , , avec la réserve que je ne comprends pas la ride du facteur .rn1

Dupliquer et chevaucher en bleu les flèches rouges de biplot()

par(mfrow = c(1,2)); par(mar=c(1.2,1.2,1.2,1.2))

biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
# R biplot with overlapping (reproduced) arrows in blue completely covering red arrows:
biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
arrows(0, 0,
       cor(X[,1], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,1], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
arrows(0, 0,
       cor(X[,2], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,2], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
arrows(0, 0,
       cor(X[,3], scr[,1]) * 0.8 * sqrt(nrow(X) - 1), 
       cor(X[,3], scr[,2]) * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)

entrez la description de l'image ici

Points d'interêts:

  • Les flèches peuvent être reproduites comme la corrélation des variables d'origine avec les scores générés par les deux premières composantes principales.
  • Alternativement, cela peut être réalisé comme dans le premier graphique de la deuxième ligne, étiqueté dans le post de @ amoeba:VS

entrez la description de l'image ici

ou en code R:

    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
    # R biplot with overlapping arrows in blue completely covering red arrows:
    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[1,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[1,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[2,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[2,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (svd(CEN)$v %*% diag(svd(CEN)$d))[3,1] * 0.8, 
       (svd(CEN)$v %*% diag(svd(CEN)$d))[3,2] * 0.8, 
       lwd = 1, angle = 30, length = 0.1, col = 4)

ou même encore ...

    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) # R biplot
    # R biplot with overlapping (reproduced) arrows in blue completely covering red arrows:
    biplot(PCA, cex = 0.6, cex.axis = .6, ann = F, tck=-0.01) 
    arrows(0, 0,
       (loaded)[1,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[1,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (loaded)[2,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[2,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)
    arrows(0, 0,
       (loaded)[3,1] * 0.8 * sqrt(nrow(X) - 1), 
       (loaded)[3,2] * 0.8 * sqrt(nrow(X) - 1), 
       lwd = 1, angle = 30, length = 0.1, col = 4)

se connecter avec l' explication géométrique des chargements par @ttnphns , ou cet autre poste informatif également par @ttnphns .

  • Il y a un facteur d'échelle:, sqrt(nrow(X) - 1)qui reste un peu mystérieux.

  • 0.8 a à voir avec la création d'espace pour l'étiquette - voir ce commentaire ici :

De plus, il faut dire que les flèches sont tracées de telle sorte que le centre de l'étiquette de texte soit là où il devrait être! Les flèches sont ensuite multipliées par 0,80,8 avant le traçage, c'est-à-dire que toutes les flèches sont plus courtes que ce qu'elles devraient être, probablement pour éviter le chevauchement avec l'étiquette de texte (voir le code pour biplot.default). Je trouve que c'est extrêmement déroutant. - amoeba 19 mars 15 à 10:06


2. Tracer le biplot()tracé des scores (et les flèches simultanément):

Les axes sont mis à l'échelle en une somme unitaire de carrés, correspondant au premier tracé de la première ligne sur le post de @ amoeba , qui peut être reproduit en traçant la matrice de la décomposition svd (plus de détails plus loin) - " Colonnes de : ce sont les principaux composants mis à l'échelle de la somme unitaire des carrés. "UU

entrez la description de l'image ici

Il y a deux échelles différentes en jeu sur les axes horizontaux inférieur et supérieur dans la construction du biplot:

entrez la description de l'image ici

Cependant, l'échelle relative n'est pas immédiatement évidente, ce qui nécessite de se plonger dans les fonctions et les méthodes:

biplot()trace les scores sous forme de colonnes de dans SVD, qui sont des vecteurs unitaires orthogonaux:U

> scr.svd = svd(CEN)$u %*% diag(svd(CEN)$d) 
> U = svd(CEN)$u
> apply(U, 2, function(x) sum(x^2))
[1] 1 1 1

Alors que la prcomp()fonction dans R renvoie les scores mis à l'échelle à leurs valeurs propres:

> apply(scr, 2, function(x) var(x))         # pr.comp() scores scaled to evals
       PC1        PC2        PC3 
2.02142986 0.90743458 0.07113557 
> evals                                     #... here is the proof:
[1] 2.02142986 0.90743458 0.07113557

Par conséquent, nous pouvons mettre la variance à en divisant par les valeurs propres:1

> scr_var_one = scr/sqrt(evals)[col(scr)]  # to scale to var = 1
> apply(scr_var_one, 2, function(x) var(x)) # proved!
[1] 1 1 1

Mais comme nous voulons que la somme des carrés soit , nous devons diviser par car:1n1

var(scr_var_one)=1=1nscr_var_onen1
> scr_sum_sqrs_one = scr_var_one / sqrt(nrow(scr) - 1) # We / by sqrt n - 1.
> apply(scr_sum_sqrs_one, 2, function(x) sum(x^2))     #... proving it...
PC1 PC2 PC3 
  1   1   1

Il convient de noter que l'utilisation du facteur d'échelle , est modifiée ultérieurement en lorsque la définition de l'explication semble résider dans le fait quen1nlan

prcomputilise : "Contrairement à princomp, les variances sont calculées avec le diviseur habituel ".n - 1n1n1


Après les avoir dépouillés de toutes les ifdéclarations et autres peluches de ménage, biplot()procède comme suit:

X   = as.matrix(iris[,1:3])                    # The original dataset
CEN = scale(X, center = T, scale = T)          # Centered and scaled
PCA = prcomp(CEN)                              # PCA analysis

par(mfrow = c(1,2))                            # Splitting the plot in 2.
biplot(PCA)                                    # In-built biplot() R func.

# Following getAnywhere(biplot.prcomp):

choices = 1:2                                  # Selecting first two PC's
scale = 1                                      # Default
scores= PCA$x                                  # The scores
lam = PCA$sdev[choices]                        # Sqrt e-vals (lambda) 2 PC's
n = nrow(scores)                               # no. rows scores
lam = lam * sqrt(n)                            # See below.

# at this point the following is called...
# biplot.default(t(t(scores[,choices])      /  lam), 
#                t(t(x$rotation[,choices]) *   lam))

# Following from now on getAnywhere(biplot.default):

x = t(t(scores[,choices])       / lam)         # scaled scores
# "Scores that you get out of prcomp are scaled to have variance equal to      
#  the eigenvalue. So dividing by the sq root of the eigenvalue (lam in 
#  biplot) will scale them to unit variance. But if you want unit sum of 
#  squares, instead of unit variance, you need to scale by sqrt(n)" (see comments).
# > colSums(x^2)
# PC1       PC2 
# 0.9933333 0.9933333    # It turns out that the it's scaled to sqrt(n/(n-1)), 
# ...rather than 1 (?) - 0.9933333=149/150

y = t(t(PCA$rotation[,choices]) * lam)         # scaled eigenvecs (loadings)


n = nrow(x)                                    # Same as dataset (150)
p = nrow(y)                                    # Three var -> 3 rows

# Names for the plotting:

xlabs = 1L:n
xlabs = as.character(xlabs)                    # no. from 1 to 150 
dimnames(x) = list(xlabs, dimnames(x)[[2L]])   # no's and PC1 / PC2

ylabs = dimnames(y)[[1L]]                      # Iris species
ylabs = as.character(ylabs)
dimnames(y) <- list(ylabs, dimnames(y)[[2L]])  # Species and PC1/PC2

# Function to get the range:
unsigned.range = function(x) c(-abs(min(x, na.rm = TRUE)), 
                                abs(max(x, na.rm = TRUE)))
rangx1 = unsigned.range(x[, 1L])               # Range first col x
# -0.1418269  0.1731236
rangx2 = unsigned.range(x[, 2L])               # Range second col x
# -0.2330564  0.2255037
rangy1 = unsigned.range(y[, 1L])               # Range 1st scaled evec
# -6.288626   11.986589
rangy2 = unsigned.range(y[, 2L])               # Range 2nd scaled evec
# -10.4776155   0.8761695

(xlim = ylim = rangx1 = rangx2 = range(rangx1, rangx2))
# range(rangx1, rangx2) = -0.2330564  0.2255037

# And the critical value is the maximum of the ratios of ranges of 
# scaled e-vectors / scaled scores:

(ratio = max(rangy1/rangx1, rangy2/rangx2)) 
# rangy1/rangx1   =   26.98328    53.15472
# rangy2/rangx2   =   44.957418   3.885388
# ratio           =   53.15472

par(pty = "s")                                 # Calling a square plot

# Plotting a box with x and y limits -0.2330564  0.2255037
# for the scaled scores:

plot(x, type = "n", xlim = xlim, ylim = ylim)  # No points
# Filling in the points as no's and the PC1 and PC2 labels:
text(x, xlabs) 
par(new = TRUE)                                # Avoids plotting what follows separately

# Setting now x and y limits for the arrows:

(xlim = xlim * ratio)  # We multiply the original limits x ratio
# -16.13617  15.61324
(ylim = ylim * ratio)  # ... for both the x and y axis
# -16.13617  15.61324

# The following doesn't change the plot intially...
plot(y, axes = FALSE, type = "n", 
     xlim = xlim, 
     ylim = ylim, xlab = "", ylab = "")

# ... but it does now by plotting the ticks and new limits...
# ... along the top margin (3) and the right margin (4)
axis(3); axis(4)
text(y, labels = ylabs, col = 2)  # This just prints the species

arrow.len = 0.1                   # Length of the arrows about to plot.

# The scaled e-vecs are further reduced to 80% of their value
arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, 
       length = arrow.len, col = 2)

qui, comme prévu, reproduit (image de droite ci-dessous) la biplot()sortie appelée directement avec biplot(PCA)(graphique de gauche ci-dessous) dans toutes ses imperfections esthétiques intactes:

entrez la description de l'image ici

Points d'interêts:

  • Les flèches sont tracées à une échelle liée au rapport maximal entre le vecteur propre à l'échelle de chacune des deux composantes principales et leurs scores à l'échelle respectifs (le ratio). AS @amoeba commente:

le nuage de points et le "diagramme de flèches" sont mis à l'échelle de telle sorte que la plus grande (en valeur absolue) des coordonnées de flèche x ou y des flèches était exactement égale à la plus grande (en valeur absolue) des coordonnées x ou y des points de données dispersés

  • Comme prévu ci-dessus, les points peuvent être directement tracés comme les scores dans la matrice du SVD:U

Antoni Parellada
la source
1
+1, bonne étude. J'ai ajouté une balise Rà votre question parce que la matière déroutante (à savoir, le coefficient de mise à l'échelle) s'est avérée être en partie spécifique à R. En général, vous avez pu constater par vous-même que le biplot PCA est un nuage de points de superposition des scores des composants (coordonnées des lignes) et des coefficients de direction des composants (coordonnées des colonnes), et puisque diverses quantités de standardisations par "inertie" (variance) peuvent être appliquées à chaque également, de sorte que différents aspects du biplot peuvent apparaître. Pour ajouter: le plus souvent (plus de sens), les chargements sont affichés sous forme de coordonnées de colonne (flèches).
ttnphns
1
(suite) Voir mon aperçu de biplot qui explique, en d'autres termes, ce que vous avez montré dans votre belle réponse.
ttnphns
2
+ 1 merci d'avoir écrit des tutoriels et avec du code reproductible pour nous!
Haitao Du
Antony, as-tu dessiné (comme à la main) ou as-tu tracé (alimenté en données) ta photo? Quel logiciel avez-vous utilisé? Ça à l'air bien.
ttnphns
@ttnphns Merci! Voici le lien vers celui-ci . Je me demandais si vous pouviez l'améliorer et tracer les chargements et les PC de manière plus didactique. N'hésitez pas à changer (c'est un programme merveilleusement convivial), et si vous le faites, veuillez le partager.
Antoni Parellada