Princomp() produce puntuaciones PCA aparentemente erróneas con el argumento de entrada cor=TRUE

Question

Dado un marco de datos 2D d centrado y escalado:

d = data.frame(x = c(1,2,3,4,5,6,7,8,9,10), y = c(3,5,6,8,3,9,3,5,7,15))
d = as.data.frame(scale(d,center=TRUE, scale = TRUE))

la matriz de correlaciones y la matriz de covarianzas son iguales:

all.equal(cov(d),cor(d)) # this equals TRUE, meaning cov(d) and cor(d) are equal

Ahora bien, si utiliza princomp para hacer PCA y comparar las cargas/vectores propios producidos cuando se especifica utilizar la matriz de covarianza ( cor=FALSE ) frente a la matriz de correlaciones ( cor=TRUE ) se obtienen las mismas cargas esperadas:

princomp(d, cor = TRUE)$loading
princomp(d, cor = FALSE)$loading

pero las puntuaciones que se producen no son las mismas:

princomp(d, cor = TRUE)$scores
princomp(d, cor = FALSE)$scores

head(princomp(d, cor = TRUE)$scores,3)
#         Comp.1     Comp.2
# [1,] 1.7950077 -0.4206391
# [2,] 1.1445987 -0.5786822
# [3,] 0.6963027 -0.5346122

head(princomp(d, cor = FALSE)$scores,3)
#         Comp.1     Comp.2
# [1,] 1.7028938 -0.3990533
# [2,] 1.0858616 -0.5489861
# [3,] 0.6605707 -0.5071777

¿Por qué?

Las puntuaciones del PCA son sólo los datos originales, d multiplicado por los vectores propios (es decir, las cargas). La dirección d es el mismo en ambas llamadas a princomp y la matriz que se descompone es la misma ya que cov(d)=cor(d) ¿a qué se debe la diferencia de puntuación?

Calculado manualmente los datos originales multiplicados por los vectores propios coincide con la salida cuando cor=FALSE :

head(as.matrix(d) %*% as.matrix(princomp(d)$loading),3)

---

Intentar comprender el código fuente

Para ver lo que está pasando mira el código de la aplicación princomp() escribiendo

getAnywhere("princomp.default")

He reconstruido el código fuente para ver qué ocurre cuando cor=TRUE vs. cuando cor=FALSE . Aquí está el código fuente cambiado de lugar para que pueda coincidir con la salida cor=TRUE vs. cor=FALSE :

z = as.matrix(d)
sd_of_original_data = apply(d,2,sd)  # this is 1 because the data was scaled
covmat <- cov.wt(z)
n.obs <- covmat$n.obs
cv <- covmat$cov * (1 - 1/n.obs)
cen <- covmat$center

# if cor = TRUE code does this: i.e CORRELATION MATRIX code
sds <- sqrt(diag(cv))
sds   # this is .9486833 for both columns
cv_cor <- cv/(sds %o% sds)
cv_cor
edc_cor <- eigen(cv_cor, symmetric = TRUE)
ev_cor <- edc$values
sc_cor = sds
scr_cor = scale(z, center = cen, scale = sc_cor) %*% edc_cor$vectors
# reconstructed scores using source code match princomp output
scr_cor
princomp(d, center = FALSE, scale = FALSE, cor= TRUE)$scores

# if cor = FALSE code does this: i.e COVARAINCE MATRIX code
edc <- eigen(cv, symmetric = TRUE)
ev <- edc$values
sdev <- sqrt(ev)
sc = rep(1, ncol(cv)) #, colnames(cv)
scr = scale(z, center=cen, scale=sc) %*% edc$vectors
# reconstructed scores using source code match princomp output
scr
princomp(d, center=FALSE, scale=FALSE, cor=FALSE)$scores

En primer lugar, observe que el código fuente tiene esta división

 covmat$cov * (1 - 1/n.obs)

y nota covmat$cov tiene unos en la diagonal, es decir sd(covmat$cov) = 1 pero una vez que multiplicas por (1 - 1/n.obs) then sd(covmat$cov) != 1 . Esta es la cuestión. Detalles a continuación.

1. En cor=TRUE el cv (matriz de covarianza) se divide por el producto exterior de la matriz que es el cv/(sds %o% sds). Nota sds es el sqrt(diag(cv)) y sds son de 0,94 para ambas columnas de datos. La dirección sd_of_original_data tiene sds de 1 para ambas columnas porque los datos originales estaban estandarizados. En cor=FALSE esta división no tiene lugar.

   Observa cuando los datos están centrados y escalados diag(cv) no serán 1s debido a la elección original de hacer esto: cv <- covmat$cov * (1 - 1/n.obs) . Arriba cuando cor= TRUE, sds= diag(cv) devuelve valores distintos de uno, mientras que cuando cor=FALSE sds es igual a un vector de 1s. Este vector "sds" que es función de cv <- covmat$cov * (1 - 1/n.obs) es la razón por la que las puntuaciones son diferentes porque abajo las puntuaciones se escalan por sds.

2. En cor=TRUE las puntuaciones se crean a partir del escalado de la matriz de datos original, z por sc_cor que son los sds . A continuación, multiplicando por los vectores propios, edc_cor que recuerde se calcularon utilizando cv_cor que tenía la división en el paso 1.

   sc_cor =sds
   scr_cor = scale(z, center=cen, scale=sc_cor) %*% edc_cor$vectors

   vs. cuando cor = FALSE las puntuaciones se crean escalando la matriz de datos original por 1 y multiplicando por los vectores propios, edc de cv. Este parece ser el cálculo de puntuación PCA de libro de texto.

   sc = rep(1, ncol(cv)) #, colnames(cv)
   scr = scale(z, center = cen, scale = sc) %*% edc$vectors

   Así que al principio pensé que si cor(d) = cov(d) las cargas y puntuaciones del PCA serían las mismas. Obviamente, la matriz de covarianza que usted cree que utiliza princomp a veces no se utiliza realmente. En su lugar se modifica.

   Parece que para evitar toda esta confusión siempre se puede pasar la matriz de covarianza que se desea utilizar en el parámetro "covmat" a princomp. Entonces usted tendrá el control sobre la matriz de covarianza utilizada en el PCA. Si no desea que el cv escalado por (1 - 1/n.obs) sólo puede pasar el cv adecuado en "covmat".

   Observa que hay diferentes formas de escalar la matriz de covarianza. Aquí es la ayuda del cov.wt() utilizado en el código fuente.

   Por defecto, método = "insesgado", La matriz de covarianza se divide por uno menos la suma de cuadrados de los pesos, por lo que si los pesos son los predeterminados (1/n) se obtiene la estimación insesgada convencional de la matriz de covarianza con divisor (n - 1). Esto difiere del comportamiento en S-PLUS que corresponde a método = "ML" y no divide.

También una pregunta. El código fuente tiene escala de la matriz de covarianza para los pesos.

 cv <- covmat$cov * (1 - 1/n.obs)

pero estoy confundido si nos fijamos en getAnywhere("cov.wt") ya está escalando.

Here is the data and the result from cov.wt() and the source code used to calculated cov.wt.   See getAnywhere("cov.wt")

d = data.frame(x = c(1,2,3,4,5,6,7,8,9,10), y=c(3,5,6,8,3,9,3,5,7,15))
d = as.data.frame(scale(d,center=TRUE, scale=TRUE))
cov.wt(d)$cov
##### here is what cov.wt is doing see getAnywhere("cov.wt")
wt = rep(.1,10)
center = colSums(wt * d)
x <- sqrt(wt) * sweep(d, 2, center, check.margin=FALSE)
crossprod(as.matrix(x))/(1 - sum( wt^2 ))  # this is cov.wt(d)$cov

Ver el denominador en la última línea (1 - sum( wt^2 )) ... así que ¿por qué la necesidad de la escala de este resultado de nuevo por lo que ya es la escala de los pesos por lo que por (1 - 1 / n.obs)?

Consulte aquí la sección "Muestras ponderadas": https://en.wikipedia.org/wiki/Sample_mean_and_covariance

"Si todas las ponderaciones son iguales .... la media ponderada y la covarianza se reducen a la media muestral y la covarianza anteriores".

Esencialmente (1 - sum( wt^2 )) en el código fuente de cov.wt() es el mismo que el término (1 - 1/n.obs) usado para multiple cov.wt() en el código fuente de princomp() porque los pesos son los mismos, 1/n (sum(wt^2) == 1/n.obs)...¿así que están escalando dos veces?

Answer 1

No recomiendo utilizar princomp() porque es una fuente de confusión constante. Utilice prcomp() en su lugar.

Princomp() calcula la matriz de covarianza utilizando la $1/n$ en lugar del factor $1/(n-1)$ lo que significa que los valores propios devueltos por princomp y por cualquier otra función R como, por ejemplo prcomp o eigen(cov()) serán, muy confusamente, diferentes: ¿Por qué las funciones de R "princomp" y "prcomp" dan valores propios diferentes?

Tu problema está muy relacionado con eso (como ha identificado correctamente @whuber en los comentarios). Cuando princomp se ejecuta con cor=FALSE (por defecto), calcula los vectores propios de la matriz de covarianza, centra los datos y los proyecta sobre los vectores propios. Cuando se ejecuta con cor=TRUE la función calcula los vectores propios de la matriz de correlación, $z$ -valora los datos utilizando las varianzas calculadas con $1/n$ factor y proyecta los datos sobre los vectores propios.

Esto significa que las puntuaciones PCA producidas por las dos líneas siguientes serán ligeramente diferentes (porque scale utiliza $n-1$ ):

 princomp(scale(data))$scores
 princomp(data, cor=T)$scores

Esto me parece contraintuitivo y confuso. Como @ttnphns señala en los comentarios, el uso de $n$ no es necesariamente peor que utilizar $n-1$ . Pero princomp es incoherente con casi todas las demás funciones de R, lo que da lugar a una gran confusión.

Si ya suministra $z$ -datos puntuados y especificar cor=TRUE la función $z$ -vuelve a puntuar los datos, a su gusto, multiplicando cada columna por $\sqrt{(n-1)/n}$ . Esta es la diferencia que se observa.

(Véase aquí sobre una diferencia algorítmica no relacionada entre princomp y prcomp .)

---

Algunos comentarios sobre el código fuente

Veamos el código que nos ha proporcionado. La función comienza calculando la matriz de covarianza y multiplicándola por (1 - 1/n.obs) . Le estás dando demasiadas vueltas a esto; simplemente se convierte $1/(n-1)$ factor en $1/n$ porque $1-1/n=(n-1)/n$ . Así que de acuerdo, la matriz de covarianza está ahora escalada para ser de máxima verosimilitud (en lugar de insesgada).

En cor=FALSE (que es la predeterminada), las puntuaciones se calculan como

 scr =  scale(z, center = cen, scale = sc) %*% edc$vectors

que son datos centrados (sin escalar: sc es igual a 1 ) multiplicado por los vectores propios.

En cor=TRUE las puntuaciones se calculan con la misma fórmula, pero ahora el sc se establece en sqrt(diag(cv)) donde cv es la matriz de covarianza ML. Como usted mismo escribió, en su caso sc=.9486833 que es $\sqrt{9/10}$ . Así que aunque $z$ -calificó sus datos, princomp se re- $z$ -score it for you because it prefers different denominator. De ahí la confusión.