Práctica PCA
Alex Rodriguez Just
22/5/2020
Librerías a usar.
library(R.matlab) # Carga de datos
library(caret) # Entranar modelosCarga de datos.
data <- readMat("/Users/alexrodriguez/Desktop/data_tower_v42m.mat")1. [10 %] Escalar los datos de la matriz data$Y tal y como se explica en la Sección 2.1 del módulo. Llamar a la matriz resultante X.
–> Escalando los datos se consigue que la media sea 0 y la desviación típica 1. Pero primero visualizaré los datos para ver que son del tipo double. Selecciono las primeras cuatro filas junto a las primeras cuatro columnas:
data$Y[1:4,1:4]## [,1] [,2] [,3] [,4]
## [1,] 0.0003183557 2.463096e-04 -9.071952e-05 2.621842e-04
## [2,] 0.0003940651 2.426462e-04 7.225464e-04 4.245931e-04
## [3,] -0.0001346798 9.428016e-05 -1.444488e-04 2.833968e-05
## [4,] 0.0001980754 2.786693e-04 6.069936e-05 3.354514e-04Ahora escalo data$Y. Con la función “center” estoy restando a todos los elementos de la misma columna su media aritmética. Con la función “scale” estoy dividiendo a todos los elementos de la misma columna por su desviación típica.
X<-scale(data$Y, center = TRUE, scale = TRUE)
X[1:4,1:4]## [,1] [,2] [,3] [,4]
## [1,] 0.51711989 0.2728762 -1.1844687 0.3114625
## [2,] 0.85812155 0.2551831 2.1621197 1.0385974
## [3,] -1.52339002 -0.4613866 -1.4055647 -0.7355023
## [4,] -0.02463274 0.4291651 -0.5613802 0.63949332. [10 %] Calcular la matriz de covarianzas tal y como se explica en la Sección 2.2 del móodulo. Llamar CX a la matriz de covarianzas. ¿Os ha sido posible el cálculo de la matriz de covarianzas? ¿O habéis tenido un error por falta de memoria? En este caso, ir al problema 4. Os haya salido o no el cálculo de la matriz de covarianzas, buscar la manera de calcular CX[2,3].
CX <- cov(X)–> Al crear la matriz de covarianza me aparece mensaje de error por límite de memoria. Esto es devido a que el almacentamiento de un valor en formato double necesita 8 bytes y al intentar crear una matriz de covarianza grande como nuestra matriz “X” de 328 filas × 58008 columnas de tipo double, nos da un resultado de 58008x58008x8(Byte) / 1024x1024x1024 = 25.07Gb y mi ordenador no está preparado para soportar tanto almacenamiento.
El enunciado pide calcular CX[2,3]. Entiendo que se refiere a calcular la covarianza entre las columnas 2 y 3 de la matriz X.
cov(X[,2], X[,3])## [1] -0.0037817883. [10 %] Si habéis podido calcular la matriz de covarianzas, calcular sus valores propios y sus vectores propios. Llamar P a la matriz que tiene por columnas los vectores propios de la matriz de covarianzas, ordenados según la magnitud de los valores propios. ¿Os ha sido posible el cáalculo de los valores y de los vectores propios? ¿O habéis tenido un error por falta de memoria? Justificar vuestra respuesta. En este caso, ir al problema 4.
–> No he podido calcular la matriz de covarianza de X, por eso no puedo calcular sus valores y vectores propios. Pero podría mostrar cómo hacerlo con la misma matriz de menos variables, por ejemplo 10 variables y las 328 observaciones. Al hacer la matriz de covarianza me quedará una matriz 10x10 donde la suma de su traza es el total de sus variables.
CX <-cov(X[,1:10])
CX## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000000 0.006801465 0.222068913 0.10729467 -0.48723965
## [2,] 0.006801465 1.000000000 -0.003781788 0.28540070 -0.20321626
## [3,] 0.222068913 -0.003781788 1.000000000 -0.08063357 0.08010690
## [4,] 0.107294674 0.285400704 -0.080633567 1.00000000 -0.02848177
## [5,] -0.487239648 -0.203216258 0.080106901 -0.02848177 1.00000000
## [6,] 0.083096883 -0.225422937 0.014336486 0.30419338 -0.01812861
## [7,] -0.331836022 -0.040832353 -0.411922307 -0.06031982 0.35488767
## [8,] -0.186524120 -0.353126440 0.010978037 -0.40348745 0.18582667
## [9,] 0.301171900 0.046649831 -0.182294170 -0.18685764 -0.30078248
## [10,] -0.068177283 0.096147745 0.055946862 -0.25027768 -0.12696582
## [,6] [,7] [,8] [,9] [,10]
## [1,] 0.08309688 -0.33183602 -0.18652412 0.30117190 -0.06817728
## [2,] -0.22542294 -0.04083235 -0.35312644 0.04664983 0.09614775
## [3,] 0.01433649 -0.41192231 0.01097804 -0.18229417 0.05594686
## [4,] 0.30419338 -0.06031982 -0.40348745 -0.18685764 -0.25027768
## [5,] -0.01812861 0.35488767 0.18582667 -0.30078248 -0.12696582
## [6,] 1.00000000 -0.11542328 0.14555958 -0.20809003 -0.52294596
## [7,] -0.11542328 1.00000000 0.00134131 0.19818780 -0.07641849
## [8,] 0.14555958 0.00134131 1.00000000 -0.13079529 0.19107089
## [9,] -0.20809003 0.19818780 -0.13079529 1.00000000 0.04357162
## [10,] -0.52294596 -0.07641849 0.19107089 0.04357162 1.00000000eigenCX <- eigen(CX)
eigenCX$values## [1] 2.1098781 1.8613346 1.6634004 1.3622762 0.8860571 0.5938595 0.5433143
## [8] 0.3855146 0.3621043 0.2322609–> Vemos los valores propios en orden decreciente. –> La matriz P_muestra será la que en cada columna contenga los vectores propios asociados a los valores propios.
P_muestra <-eigenCX$vectors
P_muestra## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.50194048 -0.04810782 0.20797756 -0.2947649 -0.1933675 -0.19373313
## [2,] -0.29606381 -0.11814277 -0.33429361 0.4447779 0.1368221 0.69265364
## [3,] -0.13772258 0.04851697 0.49343799 0.3264523 -0.5679085 0.20590456
## [4,] -0.28438838 0.41455985 -0.30124135 0.2287619 0.1716341 -0.33081024
## [5,] 0.49352918 0.20490113 -0.07000146 0.2342814 -0.4403120 -0.06038837
## [6,] -0.04624206 0.58190398 0.13907908 -0.2948229 0.2242849 0.23483291
## [7,] 0.35347976 -0.05883026 -0.50405319 -0.2028558 -0.2020964 0.02329214
## [8,] 0.39010810 -0.05528759 0.40589307 -0.1908474 0.3684363 0.34503262
## [9,] -0.18359502 -0.36780421 -0.21246011 -0.5022657 -0.2701155 0.21466324
## [10,] 0.06331608 -0.53583559 0.15021459 0.2910856 0.3133390 -0.32619429
## [,7] [,8] [,9] [,10]
## [1,] -0.3492547 0.56562790 -0.21201825 0.23197278
## [2,] -0.1317839 0.17293389 -0.09738442 0.18320128
## [3,] -0.2049701 -0.13406620 0.36599328 -0.26414417
## [4,] -0.4941488 -0.24802193 -0.14398210 -0.37264279
## [5,] -0.2423243 -0.06592099 -0.44721333 0.44053224
## [6,] -0.1646255 -0.19631651 0.38283357 0.48511073
## [7,] -0.3404943 0.37127836 0.50641796 -0.15729695
## [8,] -0.4233681 0.10661549 -0.28821244 -0.33929384
## [9,] -0.2115255 -0.58865348 -0.15229494 0.01050851
## [10,] -0.3856369 -0.17904515 0.28566063 0.36690186A continuación se muestra los valores y vectores propios de la matriz con 5 variables y las 328 observaciones.
4. [10 %] Calcular directamente el análisis de componentes principales con la instrucción prcomp,llamando a la estructura resultante pca:
pca <- prcomp (X , scale = TRUE)pca$rotation es la matriz P que tiene por columnas los vectores propios de la matriz de covarianzas, ordenados según la magnitud de los valores propios. ¿Cu ál es la dimensión de pca.rotation que obtenéis y cuál es la dimensión teórica que debería tener esta matriz? ¿Qué está pasando?
–> La matriz pasa de 328x58008 a 58008x328. –> Se está haciendo una “rotación varimax” de la matriz original
“Rotación Varimax: Método de rotación ortogonal que minimiza el número de variables que tienen cargas altas en cada factor. Simplifica la interpretación de los factores. https://www.ibm.com/support/knowledgecenter/es/SSLVMB_sub/statistics_mainhelp_ddita/spss/base/idh_fact_rot.html”
Con esto se consigue que el dataset se reduzca de 58008 a 328 columnas. Donde: • Las columnas son los autovectores, asociados a las componentes principales. • Las filas, representan a las columnas.
Como las componentes principales son una combinación lineal de todas las variables, hay 58008 filas. Los valores de esta matríz representa el peso de las variables en la construcción de las componentes principales. Por ejemplo, A[1,1], representa la relación entre la primera variable con el primer componente principal.
dim(pca$rotation)## [1] 58008 3285. [10 %] Considerar vuestra fecha de nacimiento en la forma aammdd, por ejemplo 760827. Calcular, con R lo siguiente 81 + aammdd % % 10 y llamar al resultado S. ¿Cuántas componentes principales N son necesarias para retener un porcentaje de variabilidad igual o superior al S %? ¿Qué porcentaje de reducción representa trabajar con S componentes principales en lugar de todos los datos?
S <- 81 + 760930%%10
S## [1] 81(58008)/58008*100## [1] 100#summary(pca)Según función #summary(pma) con las 47 componentes principales se retiene el 81% de la variabilidad total del DataSet original. De 58008 variables hemos pasado a 47.
47/58008*100## [1] 0.08102331(58008-47)/58008*100## [1] 99.91898El porcentaje de reducción es de: 99.919%
6. [10 %] Representar gráficamente — en un único gráfico — los 328 puntos formados por la proyección de los datos originales normalizados sobre la primera y la segunda componentes principales. Utilizar un color diferente para cada uno de los cuatro estados estructurales (saludable, daño en el nivel 1, daño en el nivel 2 y daño en el nivel 3). ¿Observáis algo? ¿En qué disposición geométrica se encuentran los puntos?
Se crea la variable objetivo. Su naturaleza es categórica, por eso, se pasa a factor. Es esencial este paso para colorear los puntos.
estado <- seq(1, nrow(X))
estado[1:82] <- "Saludable"
estado[83:164] <- "Uno"
estado[165:246] <- "Dos"
estado[247:328] <- "Tres"
estado## [1] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [7] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [13] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [19] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [25] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [31] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [37] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [43] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [49] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [55] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [61] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [67] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [73] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" "Saludable"
## [79] "Saludable" "Saludable" "Saludable" "Saludable" "Uno" "Uno"
## [85] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [91] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [97] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [103] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [109] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [115] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [121] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [127] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [133] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [139] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [145] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [151] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [157] "Uno" "Uno" "Uno" "Uno" "Uno" "Uno"
## [163] "Uno" "Uno" "Dos" "Dos" "Dos" "Dos"
## [169] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [175] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [181] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [187] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [193] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [199] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [205] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [211] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [217] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [223] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [229] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [235] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [241] "Dos" "Dos" "Dos" "Dos" "Dos" "Dos"
## [247] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [253] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [259] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [265] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [271] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [277] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [283] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [289] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [295] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [301] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [307] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [313] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [319] "Tres" "Tres" "Tres" "Tres" "Tres" "Tres"
## [325] "Tres" "Tres" "Tres" "Tres"Primera vs Segunda componente.
Se observa que en el espacio de las primeras dos componentes, se logra separar bastante bien las clases. Aunque “saludable”" y el daño del tipo tres, son similares.
La posición geométrica que forman los puntos son elipsoides.
plot(
pca$x[,1],
pca$x[,2],
xlab="Primera Componente",
ylab="Segunda Componente",
pch=16,
col = factor(estado)
)
legend(106,-85,names(table(estado)), palette(), cex=.6
)
title("PCA")
7. [10 %] Representar gráficamente - en un único gráfico — los 328 puntos formados por la proyección de los datos originales normalizados sobre la primera y la tercera componentes principales. Utilizar un color diferente para cada uno de los cuatro estados estructurales (saludable, daño en el nivel 1, daño en el nivel 2 y da daño en el nivel 3). ¿Observáis algo? ¿En qué disposición geométrica se encuentran los puntos?
Ahora hay más ruido a la hora de separar los datos. Formas más lineales.
plot(
pca$x[,1],
pca$x[,3],
xlab="Primera Componente",
ylab="Tercera Componente",
pch=16,
col = factor(estado)
)
legend(106,-85,names(table(estado)), palette(), cex=.6
)
title("PCA")
8. [10 %] Nos han enviado nuevos datos (provenientes de 10 nuevas estructuras) y tenemos que decidir si cada uno de los datos pertenece a una estructura saludable o si tiene daño.
Para proyectar los datos sobre el modelo que hemos generado en la pregunta 4, tenéis primero que normalizar los datos respecto a los datos originales, es decir, a cada elemento de cada columna de la matriz dataN.Zs hay que restarle la media de todos los datos de la misma columna de la matriz data.Y. Hay que hacer lo mismo con la desviación típica, es decir, cada elemento de cada columna de la matriz dataN.Zs se ha de dividir por la desviación típica de todos los datos de la misma columna de la matriz data.Y. Llamar Z a la nueva matriz normalizada. Recordar que Z contiene 10 filas y 58008 columnas.
Carga de los nuevos datos.
dataN <- readMat("/Users/alexrodriguez/Desktop/tower_data_new.mat")
dataN <- dataN$Zs
str(dataN)## num [1:10, 1:58008] 3.18e-04 1.31e-05 3.00e-04 1.36e-04 2.46e-04 ...Escalado de datos con respecto a los datos originales.
index <- ncol(data$Y)
for (i in 1:index){
dataN[ ,i] = dataN[,i] - mean(data$Y[,i])
dataN[ ,i] = dataN[,i] / sd(data$Y[,i])
}
Z <- dataN
str(Z)## num [1:10, 1:58008] 0.517 -0.858 0.435 -0.305 0.19 ...9. [10 %] Proyectar los datos de la matriz Z sobre el espacio vectorial generado por las componentes principales mediante el producto matricial Z·pca$rotation. Llamar al resultado T. Comprobar que T es una matriz de 10 filas y 328 columnas.
Al muliplicar los valores de las observaciones de cada columna por los autovectores generados por los datos originales, se consigue general el espacio de las 328 componentes principales a los nuevos datos. Con ello se consigue contraer la matriz Z, a 10 filas (observaciones) y 328 columnas, los componentes principales.
T <- Z%*%pca$rotation
dim(T)## [1] 10 328- [10 %] Considerar vuestra fecha de nacimiento en la forma aammdd, por ejemplo 760827. Calcular, con R, lo siguiente 1 + aammdd % % 9 y llamar al resultado U.
(U <- 1 + 760930 %% 9)## [1] 8Recuperar el gráfico bidimensional de la pregunta 6 (primera componente principal versus segunda componente principal) y añadir, uno a uno, los dos puntos de las filas U y U+1. Para cada punto, visualmente (por ejemplo, según un criterio de proximidad), intentar clasificar, si podéis, las dos estructuras.
Se va a crear una nueva matriz con los puntos originales y el nuevo punto, que quedará en la última fila. Al factor estado, se le añade una nuevo valor, que quedará en la última posición del mismo, y que se corresponderá con la última fila de la matriz creada, que es el nuevo punto.
# Solo se cogen las 2 primeras componentes, que son las necesarias para graficar.
nuevos_datos <- rbind(as.data.frame(pca$x), T[U,])
estado <- as.character(estado) # El problema que al añadir un nuevo dato a un vector factor, lo convierte a numerico.
estado <- c(estado, "New") # Añadir nuevo valor.El nuevo valor quedaría clasificado como “DAÑO DOS”.
plot(
nuevos_datos[,1],
nuevos_datos[,2],
xlab="Primera Componente",
ylab="Segunda Componente",
pch=16,
col = factor(estado, levels= c("Saludable", "Uno", "Dos", "Tres", "New"))
)
legend(106,-85, levels(factor(estado, levels= c("Saludable", "Uno", "Dos", "Tres", "New"))), palette(), cex=.6)
title("PCA")
Se añade la observación U+1.
nuevos_datos <- rbind(nuevos_datos, T[U+1,1:2])
estado <- as.character(estado) # El problema que al añadir un nuevo dato a un vector factor, lo convierte a numerico.
estado <- c(estado, "New") # Añadir nuevo valor.# Por tanto, hay que volver a convertir en factor.Se grafica. Este nuevo punto es difícil de clasificar. Estaría entre SALUDABLE y DAÑO TRES. Posiblemente sea saludable.
plot(
nuevos_datos[,1],
nuevos_datos[,2],
xlab="Primera Componente",
ylab="Segunda Componente",
pch=16,
col = factor(estado, levels= c("Saludable", "Uno", "Dos", "Tres", "New"))
)
legend(106,-85, levels(factor(estado, levels= c("Saludable", "Uno", "Dos", "Tres", "New"))), palette(), cex=.6)
title("PCA")
A la vista del problema, proponer un método sistemático para clasificar de forma automática un nuevo punto (pregunta abierta).
Es una labor teodiosa ir clasificando punto por punto. Para automatizar el proceso, se usa un modelo de clasificación que se beneficie del preprocesado realizado con componentes principales.
Creación del dataset.
data_pca <- data.frame(pca$x[,1:3]) # Devuelve matriz.
estado <- estado[-length(estado)] # Se quitan las dos nuevas filas.
estado <- estado[-length(estado)]
str(data_pca)## 'data.frame': 328 obs. of 3 variables:
## $ PC1: num 75.9 -95.2 92.1 -72.1 34.6 ...
## $ PC2: num -65.8 24.3 25.8 -69.6 92.9 ...
## $ PC3: num 0.222 26.102 -22.139 -12.569 18.609 ...str(estado)## chr [1:328] "Saludable" "Saludable" "Saludable" "Saludable" "Saludable" ...data_pca <- cbind(data_pca, estado)Se usa un algoritmo que usa distancias para clasificar los nuevos puntos. El parámetro a optimizar es el nº de vecinos. El que mejor resultado da es 5. Se hace el grid search con nº impares porque:
- KNN usa los vecinos más cercanos para clasificar, si tiene 3 vecinos saludables y dos de daño tres, le clasifica como saludable. Se se usan nº pares, puede haber empates.
knn_grid <- train(
estado ~ PC1 + PC2 + PC3,
data = data_pca,
method = "knn",
trControl = trainControl(
method = "repeatedcv",
number = 5,
repeats = 10,
classProbs = TRUE
),
tuneGrid = expand.grid(k=seq(3,50,2)),
metric = "Accuracy",
)
plot(knn_grid) # 5 vecinos
predict(knn_grid, T)## [1] Saludable Saludable Uno Dos Saludable Tres Saludable
## [8] Saludable Uno Dos
## Levels: Dos Saludable Tres Uno11. Para facilitar la corrección y la localización de los posibles errores de la programación de vuestro código en R, completar la siguiente tabla en relación a las variables creadas durante la resolución de esta práctica (para vuestro valor de S de la pregunta 5 y para vuestro valor de U de la pregunta 10):
| Estadístico | Valor |
|---|---|
| X[2,2] | 0.2551831 |
| dim(pca$rotation)[1] | 58008 |
| S | 81 |
| N | 47 |
| % Reducción | 99.91897669287 |
| T[U,1] | PC1 91.49974 |
| T[U+1,2] | PC2 -126.2022 |
| Clasificación T[U,] | Dos |
| Clasificación T[U+1,] | Saludable |