Setas de variáveis ​​subjacentes no PCA biplot em R

11

Correndo o risco de tornar a pergunta específica do software e com a desculpa de sua onipresença e idiossincrasias, quero perguntar sobre a função biplot()em R e, mais especificamente, sobre o cálculo e a plotagem de suas setas vermelhas sobrepostas padrão, correspondentes às variáveis ​​subjacentes.


[Para entender alguns dos comentários, as parcelas postadas inicialmente tinham um problema de interesse escasso e agora são apagadas.]

Antoni Parellada
fonte
Não consigo descobrir como você conseguiu suas flechas verdes. Eles não estão corretos. O fato de o comprimento verde ser aprox. duas vezes mais que o verde s.width permite suspeitar que você estava plotando os vetores pertencentes às variáveis ​​não padronizadas. Isso não pode acontecer em um biplot de PCA com base em correlações.
ttnphns
Setas vermelhas parecem corretas. Veja: eles têm o mesmo comprimento e são simétricos em relação ao PC2. Essa é a única posição possível quando você faz PCA com apenas 2 variáveis ​​e com base em correlações (isto é, variáveis ​​padronizadas). No PCA baseado em correlações, as cargas (coordenadas pelas setas) são as correlações entre os PCs e as variáveis. No seu exemplo, as cargas são (VARs por PCs): .74752, .66424; -.74752, .66424.
ttnphns
@ttnphns Sim, as setas vermelhas são o que estou tentando reproduzir (elas estão corretas) e são plotadas em R com a biplot(name_of_the_PCA)chamada, que neste caso é biplot(PCA). Eu centralizei e dimensionei os dados.
Antoni Parellada
Então, qual é a sua pergunta? Como calcular as coordenadas das setas vermelhas? Eles devem ser os carregamentos de PCA . Às vezes, os vetores próprios são plotados (seu comando R provavelmente fez isso ??), no entanto, a maneira consensual e significativa é plotar cargas .
ttnphns
@ttnphns A plotagem dos autovetores (presumo que seja a mesma das cargas) me dá a orientação correta (obrigado), mas não a mesma magnitude que as setas vermelhas (estou colando a imagem no OP).
Antoni Parellada 29/04

Respostas:

19

Considere votar novamente na publicação de @ amoeba e @ttnphns . Obrigado por sua ajuda e idéias.


O que se segue baseia-se no conjunto de dados da íris em R , e especificamente as primeiras três variáveis (colunas): Sepal.Length, Sepal.Width, Petal.Length.

Um biplot combina uma plotagem de carga (autovetores não padronizados) - no concreto, as duas primeiras cargas e uma plotagem de pontuação (pontos de dados rotacionados e dilatados plotados em relação aos componentes principais). Utilizando o mesmo conjunto de dados, o @amoeba descreve 9 combinações possíveis de biplot de PCA com base em 3 possíveis normalizações do gráfico de pontuação do primeiro e segundo componentes principais e 3 normalizações do gráfico de carregamento (setas) das variáveis ​​iniciais. Para ver como R lida com essas combinações possíveis, é interessante observar o biplot()método:


Primeiro a álgebra linear pronta para copiar e colar:

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. Reprodução do gráfico de carregamento (setas):

Aqui a interpretação geométrica neste post por @ttnphns ajuda muito. A notação do diagrama no post foi mantida: representa a variável no espaço de assunto . é a seta correspondente finalmente plotada; e as coordenadas e são o componente carrega uma variável em relação a e :h ' a 1 a 2 V PC 1 PC 2VSepal L.ha1a2VPC1PC2


insira a descrição da imagem aqui


O componente da variável Sepal L.em relação a será então:PC1

a1=hcos(ϕ)

que, se as pontuações em relação a - vamos chamá-las de - são padronizadas para que suasS 1PC1S1

VS1S1=1nscores12=1 , a equação acima é equivalente ao produto escalar :VS1

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

Como ,V=x2

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

Da mesma forma,

S1=1=var(S1)n1.

Voltando à Eq. ,(1)

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

r n - 1cos(ϕ) pode, portanto, ser considerado um coeficiente de correlação de Pearson , , com a ressalva de que eu não entendo as rugas do fator .rn1

Duplicando e sobrepondo em azul as setas vermelhas 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)

insira a descrição da imagem aqui

Pontos de interesse:

  • As setas podem ser reproduzidas como a correlação das variáveis ​​originais com as pontuações geradas pelos dois primeiros componentes principais.
  • Como alternativa, isso pode ser alcançado como no primeiro gráfico da segunda linha, rotulado no post de @ amoeba:VS

insira a descrição da imagem aqui

ou no 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)

ou ainda ...

    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-se com a explicação geométrica das cargas por @ttnphns ou com este outro post informativo também por @ttnphns .

  • Há um fator de escala:, sqrt(nrow(X) - 1)que permanece um pouco misterioso.

  • 0.8 tem a ver com a criação de espaço para o rótulo - veja este comentário aqui :

Além disso, deve-se dizer que as setas são plotadas de forma que o centro do rótulo do texto esteja onde deveria estar! As setas são então multiplicadas por 0.80.8 antes da plotagem, ou seja, todas as setas são mais curtas do que deveriam, presumivelmente para evitar sobreposição com o rótulo do texto (consulte o código para biplot.default). Acho isso extremamente confuso. # Ameba 19/03/2015 às 10:06


2. Traçando o gráfico de biplot()pontuação (e setas simultaneamente):

Os eixos são dimensionados para soma unitária dos quadrados, correspondendo ao primeiro gráfico da primeira linha no post da @ amoeba , que pode ser reproduzido plotando a matriz da decomposição do svd (mais sobre isso mais adiante) - " Colunas de : esses são os principais componentes dimensionados para a soma unitária dos quadrados " .UUU

insira a descrição da imagem aqui

Existem duas escalas diferentes em jogo nos eixos horizontal inferior e superior na construção do biplot:

insira a descrição da imagem aqui

No entanto, a escala relativa não é imediatamente óbvia, exigindo a investigação das funções e métodos:

biplot()plota pontuações como colunas de em SVD, que são vetores de unidade ortogonais: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

Enquanto a prcomp()função em R retorna as pontuações dimensionadas para seus autovalores:

> 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

Portanto, podemos escalar a variação para dividindo pelos valores próprios: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

Mas como queremos que a soma dos quadrados seja , precisamos dividir por porque: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

Observe que o uso do fator de escala é alterado posteriormente para quando a definição da explicação parece residir no fato de quen1nlan

prcompusa : "Diferentemente do princomp, as variações são calculadas com o divisor usual ".n - 1n1n1


Depois de retirar todas as ifinstruções e outros cotões de limpeza da casa, biplot()proceda da seguinte maneira:

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 esperado, reproduz (imagem à direita abaixo) a biplot()saída conforme chamado diretamente com biplot(PCA)(plot esquerdo abaixo) em todas as suas deficiências estéticas intocadas:

insira a descrição da imagem aqui

Pontos de interesse:

  • As setas são plotadas em uma escala relacionada à proporção máxima entre o vetor próprio em escala de cada um dos dois componentes principais e suas respectivas pontuações em escala (o ratio). Comentários do AS @amoeba:

o gráfico de dispersão e o "gráfico de seta" são redimensionados de modo que a maior (em valor absoluto) a coordenada das setas x ou y das setas seja exatamente igual à maior (em valor absoluto) a coordenada x ou y dos pontos de dados dispersos

  • Como previsto acima, os pontos podem ser plotados diretamente como as pontuações na matriz do SVD:U

Antoni Parellada
fonte
11
+1, bom estudo. Eu adicionei tag Rà sua pergunta porque o assunto confuso (ou seja, o coeficiente de escala) provou ser parcialmente específico para R. Em geral, você pôde ver que o biplot PCA é um gráfico de dispersão de sobreposição de pontuações de componentes (coordenadas de linha) e coeficientes de direção de componente (coordenadas de coluna) e, já que várias quantidades de padronizações por "inércia" (variação) podem ser aplicadas a cada do também, podem surgir vários aspectos do biplot. Para adicionar: mais tipicamente (mais sentido), os carregamentos são exibidos como as coordenadas da coluna (setas).
ttnphns
11
(cont.) Veja minha visão geral do biplot, que explica, em palavras diferentes, o que você mostrou em sua bela resposta.
ttnphns
2
+ 1 obrigado por escrever tutoriais e com código reproduzível para nós!
Haitao Du 1/17
Antônio, você desenhou (à mão) ou plotou (alimentou dados) sua foto? Qual software você usou? Parece legal.
ttnphns
@ttnphns Obrigado! Aqui está o link para ele . Fiquei me perguntando se você poderia melhorar e plotar as cargas e os PCs de uma maneira melhor e mais didática. Sinta-se livre para mudar (é um programa maravilhosamente amigável) e, se o fizer, compartilhe.
Antoni Parellada