Practica 4: Analisis de componentes principales

Carga de librerias

# ------------------------------------------------------------------------------
# Instalar devtools:
# ------------------------------------------------------------------------------
# Antes es necesario instalar la libreria libgit2 del sistema: 
# - MacOS: brew install libgit2
# - Arch: yay -S libgit2
#
# Luego instalamos en R las sigueintes librerias:
# install.packages("gert")
# install.packages("devtools")
# ------------------------------------------------------------------------------
#
#
#
# ------------------------------------------------------------------------------
# Instalar pacman: Gestor de paquetes para R (https://github.com/trinker/pacman)
# ------------------------------------------------------------------------------
# install.packages('pacman')
# ------------------------------------------------------------------------------
#
#
#
# Cargar librería requeridas (Se instalan si no lo estan):
library(pacman)
p_load(
  tidyverse,
  ggrepel,
  devtools,
  ggbiplot,
  nortest,
  readxl,
  WVPlots,
  plotly
)

Preparacion de datos

Lectura de base de datos:

datos_crudos <- read_xlsx(
  "rendimiento academico.xlsx",
  sheet = 1,
  col_types = "numeric"
)

Nota: Se fuerza que todas las columnas sean numéricas.

Seleccionamos el grupo de alumnos a analizar según su situación académica:

situacion_analisis <- c(1, 2, 3)

En este caso vamos a seleccionar todos los grupos, recordemos que:

1: Es Cursante.
2: Desertor.
3: Egresado.

Cuatrimestre inicial de análisis (siempre hasta el presente):

cuat_ini_analisis <- 19900 # Año 1990 y cuatrimestre 0.

Variables del dataset a analizar(Nro. de columna):

columnas <- c(1, 2, 7, 9, 12, 21)

print(paste(
  'Columnas:', 
  paste(names(datos_crudos)[columnas], collapse=', '), 
  '.', 
  sep=''
))

## [1] "Columnas:legajo, carrera, sexo, sit_ac, tot_cuats, prom_alc."

Filtramos los registros que sean posteriores a cuat_ini_analisis y también aquellos que tengas las situaciones académicas que definimos mas arriba (En este caso todas):

datos_numericos_1 <- datos_crudos[
  ( datos_crudos$acuat_ing >= cuat_ini_analisis )&
  datos_crudos$sit_ac %in% situacion_analisis,
]

datos_numericos_2 <- data.frame(datos_numericos_1[columnas])

Seleccionas las columnas a que nos interesan para este análisis:

datos_numericos_2 <- data.frame(datos_numericos_1[columnas])

Nos quedamos con las filas que no tengan valores nulos, es decir, aquellas que este completas:

datos_completos <- datos_numericos_2[complete.cases(datos_numericos_2),]
names(datos_completos)

## [1] "legajo"    "carrera"   "sexo"      "sit_ac"    "tot_cuats" "prom_alc"

complete.cases: Return a logical vector indicating which cases are complete, i.e., have no missing values.

Quitamos la columna legajo:

datos_para_agrupar <- datos_completos[c(2:ncol(datos_completos))]
names(datos_para_agrupar)

## [1] "carrera"   "sexo"      "sit_ac"    "tot_cuats" "prom_alc"

Análisis exploratorio

Distribuciones de cada variable

par(mfcol = c(2,3))
# par(mfcol = c(1,length(datos)))

for ( k in 1:length(datos_para_agrupar) ) {
  hist(
    datos_para_agrupar[, k],
    proba = T,
    main = names(datos_para_agrupar[k]),
    10,
    xlab = ""
  )
  
  x0 <- seq(
    min(datos_para_agrupar[, k]), 
    max(datos_para_agrupar[, k]), 
    le = 50
  )
  
  lines(
    x0, 
    dnorm(x0, mean(datos_para_agrupar[, k]), sd(datos_para_agrupar[, k])), 
    col = "red", 
    lwd = 2
  )

  grid()
}

Agrupamos por carrera y calculamos promedio para cada variable:

datos_agrupados <- datos_para_agrupar %>%
  group_by(carrera) %>%
  dplyr::summarise(
      cantidad = n(),
      prop_sexo2      = sum(sexo==2)/(sum(sexo==1)+sum(sexo==2)),
      prop_egresados  = sum(sit_ac==3)/(sum(sit_ac==3)+sum(sit_ac==2)),
      prom_egresados  = mean(prom_alc[sit_ac==3]),
      prom_tot_cuatr  = mean(tot_cuats[sit_ac==3])
  )
datos_agrupados

## # A tibble: 12 x 6
##    carrera cantidad prop_sexo2 prop_egresados prom_egresados prom_tot_cuatr
##      <dbl>    <int>      <dbl>          <dbl>          <dbl>          <dbl>
##  1      40     3326     0.216        0.479              6.95           15.8
##  2      41     8168     0.202        0.522              6.56           16.7
##  3      42      302     0.123        0.0179             7.68           16.8
##  4      43      162     0.247        0.0127             6.91           23  
##  5      44     1798     0.0595       0.440              6.94           18.0
##  6      45      490     0.0714       0.419              7.06           17.3
##  7      46     3644     0.0703       0.424              7.26           18.8
##  8      47     3151     0.484        0.460              7.11           15.3
##  9      48     3417     0.260        0.000327           7.35           30  
## 10      49     3556     0.128        0.310              7.10           18.6
## 11      51      176     0.273        0.478              7.32           14.4
## 12      53       48     0.229        0                NaN             NaN

Seleccionamos las columnas 2 a 6 columnas y quitamos las carreras 48 y 53:

datos_agrupados <- datos_agrupados[-c(12,9), 1:6]
datos           <- datos_agrupados[2:6]

Nota:

Se excluye la carrera 48: entiendo que no tiene practicante egresados ya que debe ser muy nueva.
La carrera 53 no tiene prom_egresados ni prom_total_cuatr (Valores Nulos).

Estandarizo datos:

datos_estandarizados <- data.frame(scale(datos))

Comparativa boxplot de las variables escaladas

p <- datos_estandarizados %>%
      pivot_longer(
        ., 
        cols = names(datos_estandarizados), 
        names_to  = "Variables", 
        values_to = "Frecuencia"
      ) %>%
      ggplot(aes(x = Variables, y = Frecuencia, fill = Variables)) +
      geom_boxplot(width=0.1) + 
      geom_jitter(shape=16, position=position_jitter(0.2))

ggplotly(p)

Dispersograma segmentado por carrera

datos_agrupados_est <- datos_agrupados %>%
  mutate_at(2:6, scale)  %>%
  mutate_at(1, as.character)

p <- PairPlot(
  datos_agrupados_est[1:6], 
  colnames(datos_agrupados_est)[2:6],
  " ",
  group_var = "carrera", 
  palette=NULL
) + 
scale_color_manual(values=unique(datos_agrupados_est$carrera))

ggplotly(p)

Matriz de correlaciones aplicada directamente a “datos”

matriz_de_correlaciones <- cor(datos)

round(matriz_de_correlaciones, 2)

##                cantidad prop_sexo2 prop_egresados prom_egresados prom_tot_cuatr
## cantidad           1.00       0.08           0.53          -0.62          -0.16
## prop_sexo2         0.08       1.00           0.12          -0.10          -0.31
## prop_egresados     0.53       0.12           1.00          -0.43          -0.62
## prom_egresados    -0.62      -0.10          -0.43           1.00          -0.22
## prom_tot_cuatr    -0.16      -0.31          -0.62          -0.22           1.00

Calcula los auto-valores desde la matrix de correlacion:

desc_mat_cor    <- eigen(matriz_de_correlaciones)
autovalores_cor <- desc_mat_cor$values

round(autovalores_cor, 2)

## [1] 2.21 1.41 0.90 0.37 0.10

¿Cuanta variabilidad concentra cada autovalor?

Lso auto-valores son iguales a las varianzas de las variables originales. Si calculamos la proporción de variabilidad de cada componente podemos compararlas y determina que componente captura mas variabilidad.

variabilidad_cor <- autovalores_cor / sum(autovalores_cor)

round(variabilidad_cor, 2)

## [1] 0.44 0.28 0.18 0.07 0.02

Nota: Cada auto-valor describe el grado de variabilidad de la componente a la que pertenece.

¿Es lo mismo calcular esta variabilidad con los auto-valores de una matriz u otra? Entiendo que no, cada matriz X va a tener sus auto-valores.

Componentes principales

Calculamos las componentes principales del dataset “datos”:

datos.pc <- prcomp(datos, scale = TRUE)

# datos.pc$sdev: Es la raiz cuadrada de los autovalores.
# Por esta cuestion elevamos al cuadrado para quedarnos con los auto-valores.
round(datos.pc$sdev^2, 2)

## [1] 2.21 1.41 0.90 0.37 0.10

Autovectores (Son las columnas):

round(datos.pc$rotation, 2)

##                  PC1   PC2   PC3   PC4   PC5
## cantidad       -0.54  0.30 -0.04  0.78 -0.11
## prop_sexo2     -0.22 -0.36  0.89  0.01 -0.18
## prop_egresados -0.59 -0.16 -0.32 -0.45 -0.57
## prom_egresados  0.44 -0.55 -0.22  0.43 -0.52
## prom_tot_cuatr  0.35  0.68  0.23 -0.09 -0.60

Nota: Las compontentes de cada auto-vector(columna de la tabla anterior) son los famosos loadings.

Vector de medias (De casualidad coinciden):

round(datos.pc$center ,2)

##       cantidad     prop_sexo2 prop_egresados prom_egresados prom_tot_cuatr 
##        2477.30           0.19           0.36           7.09          17.47

Vector de desvios:

round(datos.pc$scale, 2)

##       cantidad     prop_sexo2 prop_egresados prom_egresados prom_tot_cuatr 
##        2491.82           0.13           0.19           0.30           2.40

Loadings:

carga1 <- data.frame(
  cbind(
    X=1:length(datos),
    primeracarga=data.frame(datos.pc$rotation)[,1]
  )
)
carga2 <- data.frame(
  cbind(
    X=1:length(datos),
    segundacarga=data.frame(datos.pc$rotation)[,2]
  )
)
round(cbind(carga1,carga2),2)

##   X primeracarga X segundacarga
## 1 1        -0.54 1         0.30
## 2 2        -0.22 2        -0.36
## 3 3        -0.59 3        -0.16
## 4 4         0.44 4        -0.55
## 5 5         0.35 5         0.68

p <- ggplot(carga1, aes(X, primeracarga), fill=tramo ) + 
  geom_bar(stat="identity", position="dodge", fill ="royalblue", width =0.5 ) + 
  xlab('Variables') + 
  ylab('Primera Componente')

ggplotly(p)

p <- ggplot(carga2, aes (X, segundacarga), fill=X) +
    geom_bar(stat="identity", position="dodge", fill="royalblue", width=0.5) +
    xlab('Variables') +
    ylab('Segunda Componente')

ggplotly(p)

Cambiando el alfa?:

p <- ggbiplot(datos.pc, obs.scale=0.75 ,var.scale=0.75, alpha=1)

ggplotly(p)

Calculo de scores de cada individuo:

CP1 = as.matrix(datos_estandarizados)%*%as.matrix(carga1[2])
CP2 = as.matrix(datos_estandarizados)%*%as.matrix(carga2[2])
datos_agrupados$CP1 = CP1
datos_agrupados$CP2 = CP2

p <- ggbiplot(
  datos.pc, 
  obs.scale=0.75,
  var.scale=0.75,
  alpha=1,
  labels=datos_agrupados$carrera
)

ggplotly(p)