Flechas de variables subyacentes en PCA biplot en R

Considere la posibilidad de votar en la publicación de @ amoeba y @ttnphns . Gracias a ambos por su ayuda e ideas.

El siguiente se basa en el conjunto de datos del iris en R , y específicamente los primeros tres variables (columnas): Sepal.Length, Sepal.Width, Petal.Length.

Un biplot combina una gráfica de carga (vectores propios no estandarizados) - en concreto, las dos primeras cargas , y una gráfica de puntuación (puntos de datos rotados y dilatados trazados con respecto a los componentes principales). Utilizando el mismo conjunto de datos, @amoeba describe 9 posibles combinaciones de PCA biplot basadas en 3 posibles normalizaciones del gráfico de puntuación del primer y segundo componente principal, y 3 normalizaciones del gráfico de carga (flechas) de las variables iniciales. Para ver cómo R maneja estas posibles combinaciones, es interesante observar el biplot()método:

Primero el álgebra lineal lista para copiar y pegar:

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. Reproducción del diagrama de carga (flechas):

Aquí la interpretación geométrica en esta publicación de @ttnphns ayuda mucho. La notación del diagrama en la publicación se ha mantenido: representa la variable en el espacio temático . es la flecha correspondiente finalmente trazada; y las coordenadas y son el componente que carga una variable con respecto a y : $V$ Sepal L. $h'$ $a_1$ $a_2$ $V$ $\small \text{PC} 1$ $\small \text{PC} 2$

El componente de la variable Sepal L.con respecto a será: $\small\text{PC}1$

\begin{aligned} a_{1} & = h \cdot \cos (ϕ) \end{aligned}

$\begin{align} a_1 &= h\cdot\cos(\phi)\\[2ex] \end{align}$

que, si las puntuaciones con respecto a - llamémoslas - están estandarizadas para que su $\small\text{PC}1$ $\small\text{S}1$

$\Vert\text{S}1\Vert = \sqrt{\sum_1^n \text{scores}_1^2} = 1$ , la ecuación anterior es equivalente al producto de punto : $V\cdot \text{S}1$

\begin{aligned} a_{1} & = V \cdot S 1 \\ = ‖ V ‖ ‖ S 1 ‖ \cos (ϕ) \\ (1) & = h \times 1 \times \cdot \cos (ϕ) \end{aligned}

$\begin{align} a_1 &= V\cdot \text{S}1\\[2ex] &=\Vert V\Vert\,\Vert \text{S}1\Vert\, \cos(\phi)\\[2ex] &= h\times 1\times \cdot\cos(\phi)\tag{1} \end{align}$

Como , $\Vert V \Vert=\sqrt{\small{\sum x^2}}$

\sqrt{Var (V)} = \frac{\sqrt{\sum x^{2}}}{\sqrt{n - 1}} = \frac{‖ V ‖}{\sqrt{n - 1}} ⟹ ‖ V ‖ = h = \sqrt{var (V)} \sqrt{n - 1} .

$\sqrt{\small{\text{Var}(V)}}=\frac{\sqrt{\small{\sum x^2}}}{\sqrt{n-1}}=\frac{\Vert V \Vert}{\sqrt{n-1}} \implies \Vert V\Vert =h=\sqrt{\small{\text{var}(V)}} \sqrt {n-1}.$

Igualmente,

‖ S 1 ‖ = 1 = \sqrt{var(S 1)} \sqrt{n - 1} .

$\Vert\text{S}1\Vert=1=\sqrt{\small \text{var(S}1)}\sqrt {n-1}.$

Volviendo a la ecuación. , $(1)$

a_{1} = h \times 1 \times \cdot \cos (ϕ) = \sqrt{var (V)} \sqrt{var (S 1)} \cos (θ) (n - 1)

$a_1=h\times 1\times \cdot\cos(\phi)=\sqrt{\small{\text{var}(V)}}\,\sqrt{\small{\text{var}(\text{S}1)}}\, \cos(\theta) \;(n-1)$

$\cos(\phi)$ puede, por lo tanto, considerarse un coeficiente de correlación de Pearson , , con la advertencia de que no entiendo la arruga del factor . $r$ $n-1$

Duplicando y superponiendo en azul las flechas rojas 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)

Puntos de interés:

Las flechas se pueden reproducir como la correlación de las variables originales con los puntajes generados por los dos primeros componentes principales.
Alternativamente, esto se puede lograr como en el primer gráfico en la segunda fila, etiquetado en la publicación de @ amoeba: $\mathbf{ V*S}$

o en código 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)

o incluso aún ...

    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)

conectando con la explicación geométrica de las cargas por @ttnphns , o esta otra publicación informativa también por @ttnphns .

Hay un factor de escala sqrt(nrow(X) - 1), que sigue siendo un misterio.
$0.8$ tiene que ver con la creación de espacio para la etiqueta; vea este comentario aquí :

Además, ¡se debe decir que las flechas se trazan de manera que el centro de la etiqueta de texto esté donde debería estar! Luego, las flechas se multiplican por 0.80.8 antes de trazar, es decir, todas las flechas son más cortas de lo que deberían ser, presumiblemente para evitar la superposición con la etiqueta de texto (consulte el código para biplot.default). Me parece que esto es extremadamente confuso. - ameba Mar 19 '15 a las 10:06

2. Trazar la gráfica de `biplot()`puntajes (y flechas simultáneamente):

Los ejes se escalan a la unidad de suma de cuadrados, correspondiente a la primera gráfica de la primera fila en la publicación de @ ameeba , que se puede reproducir trazando la matriz de la descomposición de svd (más sobre esto más adelante) - " Columnas de : estos son componentes principales escalados a la unidad de suma de cuadrados " . $\mathbf U$ $\mathbf U$

Hay dos escalas diferentes en juego en los ejes horizontales inferior y superior en la construcción biplot:

Sin embargo, la escala relativa no es inmediatamente obvia, lo que requiere profundizar en las funciones y métodos:

biplot()traza los puntajes como columnas de en SVD, que son vectores unitarios ortogonales: $\mathbf 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

Mientras que la prcomp()función en R devuelve los puntajes escalados a sus valores propios:

> 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

Por lo tanto, podemos escalar la varianza a dividiendo por los valores propios: $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

Pero como queremos que la suma de los cuadrados sea , necesitaremos dividir por porque: $1$ $\sqrt{n-1}$

var (scr_var_one) = 1 = \frac{\sum_{1}^{n} scr_var_one}{n - 1}

$\small \text{var}(\text{scr_var_one})= 1 =\frac{\sum_1^n \text{scr_var_one}}{n -1}$

> 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

Cabe destacar que el uso del factor de escala , se cambia más adelante a cuando la definición de la explicación parece estar en el hecho de que $\sqrt{n-1}$ $\sqrt{n}$ lan

prcompusa : "A diferencia de princomp, las variaciones se calculan con el divisor habitual ". $n-1$ $n - 1$

Después de despojarlos de todas las ifdeclaraciones y otras pelusas de limpieza, biplot()procede de la siguiente manera:

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)

que, como se esperaba, reproduce (imagen derecha abajo) la biplot()salida como se llama directamente con biplot(PCA)(gráfico izquierdo abajo) en todas sus deficiencias estéticas intactas:

Puntos de interés:

Las flechas se trazan en una escala relacionada con la relación máxima entre el vector propio escalado de cada uno de los dos componentes principales y sus respectivos puntajes escalados (el ratio). Como comenta @amoeba:

el diagrama de dispersión y el "diagrama de flecha" se escalan de manera que la coordenada de flecha x o y más grande (en valor absoluto) de las flechas era exactamente igual a la coordenada x o y más grande (en valor absoluto) de los puntos de datos dispersos

Como se anticipó anteriormente, los puntos se pueden trazar directamente como las puntuaciones en la matriz de la SVD: $\mathbf U$

Antoni Parellada
fuente

+1, buen estudio. Agregué una etiqueta Ra su pregunta porque el asunto confuso (es decir, el coeficiente de escala) resultó ser parcialmente específico de R. En general, pudo ver que el biplot PCA es un diagrama de dispersión superpuesto de puntajes de componentes (coordenadas de fila) y coeficientes de dirección de componentes (coordenadas de columna), y dado que se pueden aplicar varias cantidades de estandarizaciones por "inercia" (varianza) a cada del también, por lo que pueden surgir varias miradas del biplot. Para agregar: generalmente (más sentido), las cargas se muestran como las coordenadas de la columna (flechas).

ttnphns

(cont.) Vea mi resumen de biplot que explica, en diferentes palabras, lo que ha mostrado en su buena respuesta.

ttnphns

¡+ 1 gracias por escribir tutoriales y con código reproducible para nosotros!

Haitao Du

Antony, ¿dibujaste (como a mano) o trazaste (alimentaste con datos) tu foto? ¿Qué software usaste? Se ve bien.

ttnphns

@ttnphns Gracias! Aquí está el enlace . Me preguntaba si podría mejorarlo y trazar las cargas y las PC de una manera mejor y más didáctica. Siéntase libre de cambiar (es un programa maravillosamente fácil de usar), y si lo hace, por favor comparta.

Antoni Parellada

Flechas de variables subyacentes en PCA biplot en R

Respuestas:

1. Reproducción del diagrama de carga (flechas):

2. Trazar la gráfica de biplot()puntajes (y flechas simultáneamente):

2. Trazar la gráfica de `biplot()`puntajes (y flechas simultáneamente):