Ejercicio 1

Trabajaremos con la base de datos USArrests . Esta base de datos contiene el porcentaje de asaltos (Assault), asesinatos (Murder) y violaciones (Rape) por cada 100,000 habitantes para cada uno de los 50 estados de USA (1973). Además, también incluye el porcentaje de la población de cada estado que vive en zonas urbanas (UrbanPoP).

# Insertamos base de datos USArrests

data1 <- USArrests
head(data1)

Ahora relizaremos los cinco pasos para PCA…

Paso 1- Estandarizar de los datos

colMeans(data1)
##   Murder  Assault UrbanPop     Rape 
##    7.788  170.760   65.540   21.232

En el vector de medias podemos ver que asalto es el más alto entre las demás variables, seguido por UrbanPop, luego Rape y por ulto Murder. Esto nos dice las frecuencia de las variables, por ende, el asalto es la variables más frecuente entre las demás. Es decir, ocurre más a menudo.

Paso 2- Calcular la matriz de covarianza

Lo calcularemos de los datos reales y de los estandarizados.

# covarianza datos reales

cov(data1)
##              Murder   Assault   UrbanPop      Rape
## Murder    18.970465  291.0624   4.386204  22.99141
## Assault  291.062367 6945.1657 312.275102 519.26906
## UrbanPop   4.386204  312.2751 209.518776  55.76808
## Rape      22.991412  519.2691  55.768082  87.72916

Al Asalto ser una cantidad mucho mayor a las otras, se lleva un gran peso en nuestra prueba. Es por eso que es ideal estandarizar antes de sacar la covarianza.

# Covarianza datos escalados

cov(scale(data1))
##              Murder   Assault   UrbanPop      Rape
## Murder   1.00000000 0.8018733 0.06957262 0.5635788
## Assault  0.80187331 1.0000000 0.25887170 0.6652412
## UrbanPop 0.06957262 0.2588717 1.00000000 0.4113412
## Rape     0.56357883 0.6652412 0.41134124 1.0000000

Se debe hacer esto para que la media se cero y desviación estándar 1 antes de realizar el estudio PCA.

Al escalar, podemos ver unas contribuciones por variables más proporcionales. Además, podemos ver la correlación perfecta en la diagonal, cosa que la covarianza de los datos originales no tenía. Esta diagonal siempre debe ser 1 porque representa la relación de una variable consigo mismo, por ejemplo Murder con Murder.

Por otro lado vemos relaciones fuertes en Murder con Asalto, indicando que donde ocurrén Asesinatos tambien suelen ocurrir Asaltos. Esa es la relación más fuerte entre las variables, las demás tienen relación moderada o baja.

Paso 3- Calcular los componentes principales

Aquí se crean los componentes principales.

# Creación de componentes principales

pca <- prcomp(data1, scale. = TRUE)
round(pca$rotation,3)
##             PC1    PC2    PC3    PC4
## Murder   -0.536 -0.418  0.341  0.649
## Assault  -0.583 -0.188  0.268 -0.743
## UrbanPop -0.278  0.873  0.378  0.134
## Rape     -0.543  0.167 -0.818  0.089

A lo que me refiero con “aquí se crean los componentes principales” es que el PC1 está compuesto por…

PC1= −0.5359 Murder−0.5832 Assault−0.2782 UrbanPop−0.5434 Rape

Y así con los demás componentes.

Podemos ver como en el primer componente destacan los asaltos y asesinatos. Podriamos decir que ese primer componente se encarga sobre la información correspondiente a delitos. Por otro lado, en el PC2 resalta el porcentaje de población que vive en zonas urbanas, esto significando que este componente tiene mayormente la información relacionado con la población en zonas urbanas.

Paso 4- Selección de los componentes principales

Despues de calcular los componentes principales, los valores propios nos diran cuales de las varianza explica o conlleva mayor información entre todos los componentes. Mientra mayor sea el valor propio, mayor explicación nos da. Mientra menor sea el valor propio, menor información nos da y se podría eliminar sin perder información relevante.

# Valores propios 

prop_varianza <- pca$sdev^2/sum(pca$sdev^2)
prop_varianza
## [1] 0.62006039 0.24744129 0.08914080 0.04335752

Podemos ver como el primer componente principal explica el 62% de la varianza observada en los datos. Es decir, al proyectar los datos sobre esta componente se retiene la mayoría de la información. El segundo componente explica el 25% de la variabilidad o varianza observada. Con esos dos componentes ya explicamos o tenemos información de aproximádamente el 87% de la varianza observada.

Esto lo podemos ver por medio de la varianza acumulada.

# Varianza acumulada

prop_varianza_acum <- cumsum(prop_varianza)
round(prop_varianza_acum,2)
## [1] 0.62 0.87 0.96 1.00

Podemos ver lo que explicamos anteriórmente. El primer componente explica el 62% de la varianza observada. Este componente unido con el segundo explican el 87% de la varianza observada.

Podemos ver estos componente y sus valores en un gráfico.

# Varianza por cada componente principal 

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
ggplot(data = data.frame(prop_varianza_acum, pc= 1:4), aes(x=pc, y=prop_varianza_acum, group = 1)) + geom_point() + geom_line() + theme_bw() + labs(x= "Componente principal", y= "Prop. varianza explicada acumulada" )

Vemos como el primer componente empieza en el 62% que mencionamos en los análisis anteriores. Luego, acumulado con el segundo componente, alcanzamos el 87% de la explicación de la varianza observada. Continuamos con 96% y luego 100%.

Paso 5- Transformar los datos al nuevo espacio dimensional

Esta es una manera de ver las variables y como influyen con las componentes principales. Veremos unas flechas que apuntarán hacia una dirección.

# Nuevo espacio dimensional

pca$rotation <- -pca$rotation
pca$x        <- -pca$x
biplot(pca, scale = 0, cex = 0.5, col = c("blue2", "red2"))

Viendo las flechas, podemos destacar que UrbanPop apunta bastante hacia PC1 y está bastante alargado. Esto nos dice que la variable UrbanPop influye bastante en PC1.

Tambien podemos ver como Murder y Asauult van hacia el lado al igual que hacia arriba. Esto lo que nos explica es que Murder y Assault influyen tanto en PC1 como en PC2. Rape de igual modo, la diferencia es que va hacia abajo tal vez significando que influye más en PC1.

Respecto a las ciudades, si un estado está cerca de una flecha, significa que tiene un valor alto en la variable que representa esa flecha. Por ejemplo, Rape esta cercano de Texas, illinois, colorado y New York. Eso quiere decir que esos estados tienen valores significante en esa variable. Otro ejemplo a destacar es como Murder apunta directamente a Georgia, esto queriendo decir que Georgia tiene un valor muy elevado en la variables Murder.

En conclusión, debemos fijarnos en la dirección de laas flechas. Si apunta hacia PC1, PC2 o una combinación de ambas. Además, debemos observar que tan larga es la flecha. Mientras más larga, más influye. Los estados agrupados tienen patrones similares en el conjunto de dato y la flecha que apunte a su dirección, es la variable con mayor relación a ese estado o grupo de estados.

Ejercicio 2

Trabajaremos con una base de datos bastante conocida denominada iris . Este famoso conjunto de datos (de Fisher) proporciona las medidas en centímetros de las variables longitud y ancho del sépalo y longitud y ancho de los pétalos, respectivamente, para 50 flores de cada una de las 3 especies de iris (setosa, versicolor y virginica).

# Insertamos base de datos iris

data2 <- iris
head(data2)

En este ejercicio se siguieron los pasos de PCA pero no todos, solo los relevante. En primer lugar se relizaron diversos grafico ya que al ser solo 3 variables, se pueden hacer gráficos 3D bastante interesante.

Diferentes visualizaciones dinámicas de los datos

# Dispersión entre variables

library(plotly)
## Warning: package 'plotly' was built under R version 4.3.3
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
colores <- c('red2','blue2','green2')
axis = list(showline=F,zeroline=F,gridcolor='#ffff',ticklen=4,titlefont=list(size=13))

plot_ly(data2) %>%
  add_trace(type = 'splom',dimensions = list(list(label='sepal length', values=~Sepal.Length),
  list(label='sepal width', values=~Sepal.Width),list(label='petal length', values=~Petal.Length),list(label='petal width', values=~Petal.Width)),
  color = ~Species, colors = colores,marker = list(size = 7,line = list(width = 1,color ='rgb(230,230,230)' ))) %>% 
  style(diagonal = list(visible = F)) %>%
  layout(hovermode='closest',dragmode= 'select',plot_bgcolor='rgba(240,240,240, 0.95)',xaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),yaxis=list(domain=NULL, showline=F, zeroline=F, gridcolor='#ffff', ticklen=4),xaxis2=axis,xaxis3=axis,xaxis4=axis,yaxis2=axis,yaxis3=axis,yaxis4=axis)

Si no hay dispersión significa que la variabilidad de la variables no es mucha, es decir que el promiedio de los datos es bastan parecido. Por otro lado, si la dispersión es mayor quiere decir que existe mucha variabilidad en el grupo. ´

# Comportamiento 3D

plot_ly(data2, x = ~Sepal.Length, y = ~Petal.Length, z = ~Sepal.Width, color = ~Species, colors = colores) %>% add_markers()

Para hacer un grafico 3D se debe tener x,y y z. Si los puntos se agrupan en ciertas regiones del espacio 3D, eso sugiere una relación entre las variables. Si los puntos están dispersos en todas direcciones, es probable que no haya relación clara entre ellas.

Paso 2- Calcular la matriz de covarianza

# Matriz de covarianza datos originales

round(cov(data2[,-5]),2)
##              Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length         0.69       -0.04         1.27        0.52
## Sepal.Width         -0.04        0.19        -0.33       -0.12
## Petal.Length         1.27       -0.33         3.12        1.30
## Petal.Width          0.52       -0.12         1.30        0.58
# Visualización del comportamiento de cada variable en sus niveles

p1 <- plot_ly(data2, x = ~Sepal.Length,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(barmode = "stack", yaxis = list(title = 'Sepal length',range = c(0,40)))
p2 <- plot_ly(data2, x = ~Sepal.Width,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Sepal width',range = c(0,40)))
p3 <- plot_ly(data2, x = ~Petal.Length,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Petal length',range = c(0,50)))
p4 <- plot_ly(data2, x = ~Petal.Width,color = ~Species, colors= colores, type = 'histogram') %>% 
layout(yaxis = list(title = 'Petal width',range =c(0,50)))

subplot(p1,style(p2, showlegend = F),style(p3, showlegend = F),style(p4, showlegend = F),nrows = 2,titleY = T,shareX=F,margin = 0.1)

Paso 3- Calcular los componentes principales

# Calcular componentes principales 

library(stats)
iris.pca <- prcomp(data2[-5], scale=T)
iris.pca$rotation
##                     PC1         PC2        PC3        PC4
## Sepal.Length  0.5210659 -0.37741762  0.7195664  0.2612863
## Sepal.Width  -0.2693474 -0.92329566 -0.2443818 -0.1235096
## Petal.Length  0.5804131 -0.02449161 -0.1421264 -0.8014492
## Petal.Width   0.5648565 -0.06694199 -0.6342727  0.5235971

Paso 4- Selección de los componentes principales

Este código hace ambas la proporción de varianza y la acumulación de proporción a la vez

# Proporción y acumulación 

summary(iris.pca)
## Importance of components:
##                           PC1    PC2     PC3     PC4
## Standard deviation     1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

Veamos la gráfica.

# Gráfica de CODO

prop_varianza_acum <- summary(iris.pca)$importance[3,]*100

ggplot(data = data.frame(prop_varianza_acum, pc = 1:4),
       aes(x = pc, y = prop_varianza_acum, group = 1)) +
  geom_point() +
  geom_line() +
  geom_label(aes(label = round(prop_varianza_acum,2))) +
  theme_bw() +
  labs(x = "Componente principal",y = "Prop. varianza explicada acumulada")

Podemos de igualmodo graficar las componentes de manera 3D.

#PCA 3D

pca_3          <- prcomp(subset(data2, select = -c(Species)), rank = 3)
components     <- data.frame(pca_3[["x"]])
components     <- cbind(components, data2$Species)

tit = 'Total Varianza explicada = 99.48'

plot_ly(components, x = ~PC1, y = ~PC2, z = ~PC3, color = ~data2$Species, colors = colores ) %>%
  add_markers()

Paso 5- Transformar los datos al nuevo espacio dimensional

# Biplot
pca_2              <- prcomp(subset(data2, select = -c(Species)), rank = 2)
components         <- data.frame(pca_2[["x"]])
exp_variance       <- summary(pca_2)$importance[1,][1:2]
comp               <- -pca_2$rotation
loadings           <- -comp
for (i in seq(exp_variance)){ loadings[,i] <- comp[,i] * exp_variance[i]}

features = c('sepal_length', 'sepal_width', 'petal_length', 'petal_width')

fig <- plot_ly(components, x = ~PC1, y = ~PC2, color = ~data2$Species, colors = colores, type = 'scatter', mode = 'markers') %>%
  layout(legend=list(title=list(text='color')), plot_bgcolor = "#e5ecf6",
    xaxis = list(title = "0"),yaxis = list(title = "1"))
for (i in seq(4)){
fig <- fig %>%
    add_segments(x = 0, xend = loadings[i, 1], y = 0, yend = loadings[i, 2], line = list(color = 'purple'),inherit = FALSE, showlegend = FALSE) %>%
    add_annotations(x=loadings[i, 1], y=loadings[i, 2], ax = 0, ay = 0,text = features[i], xanchor = 'center', yanchor= 'bottom')}

fig