Jak obliczyć główne składniki obrócone varimax w R?

13

Uruchomiłem PCA na 25 zmiennych i wybrałem 7 najlepszych komputerów za pomocą prcomp.

prc <- prcomp(pollutions, center=T, scale=T, retx=T)

Następnie wykonałem obrót varimax na tych elementach.

varimax7 <- varimax(prc$rotation[,1:7])

A teraz chcę varimax obrócić dane obrócone PCA (ponieważ nie jest to część obiektu varimax - tylko macierz obciążeń i macierz obrotu). Przeczytałem, że w tym celu mnożymy transpozycję macierzy rotacji przez transpozycję danych, więc zrobiłbym to:

newData <- t(varimax7$rotmat) %*% t(prc$x[,1:7])

Ale to nie ma sensu, ponieważ powyższe wymiary macierzy transponują odpowiednio i więc pozostanie mi macierz składająca się tylko z wierszy, a nie wierszy ... czy ktoś wie co robię źle tutaj lub jaka powinna być moja ostatnia linia? Czy po prostu muszę później dokonać transpozycji?7 × 16933 7 169337×77×16933716933

Scott
źródło

Odpowiedzi:

22

„Obroty” to podejście opracowane w analizie czynnikowej; tam obroty (takie jak np. varimax) są stosowane do ładunków , a nie do wektorów własnych macierzy kowariancji. Ładunki są wektorami własnymi skalowanymi przez pierwiastki kwadratowe odpowiednich wartości własnych. Po rotacji varimax wektory ładujące nie są już ortogonalne (nawet jeśli obrót nazywany jest „ortogonalnym”), więc nie można po prostu obliczyć rzutów ortogonalnych danych na obrócone kierunki ładowania.

Odpowiedź FTusella zakłada, że ​​rotacja varimax jest stosowana do wektorów własnych (nie do obciążeń). To byłoby dość niekonwencjonalne. Szczegółowe informacje znajdują się na moim szczegółowym koncie PCA + varimax: Czy po PCA następuje obrót (np. Varimax) nadal PCA? W skrócie, jeśli spojrzymy na SVD macierzy danych , to obrócenie ładunków oznacza wstawienie dla niektórych macierzy rotacji w następujący sposób: R R R X = ( U R ) ( R S V ) .X=USVRRRX=(UR)(RSV).

Jeśli obrót zostanie zastosowany do obciążeń (jak to zwykle bywa), istnieją co najmniej trzy proste sposoby obliczania komputerów z rotacją varimax w R:

  1. Są one łatwo dostępne za pośrednictwem funkcji psych::principal(co pokazuje, że jest to rzeczywiście standardowe podejście). Pamiętaj, że zwraca standardowe wyniki , tzn. Wszystkie komputery mają wariancję jednostkową.

  2. Można ręcznie użyć varimaxfunkcji do obracania obciążeń, a następnie użyć nowych obróconych obciążeń, aby uzyskać wyniki; należy pomnożyć dane z transponowanym pseudo-odwrotnością obróconych ładunków (patrz wzory w tej odpowiedzi przez @ttnphns ). To da również znormalizowane wyniki.

  3. Można użyć varimaxfunkcji, aby obrócić ładunki, a następnie użyć $rotmatmacierzy obrotu, aby obrócić znormalizowane wyniki uzyskane za pomocą prcomp.

Wszystkie trzy metody dają ten sam wynik:

irisX <- iris[,1:4]      # Iris data
ncomp <- 2

pca_iris_rotated <- psych::principal(irisX, rotate="varimax", nfactors=ncomp, scores=TRUE)
print(pca_iris_rotated$scores[1:5,])  # Scores returned by principal()

pca_iris        <- prcomp(irisX, center=T, scale=T)
rawLoadings     <- pca_iris$rotation[,1:ncomp] %*% diag(pca_iris$sdev, ncomp, ncomp)
rotatedLoadings <- varimax(rawLoadings)$loadings
invLoadings     <- t(pracma::pinv(rotatedLoadings))
scores          <- scale(irisX) %*% invLoadings
print(scores[1:5,])                   # Scores computed via rotated loadings

scores <- scale(pca_iris$x[,1:2]) %*% varimax(rawLoadings)$rotmat
print(scores[1:5,])                   # Scores computed via rotating the scores

Daje to trzy identyczne wyniki:

1 -1.083475  0.9067262
2 -1.377536 -0.2648876
3 -1.419832  0.1165198
4 -1.471607 -0.1474634
5 -1.095296  1.0949536

Uwaga:varimax funkcja w R wykorzystuje normalize = TRUE, eps = 1e-5parametry domyślnie ( patrz dokumentacja ). Można chcieć zmienić te parametry (zmniejszyć epstolerancję i zadbać o normalizację Kaiser) podczas porównywania wyników z innym oprogramowaniem, takim jak SPSS. Dziękuję @GottfriedHelms za zwrócenie mojej uwagi na to. [Uwaga: te parametry działają po przekazaniu do varimaxfunkcji, ale nie działają po przekazaniu do psych::principalfunkcji. Wygląda to na błąd, który zostanie naprawiony.]

ameba mówi Przywróć Monikę
źródło
1
Widzę to teraz i myślę, że masz rację. Zmienię oryginalną odpowiedź (lub dodam inną), aby prześledzić źródło rozbieżności. Podobały mi się bardzo kompletne i oświecające odpowiedzi na twoje i @ttnphns, zawierające szczegółowe wyjaśnienia, których zwykle nie ma w książkach.
F. Tusell,
@amoeba Próbuję zrobić PCA + VariMAX użyciu principal, prcompa princomp, ale wynikające z obciążeń / wnioski badań są bardzo różne od siebie. Z tego, co rozumiem, prcomp i princomp nie zwracają znormalizowanych wyników ani ładowań. Moje pytanie brzmi: jakie jest najlepsze podejście? Czy naprawdę chcę znormalizować wyniki? Jest mój kod nie pca_iris <- prcomp(irisX, center=T, scale=T)następuje varimax(pca_iris$rotation)$loadingsprawidłowe jak ty jak powyżej?
JMarcelino,
@JMarcelino, nie, twój kod zmienia rotację varimax na wektorach własnych, a nie na ładunki. Nie tak rozumie się lub stosuje rotację varimax.
ameba mówi Przywróć Monikę
1
@JMarcelino, pytasz, dlaczego matematyka działa tak, jak mówię w metodzie nr 2? To proste, jeśli znasz ten rodzaj algebry liniowej. PCA to rozkład SVD . Zastosowanie obrotu, takiego jak varimax, oznacza wstawienie do macierzy obrotu w następujący sposób: . Obrócone ładunki to , obrócone standardowe wyniki to , więcZnasz i ; jak zdobyć ? Odpowiedź brzmi:X=USVRRRX=URRSVL=VSR/n1T=URn1
X=TL.
XLT
T=X(L)+=X(L+).
Ameba mówi Przywróć Monikę
1
Otrzymałem odpowiedź opiekuna pakietu Prof. Revelle. Wydaje się, że jest to błąd w obsłudze parametrów w principalprocedurze, który zawsze oblicza się z normalizacją Kaisera i eps = 1e-5. Jak dotąd nie ma informacji, dlaczego na r-fiddle.org wersja działa poprawnie. Powinniśmy więc czekać na aktualizacje - i powinienem usunąć wszystkie przestarzałe komentarze. ameba - dobrze byłoby odpowiednio zaktualizować uwagę w odpowiedzi. Dziękuję za całą współpracę!
Gottfried Helms
9

Musisz użyć matrycy $loadings, a nie $rotmat:

 x <- matrix(rnorm(600),60,10)
 prc <- prcomp(x, center=TRUE, scale=TRUE)
 varimax7 <- varimax(prc$rotation[,1:7])
 newData <- scale(x) %*% varimax7$loadings

Macierz $rotmatjest macierzą ortogonalną, która wytwarza nowe ładunki z nieobrotowych.

EDYTUJ od 12 lutego 2015 r .:

Jak słusznie wskazał poniżej @amoeba (patrz również jego / jej poprzedni post oraz inny post z @ttnphns ), ta odpowiedź jest nieprawidłowa. Rozważmy macierzy danych . Rozkład wartości w liczbie pojedynczej to gdzie ma za kolumny (znormalizowane) wektory własne . Teraz obrót jest zmianą współrzędnych i sprowadza się do zapisania powyższej równości jako: przy czym jest macierzą ortogonalną wybraną do uzyskania bliski rzadkości (maksymalny kontrast między wpisami, luźno mówiąc). Teraz, jeśli to wszystkon×mX

X=USVT
VXX
X=(UST)(TTVT)=UV
TV, co nie jest, można pomnożyć powyższą równość przez aby uzyskać wyniki jako , ale oczywiście nigdy nie obracamy całego komputera. Przeciwnie, uważa się, że podzbiór , który zapewnia jeszcze godnego rank- aproksymacji , tak więc rozwiązanie jest obrócony gdzie teraz jest macierzą . Nie możemy już po prostu pomnożyć przez transpozycję , ale raczej musimy skorzystać z jednego z rozwiązań opisanych przez @amoeba.VUX(V)Tk<mkXX(UkSkTk)(T T k V T k )=Uk Vk Vk k×nXVk
X(UkSk)(VkT)
X(UkSkTk)(TkTVkT)=UkVk
Vkk×nXVk

Innymi słowy, zaproponowane przeze mnie rozwiązanie jest poprawne tylko w tym konkretnym przypadku, w którym byłoby bezużyteczne i bezsensowne.

Serdeczne podziękowania dla @amoeba za wyjaśnienie mi tej sprawy; Od lat żyję z tym nieporozumieniem.

Jeden punkt, w którym uwaga powyżej odbiega od użytkownika @ ameba postu jest to, że ona / on wydaje skojarzyć z w . Myślę, że w PCA bardziej powszechne jest posiadanie kolumn normy 1 i pochłanianie w wartościach głównego składnika. W rzeczywistości zwykle są one przedstawiane jako kombinacje liniowe oryginalnych (wyśrodkowanych, być może skalowanych) zmiennych podlegających . Tak czy inaczej jest do zaakceptowania, jak sądzę, i wszystko pomiędzy (jak w analizie biplot).SVLVSviTX (i=1,,m)vi=1

DALSZA EDYCJA 12 lutego 2015 r

Jak zauważył @amoeba, mimo że jest prostokątny, zaproponowane mnie rozwiązanie może być nadal akceptowalne: dałoby macierz jednostkową i . Wszystko wydaje się więc zależeć od definicji wyników, które preferuje się. V k ( V k ) T X ( V k ) TU kVkVk(Vk)TX(Vk)TUk

F. Tusell
źródło
1
Ach racja. Byłem zdezorientowany, ponieważ ładunki dla prcomp nazywane są „rotacją”, powinien był lepiej przeczytać pomoc. Skoro używam metody „centrum = PRAWDA, skala = PRAWDA” w metodzie prcomp, czy to oznacza, że ​​naprawdę powinienem wyśrodkować i skalować moje dane przed pomnożeniem ich przez ładowanie varimax $?
Scott
1
Tak, dobra uwaga, mój błąd. Centrowanie nie miałoby znaczenia, jakby tylko przesunęło punkty, ale skala powinna być taka sama, jak używana do obliczania głównych składników, które nie są niezmienne dla skalowania.
F. Tusell
2
Zapomniałem wspomnieć, że możesz chcieć spojrzeć na funkcję factanal, jeśli jeszcze tego nie zrobiłeś. Dokonuje analizy czynnikowej, a nie głównych składników, ale zwraca wyniki bezpośrednio.
F. Tusell
2
-1. Uważam, że ta odpowiedź jest nieprawidłowa i opublikowałem własną odpowiedź, aby ją zademonstrować. Nie można uzyskać wyników obróconych za pomocą rzutowania ortogonalnego na obrócone obciążenia (ponieważ nie są one już ortogonalne). Najprostszym sposobem na uzyskanie poprawnych wyników jest użycie psych::principal. [Poza tym zredagowałem twoją odpowiedź, aby wstawić skalowanie, jak omówiono w komentarzach powyżej.]
amoeba mówi Przywróć Monikę
1
Przepraszam moja wina. Miałem na myśli, że to . Naprawię to teraz. I ... tak, teraz, kiedy na to patrzę, ma kolumny ortogonalne, więc nadal otrzyma macierz jednostek, prawda? Jeśli tak, to nie wprowadziłem w błąd oryginalnego plakatu, zdejmujesz ładunek z mojej duszy! k × n V ( T T k V T k ) ( V k T k )Vkk×nV(TkTVkT)(VkTk)
F. Tusell,
0

Szukałem rozwiązania, które działa na PCA wykonywane przy użyciu ade4 .

Poniżej znajdziesz funkcję:

library(ade4)

irisX <- iris[,1:4]      # Iris data
ncomp <- 2
# With ade4
dudi_iris <- dudi.pca(irisX, scannf = FALSE, nf = ncomp)

rotate_dudi.pca <- function(pca, ncomp = 2) {

  rawLoadings <- as.matrix(pca$c1[,1:ncomp]) %*% diag(sqrt(pca$eig), ncomp, ncomp)
  pca$c1 <- rawLoadings
  pca$li <- scale(pca$li[,1:ncomp]) %*% varimax(rawLoadings)$rotmat

  return(pca)
} 
rot_iris <- rotate_dudi.pca(pca = dudi_iris, ncomp = ncomp)
print(rot_iris$li[1:5,])                   # Scores computed via rotating the scores
#>        [,1]       [,2]
#> 1 -1.083475 -0.9067262
#> 2 -1.377536  0.2648876
#> 3 -1.419832 -0.1165198
#> 4 -1.471607  0.1474634
#> 5 -1.095296 -1.0949536

Utworzono 14.01.2020 przez pakiet reprezentx (v0.3.0)

Mam nadzieję, że to pomoże!

Alain Danet
źródło
Musisz użyć tego miejsca na odpowiedź.
Michael R. Chernick
Wydawało mi się, że poprawne jest dodanie odpowiedzi. Podobnie jak w przypadku tego pytania: stackoverflow.com/questions/6862742/draw-a-circle-w--ggplot2 . W razie potrzeby chętnie przedstawię swoją propozycję.
Alain Danet
Źle zrozumiałem, ponieważ brzmiało to tak, jakbyś poprawiał błąd w jednej z odpowiedzi. Widzę, że jest to dodatek do konkretnego pakietu oprogramowania ad4. Cross Validated nie analizuje pytań i odpowiedzi, które dotyczą wyłącznie kodu. Przepełnienie stosu polega na rozwiązywaniu problemów z oprogramowaniem.
Michael R. Chernick