R Notebook

polucion <- read.table("R/visualizacion/practica 1/datos_practica_1/data-usair.txt", header=TRUE)
rownames(polucion) <- polucion$City
polucion <- polucion %>% select(-c("SO2", "City"))

polucion 

Dias son precipitaciones al año Precip: precipitacon media anual en pulgadas

Wind: velocidad

Manuf: canitdad de empresas Pop: Population

resultado <- prcomp(polucion, center = TRUE, scale. = TRUE)
resultado 

Standard deviations (1, .., p=6):
[1] 1.4819456 1.2247218 1.1809526 0.8719099 0.3384829 0.1855998

Rotation (n x k) = (6 x 6):
               PC1        PC2         PC3         PC4         PC5         PC6
Temp   -0.32964613 -0.1275974 -0.67168611  0.30645728  0.55805638 -0.13618780
Manuf   0.61154243 -0.1680577 -0.27288633 -0.13684076 -0.10204211 -0.70297051
Pop     0.57782195 -0.2224533 -0.35037413 -0.07248126  0.07806551  0.69464131
Wind    0.35383877  0.1307915  0.29725334  0.86942583  0.11326688 -0.02452501
Precip -0.04080701  0.6228578 -0.50456294  0.17114826 -0.56818342  0.06062222
Days    0.23791593  0.7077653  0.09308852 -0.31130693  0.58000387 -0.02196062

Pareciera que en este caso necesito la tercer componente,

cumsum(resultado$sdev) / sum (resultado$sdev)

[1] 0.2804796 0.5122759 0.7357882 0.9008098 0.9648726 1.0000000

Osea necesito 3 components para el 70% y 40 para el 90%.

El ejercicio pretende que pensemos si hacer PCA sobre la matriz de covarianza o la de correlacion, en realidad al standarizar es lo mismo. El unico caso donde podemos no estandarizar es si las variables estan en la misma escala y nos importan las varianzas. Es decir, que mas varianza pueda afectar mas el resultado.

Componentes principales: miro su covarianza con las variables originales (recordemos componente distinto de direccion principal). Componente es una vez proyectado

La primer componente contrapone Temperatura con Manuf, Pop y Wind La segunda tiene una influencia positiva de las variables de precipitacion y negativa de todo lo otro. La tercera se ve influenciada negativamente por precipitacion y temperatura y positivamente por el viento.

Siempre recordemos esta intepretacion es siempre relativa al promedio de las variables.

library(factoextra)
fviz_pca_biplot(
  resultado,
 axes = c(1,2),
  repel    = TRUE,          # evita solapamiento de etiquetas
  legend.title = list(fill = "Observ./Var.")
) +
  ggtitle("Biplot") +
  theme_minimal()

library(factoextra)
fviz_pca_biplot(
  resultado,
 axes = c(2, 3),
  repel    = TRUE,          # evita solapamiento de etiquetas
  legend.title = list(fill = "Observ./Var.")
) +
  ggtitle("Biplot") +
  theme_minimal()

Reconstruccion de las variables originales

proyectar <- function(resultado, k, ciudad) {
  a <- t(resultado$rotation[,1:k])
  b <- t(as.matrix(polucion[ciudad,]))
  a%*% b  
}


pr <- proyectar(resultado, 2, "Buffalo")
reconstruir <- function(resultado, proyectado, k) {
  resultado$rotation[,1:k] %*% proyectado 
}

error_reconstruccion <- function(resultado, k, ciudad) {
  norm(reconstruir(resultado, proyectar(resultado, k, ciudad), k) - t( as.matrix(polucion[ciudad,])))  
}

error_reconstruccion(resultado, 2, "Buffalo")

[1] 738.3217

error_reconstruccion(resultado, 3, "Buffalo")

[1] 328.5152

error_reconstruccion(resultado, 4, "Buffalo")

[1] 243.8951

error_reconstruccion(resultado, 5, "Buffalo")

[1] 63.31033

Una extrema:

ciudad <- "Phoenix"
error_reconstruccion(resultado, 2,  ciudad)

[1] 858.4979

error_reconstruccion(resultado, 3,  ciudad)

[1] 500.6273

error_reconstruccion(resultado, 4,  ciudad)

[1] 443.1274

error_reconstruccion(resultado, 5,  ciudad)

[1] 401.1431

Una promedio

ciudad <- "Cincinnati"
error_reconstruccion(resultado, 2,  ciudad)

[1] 785.5298

error_reconstruccion(resultado, 3,  ciudad)

[1] 254.251

error_reconstruccion(resultado, 4,  ciudad)

[1] 150.3578

error_reconstruccion(resultado, 5,  ciudad)

[1] 29.79927

rownames(polucion) %>% purrr::map_vec(function(ciudad) error_reconstruccion(resultado, 2,  ciudad)) %>% var()

[1] 638187.9

rownames(polucion) %>% purrr::map_vec(function(ciudad) error_reconstruccion(resultado, 3,  ciudad)) %>% var()

[1] 53397.83

rownames(polucion) %>% purrr::map_vec(function(ciudad) error_reconstruccion(resultado, 4,  ciudad)) %>% var()

[1] 13752.1