Strzały zmiennych podstawowych w biplocie PCA w R.

11

Ryzykując, że pytanie będzie specyficzne dla oprogramowania i pod pretekstem jego wszechobecności i osobliwości, chcę zapytać o funkcję biplot()w R, a dokładniej o obliczenie i wykreślenie domyślnych, nałożonych na siebie czerwonych strzałek, odpowiadających do podstawowych zmiennych.


[Aby nadać sens niektórym komentarzom, początkowo opublikowane wykresy miały problem z brakiem zainteresowania, a teraz zostały usunięte.]

Antoni Parellada
źródło
Nie mogę rozgryźć, jak naprawdę dostałeś swoje zielone strzały. Nie są poprawne. Fakt, że s.length zielony wynosi ok. dwa razy dłużej niż s.width green pozwala podejrzewać, że kreślisz wektory odnoszące się do zmiennych niestandaryzowanych. Nie może się to zdarzyć na dwójkowej PCA opartej na korelacjach.
ttnphns
Czerwone strzałki wydają się prawidłowe. Zobacz: mają tę samą długość i są symetryczne w stosunku do PC2. Jest to jedyna możliwa pozycja, gdy wykonujesz PCA przy użyciu tylko 2 zmiennych i opartych na korelacjach (tj. Zmiennych standaryzowanych). W PCA opartym na korelacjach obciążenia (koordynowane strzałkami) są korelacjami między komputerami PC a zmiennymi. W przykładzie, aby obciążenia są (vars przez komputery) .74752, .66424; -.74752, .66424.
ttnphns
@ttnphns Tak, czerwone strzałki są tym, co próbuję odtworzyć (są poprawne) i są narysowane w R wraz z biplot(name_of_the_PCA)wywołaniem, co w tym przypadku jest biplot(PCA). Wyśrodkowałem i skalowałem dane.
Antoni Parellada
Więc jakie jest twoje pytanie? Jak obliczyć współrzędne czerwonych strzałek? Powinny to być ładunki PCA . Czasami wykreślane są wektory własne (twoje polecenie R prawdopodobnie to zrobiło?), Jednak zgodnym, znaczącym sposobem jest wykreślanie ładunków .
ttnphns
@ttnphns Rysowanie wektorów własnych (zakładam, że jest takie samo jak ładunki) daje mi odpowiednią orientację (dziękuję), ale nie taką samą jasność jak czerwone strzałki (wklejam obraz w OP).
Antoni Parellada

Odpowiedzi:

19

Rozważ upvoting posta @ amoeba i @ttnphns . Dziękujemy zarówno za pomoc, jak i za pomysły.


Poniższa opiera się na zbiorze Iris R , a zwłaszcza pierwsze trzy zmienne (kolumny) Sepal.Length, Sepal.Width, Petal.Length.

Biplot łączy fabułę ładowania (unstandardized wektory własne) - w betonie, pierwsze dwa obciążenia i działki wynik (obrócone i rozszerzone Punkty danych wykreślono w stosunku do głównych składników). Korzystając z tego samego zestawu danych, @amoeba opisuje 9 możliwych kombinacji bipotu PCA w oparciu o 3 możliwe normalizacje wykresu punktowego pierwszego i drugiego głównego składnika oraz 3 normalizacje wykresu obciążenia (strzałki) zmiennych początkowych. Aby zobaczyć, jak R obsługuje te możliwe kombinacje, warto przyjrzeć się biplot()metodzie:


Najpierw algebra liniowa gotowa do skopiowania i wklejenia:

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. Odtwarzanie wykresu załadunku (strzałki):

Tutaj bardzo pomaga interpretacja geometryczna tego postu autorstwa @ttnphns . Zachowano zapis diagramu w poście: oznacza zmienną w przestrzeni tematycznej . jest ostatecznie wykreśloną odpowiednią strzałką; a współrzędne i że składnik ładuje zmienną w odniesieniu do i :h ' 1 2 V PC 1 PC 2VSepal L.ha1a2VPC1PC2


wprowadź opis zdjęcia tutaj


Składnikiem zmiennej Sepal L.w odniesieniu do będzie wówczas:PC1

a1=hcos(ϕ)

które, jeśli wyniki w odniesieniu do - nazwijmy je - są znormalizowane, aby ichPC1S1

S1=1nscores12=1 , powyższe równanie jest równoważne iloczynowi :VS1

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

Ponieważ ,V=x2

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

Również,

S1=1=var(S1)n1.

Wracając do Eq. ,(1)

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

r n - 1cos(ϕ) może zatem być uważana za współczynnik korelacji Pearsona , , z zastrzeżeniem, że nie rozumie się zmarszczek w czynnik.rn1

Kopiowanie i nakładanie się na niebiesko czerwonych strzałek 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)

wprowadź opis zdjęcia tutaj

Punkty zainteresowania:

  • Strzałki można odtworzyć jako korelację pierwotnych zmiennych z wynikami wygenerowanymi przez dwa pierwsze główne składniki.
  • Alternatywnie można to osiągnąć jak na pierwszym wykresie w drugim rzędzie, oznaczonym w poście @ amoeba:VS

wprowadź opis zdjęcia tutaj

lub w kodzie 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)

a nawet jeszcze ...

    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)

łącząc się z geometrycznym objaśnieniem ładunków przez @ttnphns lub innym postem informacyjnym również przez @ttnphns .

  • Istnieje czynnik skalowania: sqrt(nrow(X) - 1)który pozostaje nieco tajemnicą.

  • 0.8 ma związek z tworzeniem miejsca na etykietę - zobacz ten komentarz tutaj :

Ponadto należy powiedzieć, że strzałki są tak wykreślone, że środek etykiety tekstowej jest tam, gdzie powinien być! Strzałki są następnie mnożone przez 0,80,8 przed wykreśleniem, tzn. Wszystkie strzałki są krótsze niż powinny, prawdopodobnie po to, aby zapobiec nakładaniu się etykiety tekstowej (patrz kod dla biplot.default). Uważam, że jest to bardzo mylące. - ameba 19 marca 2015 o 10:06


2. Wykreślenie biplot()wykresu wyników (i strzałek jednocześnie):

Osie są skalowane do sumy jednostkowej kwadratów, odpowiadającej pierwszemu wykresowi pierwszego rzędu na słupku @ amoeba , który można odtworzyć wykreślając macierz rozkładu svd (więcej na ten temat później) - " Kolumny : są to główne składniki skalowane do jednostkowej sumy kwadratów.UUU

wprowadź opis zdjęcia tutaj

Istnieją dwie różne skale podczas gry na dolnej i górnej poziomej osi w konstrukcji biplota:

wprowadź opis zdjęcia tutaj

Jednak skala względna nie jest od razu oczywista, co wymaga zagłębienia się w funkcje i metody:

biplot()wykreśla wyniki jako kolumny w SVD, które są ortogonalnymi wektorami jednostek: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

Podczas gdy prcomp()funkcja w R zwraca wyniki skalowane do ich wartości własnych:

> 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

Dlatego możemy przeskalować wariancję do , dzieląc przez wartości własne: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

Ale ponieważ chcemy, aby suma kwadratów wynosiła , musimy podzielić przez ponieważ: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

Warto zauważyć, że użycie współczynnika skalowania zostało później zmienione na gdy zdefiniowanie wyjaśnienia wydaje się polegać na tym, żen1nlan

prcompużywa : „W przeciwieństwie do princomp, wariancje są obliczane za pomocą zwykłego dzielnika ”.n - 1n1n1


Po usunięciu ich ze wszystkich ifoświadczeń i innych puchów do czyszczenia domu biplot()postępuje w następujący sposób:

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)

który, zgodnie z oczekiwaniami, odtwarza (prawy obraz poniżej) dane biplot()wyjściowe wywoływane bezpośrednio biplot(PCA)(lewy wykres poniżej) we wszystkich jego nietkniętych niedociągnięciach estetycznych:

wprowadź opis zdjęcia tutaj

Punkty zainteresowania:

  • Strzałki są wykreślane w skali odnoszącej się do maksymalnego stosunku między skalowanym wektorem własnym każdego z dwóch głównych składników i odpowiadającymi im skalowanymi punktami (the ratio). Komentarze AS @amoeba:

wykres punktowy i „wykres strzałkowy” są skalowane w taki sposób, że największa (w wartości bezwzględnej) współrzędna strzałki x lub y strzałek była dokładnie równa największej (w wartości bezwzględnej) współrzędnej x lub y rozrzuconych punktów danych

  • Zgodnie z przewidywaniami powyżej, punkty można narysować bezpośrednio jako wyniki w macierzy SVD:U

Antoni Parellada
źródło
1
+1, dobre badanie. Dodałem tag Rdo twojego pytania, ponieważ myląca materia (mianowicie współczynnik skalowania) okazała się częściowo specyficzna dla R. Zasadniczo można było zobaczyć, że dwupłat PCA jest nakładkowym wykresem rozproszenia wyników składowych (współrzędnych wiersza) i współczynników kierunku składowych (współrzędnych kolumn), a ponieważ do każdego z nich można zastosować różne wielkości standaryzacji według „bezwładności” (wariancji) z tego, więc mogą powstać różne spojrzenia na dwupłat. Aby dodać: najczęściej (bardziej sensownie), ładunki są wyświetlane jako współrzędne kolumny (strzałki).
ttnphns
1
(cd.) Zobacz mój przegląd biplota, który wyjaśnia innymi słowy, co pokazałeś w swojej miłej odpowiedzi.
ttnphns
2
+ 1 dzięki za napisanie samouczków i powtarzalny kod dla nas!
Haitao Du
Antony, czy narysowałeś (ręcznie) czy nakreśliłeś (wprowadziłeś dane) swoje zdjęcie? Z jakiego oprogramowania korzystałeś? Wygląda ładnie.
ttnphns
@ttnphns Dziękujemy! Oto link do tego . Zastanawiałem się, czy mógłbyś to poprawić, a także przedstawić ładunki i komputery w lepszy, bardziej dydaktyczny sposób. Zapraszam do zmiany (jest to niezwykle przyjazny dla użytkownika program), a jeśli tak, prosimy o udostępnienie.
Antoni Parellada,