Por que as funções R 'princomp' e 'prcomp' fornecem diferentes autovalores?

22

Você pode usar o conjunto de dados de decatlo {FactoMineR} para reproduzir isso. A questão é por que os autovalores computados diferem daqueles da matriz de covariância.

Aqui estão os autovalores usando princomp:

> library(FactoMineR);data(decathlon)
> pr <- princomp(decathlon[1:10], cor=F)
> pr$sd^2
      Comp.1       Comp.2       Comp.3       Comp.4       Comp.5       Comp.6 
1.348073e+02 2.293556e+01 9.747263e+00 1.117215e+00 3.477705e-01 1.326819e-01 
      Comp.7       Comp.8       Comp.9      Comp.10 
6.208630e-02 4.938498e-02 2.504308e-02 4.908785e-03 

E o mesmo usando PCA:

> res<-PCA(decathlon[1:10], scale.unit=FALSE, ncp=5, graph = FALSE)
> res$eig
          eigenvalue percentage of variance cumulative percentage of variance
comp 1  1.348073e+02           79.659589641                          79.65959
comp 2  2.293556e+01           13.552956464                          93.21255
comp 3  9.747263e+00            5.759799777                          98.97235
comp 4  1.117215e+00            0.660178830                          99.63252
comp 5  3.477705e-01            0.205502637                          99.83803
comp 6  1.326819e-01            0.078403653                          99.91643
comp 7  6.208630e-02            0.036687700                          99.95312
comp 8  4.938498e-02            0.029182305                          99.98230
comp 9  2.504308e-02            0.014798320                          99.99710
comp 10 4.908785e-03            0.002900673                         100.00000

Você pode me explicar por que os autovalores calculados diretamente diferem desses? (os autovetores são os mesmos):

> eigen(cov(decathlon[1:10]))$values
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

Além disso, o prcompmétodo alternativo fornece os mesmos valores próprios da computação direta:

> prc <- prcomp(decathlon[1:10])
> prc$sd^2
 [1] 1.381775e+02 2.350895e+01 9.990945e+00 1.145146e+00 3.564647e-01
 [6] 1.359989e-01 6.363846e-02 5.061961e-02 2.566916e-02 5.031505e-03

Por que PCA/ princompe prcompfornecer valores próprios diferentes?

George Dontas
fonte
O PCA fornecerá resultados diferentes dependendo de você usar a matriz de covariância ou a matriz de correlação.
charles.y.zheng
7
As diferenças parecem relativamente pequenas, embora provavelmente grandes demais para serem simples questões numéricas. Poderia ser a diferença entre normalizar por ou n - 1 , por exemplo, ao calcular uma estimativa da covariância antes de calcular a decomposição de SVD ou autovalores? nn1
cardeal
7
@ Cardinal Bom palpite! Observe que as duas seqüências diferentes de valores próprios têm proporções sucessivas idênticas. Assim, um conjunto é um múltiplo constante do outro. O múltiplo é 1.025 = 41/40 ( exatamente ). Não está claro para mim de onde isso vem. Talvez o conjunto de dados tenha 41 elementos e o OP esteja revelando apenas os 10 primeiros?
whuber
7
@ cardinal De fato: Página de ajuda para princomp: "Observe que o cálculo padrão usa o divisor N para a matriz de covariância". Página de ajuda para prcomp: "Ao contrário do princomp, as variações são calculadas com o divisor usual N-1".
caracal
2
@caracal, copie seu comentário em uma resposta (e talvez faça CW) para que possa ser aceito e a pergunta possa ser marcada como resolvida.
cardinal

Respostas:

16

princompNprcompcovN-1N

Isso é mencionado na seção Detalhes de help(princomp):

Observe que o cálculo padrão usa o divisor 'N' para a matriz de covariância.

e a seção Detalhes de help(prcomp):

Diferentemente princomp, as variações são calculadas com o divisor usual N - 1.

princompNn.obscv

else if (is.null(covmat)) {
    dn <- dim(z)
    if (dn[1L] < dn[2L]) 
        stop("'princomp' can only be used with more units than variables")
    covmat <- cov.wt(z)
    n.obs <- covmat$n.obs
    cv <- covmat$cov * (1 - 1/n.obs)
    cen <- covmat$center
}

Você pode evitar essa multiplicação especificando o covmatargumento em vez do xargumento.

princomp(covmat = cov(iris[,1:4]))$sd^2

Atualização sobre as pontuações do PCA:

Você pode configurar cor = TRUEem sua chamada para princompexecutar o PCA na matriz de correlação (em vez da matriz de covariância). Isso fará com princompquez- pontue os dados, mas ele ainda usará N para o denominador.

Como resultado, princomp(scale(data))$scorese princomp(data, cor = TRUE)$scoresdiferirá pelo fator(N-1)/N.

Joshua Ulrich
fonte
1
Você pode substituir "adivinhado" por "já confirmado" (consulte o fluxo de comentários acima.) Você também pode editar sua resposta para torná-la CW. Felicidades.
cardeal
@ cardinal Eu não vi esses comentários. Eu só vi aqueles que foram votados. Obrigado. Além disso, você poderia explicar a lógica por trás da resposta CW? Quais são as regras / diretrizes para isso?
21711 Joshua Ulrich
Alguém pode adivinhar por que o código não está simplesmente cv <- cov.wt(z, method="ML")tornando as duas linhas seguintes desnecessárias?
Caracal
2
@ Josué: Minha sugestão em relação à resposta CW deveu-se ao fato de que a resposta apareceu por meio de um fluxo de comentários e foi gerada por uma discussão de "comunidade". Como foi resolvido nos comentários, penso que faz mais sentido reformulá-lo como resposta, marcado como CW para indicar essa colaboração, e isso permite que a resposta seja aceita e a pergunta seja marcada como resolvida. (Caso contrário, ele será automaticamente ter colidido de volta pelo software depois de um determinado período de tempo.)
cardeal
2
@amoeba, seria útil mencionar isso no seu comentário de edição. "Adicionado 860 caracteres ao corpo" para uma resposta de ~ 450 caracteres não ajuda ninguém a avaliar se a edição é razoável.
27516 Joshua Ulrich