Como calcular os principais componentes rotacionados em varimax em R?

Executei o PCA em 25 variáveis e selecionei os 7 principais PCs usando prcomp.

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

Eu fiz a rotação varimax nesses componentes.

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

E agora desejo varimax rotacionar os dados rotacionados com PCA (como ele não faz parte do objeto varimax - apenas a matriz de cargas e a matriz de rotação). Eu li que, para fazer isso, você multiplica a transposição da matriz de rotação pela transposição dos dados, então eu teria feito isso:

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

Mas isso não faz sentido, pois as dimensões da matriz transposta acima são e respectivamente, e por isso ficarei com uma matriz de apenas linhas, em vez de linhas ... alguém sabe o que eu estou fazendo errado aqui ou qual deve ser minha linha final? Eu só preciso transpor depois? $7\times 7$ $7 \times 16933$ $7$ $16933$

r pca factor-rotation Scott
fonte

Respostas:

"Rotações" é uma abordagem desenvolvida na análise fatorial; rotações (como, por exemplo, varimax) são aplicadas a cargas , não a autovetores da matriz de covariância. As cargas são vetores próprios dimensionados pelas raízes quadradas dos respectivos valores próprios. Após a rotação varimax, os vetores de carregamento não são mais ortogonais (mesmo que a rotação seja chamada "ortogonal"), portanto, não é possível simplesmente calcular projeções ortogonais dos dados nas direções de carregamento giradas.

A resposta da @ FTusell assume que a rotação varimax é aplicada aos autovetores (não às cargas). Isso seria bem pouco convencional. Consulte a minha conta detalhada do PCA + varimax para obter detalhes: O PCA seguido de uma rotação (como o varimax) ainda é PCA? Resumidamente, se olharmos para o SVD da matriz de dados , girar as cargas significa inserir para alguma matriz de rotação seguinte maneira: $X=USV^\top$ $RR^\top$ $R$ $X=(UR)(R^\top SV^\top).$

Se a rotação é aplicada às cargas (como normalmente é), há pelo menos três maneiras fáceis de calcular PCs com rotação varimax no R:

Eles estão prontamente disponíveis via função psych::principal(demonstrando que essa é realmente a abordagem padrão). Observe que ele retorna pontuações padronizadas , ou seja, todos os PCs têm variação de unidade.
Pode-se usar manualmente a varimaxfunção para girar as cargas e, em seguida, usar as novas cargas giradas para obter as pontuações; é necessário multiplicar os dados com o pseudo-inverso transposto das cargas rotacionadas (consulte as fórmulas nesta resposta por @ttnphns ). Isso também produzirá pontuações padronizadas.
Pode-se usar a varimaxfunção para girar as cargas e, em seguida, usar a $rotmatmatriz de rotação para girar as pontuações padronizadas obtidas com prcomp.

Todos os três métodos produzem o mesmo resultado:

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

Isso produz três saídas idênticas:

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

Nota: A varimaxfunção em R usa normalize = TRUE, eps = 1e-5parâmetros por padrão ( consulte a documentação ). Pode-se alterar esses parâmetros (diminuir a epstolerância e cuidar da normalização do Kaiser) ao comparar os resultados com outros softwares, como o SPSS. Agradeço a @GottfriedHelms por trazer isso à minha atenção. [Nota: esses parâmetros funcionam quando passados para a varimaxfunção, mas não funcionam quando passados para a psych::principalfunção. Parece ser um bug que será corrigido.]

ameba diz Restabelecer Monica
fonte

Eu vejo isso agora e acho que você está correto. Vou editar minha resposta original (ou adicionar outra) para rastrear a fonte da discrepância. Gostei das suas respostas @ttnphns e muito completas e esclarecedoras, fornecendo explicações detalhadas que geralmente não são encontradas em livros.

F. Tusell 12/02

@amoeba Estou tentando fazer um PCA + varimax usando principal, prcompe princomp, mas as conclusões resultantes do estudo / carregamento são muito diferentes umas das outras. Pelo que entendi, prcomp e princomp não retornam pontuações nem cargas padronizadas. Minha pergunta é: qual é a melhor abordagem? Eu realmente quero resultados padronizados? O meu código não é pca_iris <- prcomp(irisX, center=T, scale=T)seguido varimax(pca_iris$rotation)$loadingstão correto quanto o seu acima?

JCelebrino1

@JMarcelino, não, seu código faz a rotação varimax nos autovetores, não nas cargas. Não é assim que a rotação varimax é geralmente entendida ou aplicada.

Ameba diz Reinstate Monica

@JMarcelino, você está perguntando por que a matemática funciona como eu digo no método 2? É simples se você estiver familiarizado com esse tipo de álgebra linear. PCA é decomposição SVD

. Aplicando uma rotação tais como meios Varimax inserindo

para uma rotação matriz

como se segue:

. Cargas rotadas são

X = U S V^{⊤}

$X=USV^\top$

R R^{⊤}

$RR^\top$

R

$R$

X = U R R^{⊤} S V^{⊤}

$X=URR^\top SV^\top$

, os escores padronizados rotacionados são

L = V S R / \sqrt{n - 1}

$L=VSR/\sqrt{n-1}$

, então

Você conhece

; como obter

? Bem, a resposta é

T = U R \sqrt{n - 1}

$T=UR\sqrt{n-1}$

X = T {eu}^{⊤} .

$X=TL^\top.$

X

$X$

L

$L$

T

$T$

T = X ({eu}^{⊤})^{+} = X ({eu}^{+})^{⊤} .

$T=X(L^\top)^+ = X(L^+)^\top.$

Ameba diz Reinstate Monica

Eu recebi uma resposta do mantenedor do pacote Prof. Revelle. Parece haver um erro no manuseio dos parâmetros no principalprocedimento, que sempre calcula com normalização de Kaiser e eps = 1e-5. Até o momento, não há informações, por que no r-fiddle.org a versão funciona corretamente. Portanto, devemos aguardar atualizações - e eu devo excluir todos os comentários agora obsoletos. ameba - seria bom atualizar a observação em sua resposta de acordo. Obrigado por toda a cooperação!

Gottfried Helms

Você precisa usar a matriz $loadings, não $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

A matriz $rotmaté a matriz ortogonal que produz as novas cargas a partir das não rotacionadas.

EDIT em 12 de fevereiro de 2015:

$n\times m$ $X$

X = U S V^{T}

$X = USV^T$

V

$V$

X^{'} X

$X'X$

X = (U S T) (T^{T} V^{T}) = U^{*} V^{*}

$X = (UST)(T^TV^T) = U^*V^*$

T

$T$

V^{*}

$V^*$

V^{*}

$V^*$

U^{*}

$U^*$

X (V^{*})^{T}

$X(V^*)^T$

k < m

$k<m$

k

$k$

X

$X$

X \approx (U_{k} S_{k}) (V_{k}^{T})

$X \approx (U_kS_k)(V_k^T)$

X \approx (U_{k} S_{k} T_{k}) (T_{k}^{T} V_{k}^{T}) = U_{k}^{*} V_{k}^{*}

$X \approx (U_kS_kT_k)(T_k^TV_k^T) = U_k^*V_k^*$

V_{k}^{*}

$V_k^*$

k \times n

$k\times n$

X

$X$

V_{k}^{*}

$V_k^*$ , but rather we need to resort to one of the solutions described by @amoeba.

In other words, the solution I proposed is only correct in the particular case where it would be useless and nonsensical.

Heartfelt thanks go to @amoeba for making clear this matter to me; I have been living with this misconception for years.

One point where the note above departs from @amoeba's post is that she/he seems to associate $S$ with $V$ in $L$ . I think in PCA it is more common to have $V$ 's columns of norm 1 and absorb $S$ in the principal component's values. In fact, usually those are presented as linear combinations $v_i^TX$ $(i=1,\ldots,m)$ of the original (centered, perhaps scaled) variables subject to $\|v_i\|=1$ . Either way is acceptable I think, and everything in between (as in biplot analysis).

FURTHER EDIT Feb. 12, 2015

As pointed out by @amoeba, even though $V_k^*$ is rectangular, the solution I proposed might still be acceptable: $V_k^*(V_k^*)^T$ would give a unit matrix and $X(V_k^*)^T \approx U_k^*$ . So it all seems to hinge on the definition of scores that one prefers.

F. Tusell
fonte

Ah right grand. I got confused because the loadings for the prcomp are called "rotation", should have read the help better. Since I am using "center=TRUE,scale=TRUE" in the prcomp method does that mean that really I ought to be centering and scaling my data before multiplying it by my varimax$loadings?

Scott

Yes, good point, my mistake. Centering would not matter, as if only would shift the points, but the scale should be the same used to compute the principal components, which are not invariant to the scaling.

F. Tusell

I forgot to mention that you might want to look at function factanal, if you have not done it already. It does factor analysis rather than principal components, but will return the scores directly.

F. Tusell

-1. I believe that this answer is not correct and I posted my own answer to demonstrate it. One cannot get rotated scores by orthogonal projection on the rotated loadings (because they are not orthogonal anymore). The simplest way to obtain the correct scores is to use psych::principal. [Apart from that, I edited your answer to insert the scaling, as discussed in the comments above.]

amoeba says Reinstate Monica

Sorry, my bad. I meant

V_{k}^{*}

$V_k^*$ is

k \times n

$k\times n$ . I will correct it now. And... yes, now that I look at it,

V

$V$ has orthogonal columns so

(T_{k}^{T} V_{k}^{T}) (V_{k} T_{k})

$(T_k^TV_k^T)(V_kT_k)$ would still get us a unit matrix, right? If so, I did not misled the original poster, you lift a load from my soul!

F. Tusell

I was looking for a solution that works for PCA performed using ade4.

Please find the function below:

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

^{Created on 2020-01-14 by the reprex package (v0.3.0)}

Hope this help!

Alain Danet
fonte

You need to use this space for an answer.

Michael R. Chernick

It seemed to me that it is valid to add an answer for completeness. Like for this question: stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2. I will be happy to move my proposition if necessary.

Alain Danet

I misunderstood because it sounded like you were making a correction to an error in one of the answers. I see that it is an addition for a particular software package ad4. Cross Validated doesn't look at questions or answers that are strictly about code. Stack Overflow is where software issues are addressed.

Michael R. Chernick