Bagaimana cara menghitung komponen utama yang dirotasi varimax dalam R?

13

Saya menjalankan PCA pada 25 variabel dan memilih 7 PC teratas menggunakan prcomp.

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

Saya kemudian melakukan rotasi varimax pada komponen-komponen itu.

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

Dan sekarang saya ingin varimax memutar data yang diputar PCA (karena ini bukan bagian dari objek varimax - hanya matriks pemuatan dan matriks rotasi). Saya membaca bahwa untuk melakukan ini, Anda mengalikan transpos dari matriks rotasi dengan transpos data sehingga saya akan melakukan ini:

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

Tapi itu tidak masuk akal karena dimensi matriks transposisi di atas adalah 7×7 dan masing-masing dan jadi saya akan dibiarkan dengan matriks hanya baris, daripada baris ... apakah ada yang tahu apa Saya melakukan kesalahan di sini atau apa garis akhir saya seharusnya? Apakah saya hanya perlu memindahkan kembali setelah itu?7×16933716933

Scott
sumber

Jawaban:

22

"Rotasi" adalah pendekatan yang dikembangkan dalam analisis faktor; ada rotasi (seperti misalnya varimax) diterapkan pada pemuatan , bukan vektor eigen dari matriks kovarians. Loading adalah vektor eigen yang diskalakan oleh akar kuadrat dari masing-masing nilai eigen. Setelah rotasi varimax, vektor pemuatan tidak ortogonal lagi (meskipun rotasi disebut "ortogonal"), jadi orang tidak bisa hanya menghitung proyeksi ortogonal data ke arah pemuatan yang diputar.

@ Jawaban FTusell mengasumsikan bahwa rotasi varimax diterapkan pada vektor eigen (bukan untuk memuat). Ini akan sangat tidak konvensional. Silakan lihat akun terperinci saya PCA + varimax untuk detail: Apakah PCA diikuti oleh rotasi (seperti varimax) masih PCA? Secara singkat, jika kita melihat SVD dari matriks data X=USV , maka untuk memutar beban berarti memasukkan RR untuk beberapa matriks rotasi R sebagai berikut: X=(UR)(RSV).

Jika rotasi diterapkan ke pemuatan (seperti biasanya), maka setidaknya ada tiga cara mudah untuk menghitung PC yang dirotasi dengan varimax di R:

  1. Mereka sudah tersedia melalui fungsi psych::principal(menunjukkan bahwa ini memang pendekatan standar). Perhatikan bahwa ia mengembalikan skor terstandarisasi , yaitu semua PC memiliki varian unit.

  2. Satu dapat secara manual menggunakan varimaxfungsi untuk memutar beban, dan kemudian menggunakan beban diputar baru untuk mendapatkan skor; kita perlu melipatgandakan data dengan pseudo-invers dari pemuatan yang diputar (lihat rumus dalam jawaban ini oleh @ttnphns ). Ini juga akan menghasilkan skor standar.

  3. Seseorang dapat menggunakan varimaxfungsi untuk memutar beban, dan kemudian menggunakan $rotmatmatriks rotasi untuk memutar skor standar yang diperoleh prcomp.

Ketiga metode menghasilkan hasil yang sama:

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

Ini menghasilkan tiga keluaran identik:

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

Catatan: The varimaxfungsi dalam R menggunakan normalize = TRUE, eps = 1e-5parameter secara default ( lihat dokumentasi ). Seseorang mungkin ingin mengubah parameter ini (mengurangi epstoleransi dan menjaga normalisasi Kaiser) ketika membandingkan hasilnya dengan perangkat lunak lain seperti SPSS. Saya berterima kasih kepada @GottfriedHelms karena membawa ini menjadi perhatian saya. [Catatan: parameter ini berfungsi saat diteruskan ke varimaxfungsi, tetapi tidak berfungsi saat diteruskan ke psych::principalfungsi. Tampaknya ini adalah bug yang akan diperbaiki.]

amuba kata Reinstate Monica
sumber
1
Saya melihat ini sekarang, dan saya pikir Anda benar. Saya akan mengedit jawaban asli saya (atau menambahkan yang lain) untuk melacak sumber perbedaan tersebut. Saya menyukai jawaban Anda dan @ttnphns sangat lengkap dan mencerahkan, memberikan penjelasan terperinci yang biasanya tidak ditemukan dalam buku.
F. Tusell
@amoeba Saya mencoba melakukan varimax menggunakan PCA principal, prcompdan princomp, tetapi hasil yang diambil / kesimpulan studi sangat berbeda satu sama lain. Untuk apa yang saya mengerti, prcomp dan princomp tidak mengembalikan skor standar atau memuat. Pertanyaan saya adalah: apa pendekatan terbaik? Apakah saya benar-benar menginginkan hasil standar? Bukankah kode saya pca_iris <- prcomp(irisX, center=T, scale=T)diikuti oleh varimax(pca_iris$rotation)$loadingsbenar seperti Anda di atas?
JMarcelino
@ JMarcelino, tidak, kode Anda melakukan rotasi varimax pada vektor eigen, bukan pada pemuatan. Ini bukan bagaimana rotasi varimax biasanya dipahami atau diterapkan.
Amuba kata Reinstate Monica
1
@ JMarcelino, apakah Anda bertanya mengapa matematika bekerja seperti yang saya katakan dalam metode # 2? Sederhana jika Anda terbiasa dengan aljabar linier semacam ini. PCA adalah dekomposisi SVD . Menerapkan rotasi seperti varimax berarti menyisipkan R R untuk matriks rotasi R sebagai berikut: X = U R R S V . Pemuatan yang diputar adalah L = V S R /X=USVRRRX=URRSV , skor standar yang dirotasi adalahT=URL=VSR/n-1 , jadiX=TL. Anda tahuXdanL; bagaimana cara mendapatkanT? Nah, jawabannya adalahT=X(L)+=X(L+). T=URn-1
X=TL.
XLT
T=X(L)+=X(L+).
Amoeba berkata Reinstate Monica
1
Saya mendapat jawaban dari pengelola paket Prof. Revelle. Tampaknya menjadi bug dalam penanganan parameter dalam principalprosedur, yang selalu dihitung dengan Kaiser-normalisasi dan eps = 1e-5. Tidak ada informasi sejauh ini, mengapa di r-fiddle.org versi ini berfungsi dengan benar. Jadi kita harus menunggu pembaruan - dan saya harus menghapus semua komentar yang sudah usang. amoeba - akan lebih baik untuk memperbarui komentar dalam jawaban Anda sesuai. Terima kasih atas semua kerjasamanya!
Gottfried Helms
9

Anda harus menggunakan matriks $loadings, bukan$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

Matriks $rotmat adalah ortogonal yang menghasilkan pemuatan baru dari yang tidak diputar.

EDIT pada 12 Februari 2015:

Sebagai benar menunjukkan bawah oleh @amoeba (lihat juga / nya posting sebelumnya serta pos lain dari @ttnphns ) jawaban ini tidak benar. Pertimbangkan matriks data X . Dekomposisi nilai singular adalah X = U S V T dimana V memiliki sebagai kolom nya yang (dinormalisasi) vektor eigen dari X ' X . Sekarang, rotasi adalah perubahan koordinat dan jumlah untuk menulis persamaan di atas sebagai: X = ( U S T ) ( T T V T =n×mX

X=USVT
VXX dengan T menjadi matriks ortogonal yang dipilih untuk mencapai V mendekati sparse (kontras maksimum antara entri, secara longgar). Sekarang,jika hanya itu, yang bukan, kita bisa mengalikan kesetaraan di atas dengan V untuk mendapatkan skor U sebagai X ( V ) T , Tapi tentu saja kita tidak pernah merotasi semua PC. Sebaliknya, kami menganggap himpunan bagian dari k < m yang masih memberikanperkiraanperingkat k yang layakdari X , X
X=(UST)(TTVT)=UV
TVVUX(V)Tk<mkX sehingga solusi yang diputar sekarang X ( U k S k T k ) ( T T k V T k ) = U k V k di mana sekarang V k adalah k × n matriks. Kita tidak bisa lagi hanya kalikan X dengan transpos dari V * k
X(UkSk)(VkT)
X(UkSkTk)(TkTVkT)=UkVk
Vkk×nXVk, tetapi kita perlu menggunakan salah satu solusi yang dijelaskan oleh @amoeba.

Dengan kata lain, solusi yang saya usulkan hanya benar dalam kasus tertentu di mana itu akan sia-sia dan tidak masuk akal.

Terima kasih yang tulus pergi ke @amoeba untuk menjelaskan masalah ini kepada saya; Saya telah hidup dengan kesalahpahaman ini selama bertahun-tahun.

SVLVSviTX (i=1,,m)vi=1. Either way diterima saya pikir, dan segala sesuatu di antaranya (seperti dalam analisis biplot).

EDIT LEBIH LANJUT 12 Februari 2015

VkVk(Vk)TX(Vk)TUk

F. Tusell
sumber
1
Ah benar besar. Saya jadi bingung karena memuat untuk prcomp disebut "rotasi", seharusnya membaca bantuan lebih baik. Karena saya menggunakan "center = TRUE, skala = TRUE" dalam metode prcomp apakah itu berarti bahwa saya benar-benar harus memusatkan dan menskala data saya sebelum mengalikannya dengan $ loadm varimax saya?
Scott
1
Ya, poin bagus, kesalahanku. Pemusatan tidak akan masalah, seolah-olah hanya akan menggeser titik, tetapi skala harus sama dengan yang digunakan untuk menghitung komponen utama, yang tidak berbeda dengan penskalaan.
F. Tusell
2
Saya lupa menyebutkan bahwa Anda mungkin ingin melihat fungsi factanal, jika Anda belum melakukannya. Ia melakukan analisis faktor daripada komponen utama, tetapi akan mengembalikan skor secara langsung.
F. Tusell
2
-1. Saya percaya bahwa jawaban ini tidak benar dan saya memposting jawaban saya sendiri untuk menunjukkannya. Seseorang tidak bisa mendapatkan skor yang diputar dengan proyeksi ortogonal pada pemuatan yang diputar (karena mereka tidak ortogonal lagi). Cara paling sederhana untuk mendapatkan skor yang benar adalah menggunakan psych::principal. [Selain itu, saya mengedit jawaban Anda untuk menyisipkan penskalaan, seperti yang dibahas dalam komentar di atas.]
amoeba berkata Reinstate Monica
1
Vkk×nV(TkTVkT)(VkTk)
F. Tusell
0

Saya sedang mencari solusi yang berfungsi untuk PCA dilakukan menggunakan ade4 .

Silakan temukan fungsi di bawah ini:

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

Dibuat pada 2020-01-14 oleh paket reprex (v0.3.0)

Semoga bantuan ini!

Alain Danet
sumber
Anda perlu menggunakan ruang ini untuk mendapat jawaban.
Michael R. Chernick
Bagi saya tampaknya sah untuk menambahkan jawaban untuk kelengkapan. Seperti untuk pertanyaan ini: stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2 . Saya akan senang untuk memindahkan proposisi saya jika perlu.
Alain Danet
Saya salah paham karena sepertinya Anda melakukan koreksi terhadap salah satu jawaban. Saya melihat bahwa itu adalah tambahan untuk paket perangkat lunak tertentu ad4. Cross Validated tidak melihat pertanyaan atau jawaban yang sepenuhnya tentang kode. Stack Overflow adalah tempat masalah perangkat lunak ditangani.
Michael R. Chernick