Analitica de Datos, Maestria en Ciencias del Comportamiento, Universidad de San Andres, Primavera 2026
Lenguaje: R. Esta es la version R, equivalente y validada, de la notebook Python de la clase: mismo dataset, mismos pasos, mismos resultados. Usala si te sentis mas comodo/a en R; el dictado sigue la version Python.
Abrir en Colab
Objetivo. Introducir la reduccion de dimensionalidad, con foco en el analisis de componentes principales (PCA) y su utilidad para descubrir estructura en problemas de Ciencias del Comportamiento.
Al terminar esta clase vas a poder:
explicar por que muchas variables correlacionadas se pueden resumir en pocas dimensiones;
construir la intuicion de PCA en 2D y llevarla a 3 y mas dimensiones;
calcular una PCA con prcomp e interpretar sus loadings y scores en un biplot;
decidir cuando estandarizar y cuantos componentes retener (scree plot, varianza acumulada);
aplicar PCA (y combinarlo con clustering) para descubrir perfiles latentes de comportamiento.
Lectura obligatoria: James et al. (2023), An Introduction to Statistical Learning, capitulo 12 (secciones 12.1 y 12.2). La intuicion inicial retoma tambien la Seccion 6.3.1.
Tip: hace click en cualquier tema para saltar directo. Los ejercicios estan intercalados en la seccion que les corresponde.
Donde se ubica esta clase en la materia
Seguimos en el Eje III, aprendizaje no supervisado: datos sin variable a predecir, donde el objetivo es descubrir estructura. En la Clase 9 agrupamos observaciones con clustering; hoy reducimos variables, y al final combinamos las dos ideas.
Sobre que se apoya. La nocion de varianza y correlacion de Estadistica; las combinaciones lineales y proyecciones del algebra lineal; y el K-means de la Clase 9, que reaparece en el ejercicio final.
A que habilita. Representar datos de alta dimension en pocas coordenadas es la base de los embeddings de las Clases 11 y 12 (NLP), y una herramienta central para el proyecto final.
se apoya en
esta clase
habilita
Varianza y correlacion; proyecciones lineales; K-means (Clase 9)
PCA: de la intuicion 2D a datos de alta dimension
Embeddings en NLP (Clases 11 a 12); proyecto final
El hilo de hoy
Muchas variables academicas correlacionadas de cada estudiante \(\to\) mirarlas todas es imposible (hay \(p(p-1)/2\) scatterplots) \(\to\) en 2D, PCA encuentra la recta de maxima varianza\(\to\) la misma idea escala a 3 y mas dimensiones, y la leemos en un biplot\(\to\) decidimos estandarizar y cuantos componentes retener \(\to\) vemos que t-SNE y UMAP capturan estructura no lineal \(\to\) aplicamos todo a un dataset que ya conoces (attrition) para descubrir perfiles latentes.
A lo largo de la clase nos acompana un dataset real de desercion y exito academico: 4424 estudiantes de educacion superior. Vamos a ver que PCA encuentra, casi sola, un eje que separa a quienes desertan de quienes se graduan.
## 1. Motivacion: demasiadas variables
La idea. El aprendizaje no supervisado no tiene una variable a predecir: su objetivo es descubrir estructura en un conjunto de features, casi siempre como parte de un analisis exploratorio. Un sitio de compras online, por ejemplo, agrupa a sus clientes segun su historial de navegacion y compra para recomendarle a cada uno lo que le interesa, sin ninguna etiqueta que diga de antemano a que grupo pertenece cada persona.
PCA ataca una version de ese problema: cuando hay muchas variables correlacionadas, permite resumirlas con unas pocas variables representativas que conservan casi toda la variabilidad. El problema practico es doble. Primero, no se pueden mirar: con \(p\) variables hay \(p(p-1)/2\) diagramas de dispersion posibles, y con \(p=10\) ya son 45. Segundo, muchas variables son en parte ruido. Empezamos midiendo cuanta redundancia hay en un dataset real de desercion academica.
# Dataset que nos acompana toda la clase: desercion y exito academico (UCI, CC BY 4.0)est <-cargar_desercion()# Bloque de variables academicas numericas, con sentido sustantivo comunacademico <-c("Admission grade", "Age at enrollment","Curricular units 1st sem (enrolled)", "Curricular units 1st sem (approved)", "Curricular units 1st sem (grade)","Curricular units 2nd sem (enrolled)", "Curricular units 2nd sem (approved)", "Curricular units 2nd sem (grade)")etq <-c("Nota admision", "Edad ingreso", "Inscriptas 1S", "Aprobadas 1S", "Nota 1S","Inscriptas 2S", "Aprobadas 2S", "Nota 2S")names(etq) <- academicoX <- est[, academico]p <-ncol(X)n_scatter <- p * (p -1) %/%2cat("Estudiantes:", nrow(est), "| variables del bloque:", p, "| scatterplots posibles:", n_scatter, "\n")print(table(est$Target))Rc <-cor(X)# Verificacion: las materias aprobadas en 1er y 2do semestre estan muy correlacionadasstopifnot(Rc["Curricular units 1st sem (approved)", "Curricular units 2nd sem (approved)"] >0.85)Rl <-as.data.frame(as.table(Rc)); names(Rl) <-c("v1", "v2", "cor")Rl$v1 <-factor(etq[as.character(Rl$v1)], levels =unname(etq))Rl$v2 <-factor(etq[as.character(Rl$v2)], levels =unname(etq))g <-ggplot(Rl, aes(v1, v2, fill = cor)) +geom_tile() +geom_text(aes(label =sprintf("%.2f", cor)), size =2.3, colour = COL[["ink"]]) +scale_fill_gradient2(low ="#B4232E", mid ="white", high ="#00529B", midpoint =0, limits =c(-1, 1)) +labs(x =NULL, y =NULL, title ="Varias variables academicas se mueven juntas: hay redundancia") +theme(axis.text.x =element_text(angle =45, hjust =1, size =8), axis.text.y =element_text(size =8))print(g); save_r("clase10_correlaciones", g, w =7, h =6)callout(paste0("Con solo ", p, " variables ya hay ", n_scatter," scatterplots posibles, y ninguno muestra la foto completa. PCA busca, en cambio, ","unas pocas direcciones que resuman toda la nube. Empecemos con la version mas simple: dos variables."),tipo ="info", titulo ="Por que necesitamos reducir")
Con solo 8 variables ya hay 24 scatterplots posibles, y ninguno muestra la foto completa. PCA busca, en cambio, unas pocas direcciones que resuman toda la nube. Empecemos con la version mas simple: dos variables.
## 2. PCA en dos dimensiones
La idea. Antes de ir a muchas variables, la intuicion se ve entera en 2D. Tomamos dos variables muy correlacionadas: las materias aprobadas en el primer y el segundo semestre. Los puntos forman una nube alargada en diagonal.
La primera componente principal es la direccion a lo largo de la cual los puntos varian mas. Tiene una segunda lectura equivalente: es tambien la recta mas cercana a todos los puntos (la que minimiza las distancias). Proyectar cada estudiante sobre esa recta da su score: un unico numero que resume sus dos variables. Aca ese numero es, en la practica, “cuantas materias aprueba”: un indice de avance academico.
# Dos variables muy correlacionadas: aprobadas en 1er y 2do semestrepar_vars <-c("Curricular units 1st sem (approved)", "Curricular units 2nd sem (approved)")D <-as.matrix(X[, par_vars])Ds <-scale(D)pca2 <-prcomp(Ds, center =FALSE, scale. =FALSE)v1 <- pca2$rotation[, 1]; if (sum(v1) <0) v1 <--v1pve2 <- pca2$sdev^2/sum(pca2$sdev^2)stopifnot(pve2[1] >0.9) # con correlacion ~0.9, PC1 concentra casi toda la varianzadf2 <-data.frame(x = Ds[, 1], y = Ds[, 2], Target = est$Target)linea <-data.frame(x =c(-3, 3) * v1[1], y =c(-3, 3) * v1[2])g <-ggplot(df2, aes(x, y)) +geom_point(aes(colour = Target), size =0.8, alpha =0.4) +scale_colour_manual(values = COL_TARGET, name =NULL) +geom_line(data = linea, aes(x, y), colour = COL[["ink"]], linewidth =1) +coord_equal() +labs(x ="Aprobadas 1er sem (estand.)", y ="Aprobadas 2do sem (estand.)",title ="PC1 = la direccion de maxima varianza (y la recta mas cercana a los puntos)")print(g); save_r("clase10_pca2d", g, w =6.5, h =6)tiles(list(c("PC1", sprintf("%.1f%%", pve2[1] *100), "de la varianza en 1 eje"),c("PC2", sprintf("%.1f%%", pve2[2] *100), "lo que queda"),c("Correlacion", sprintf("%.2f", cor(D[, 1], D[, 2])), "entre las 2 variables")))
PC1
95.2%
de la varianza en 1 eje
PC2
4.8%
lo que queda
Correlacion
0.90
entre las 2 variables
# Por que "maxima varianza": proyectamos la nube sobre una direccion que rotaangs <-seq(0, 180, by =1)var_ang <-sapply(angs, function(a) { u <-c(cos(a * pi /180), sin(a * pi /180)); var(as.numeric(Ds %*% u))})best <- angs[which.max(var_ang)]g <-ggplot(data.frame(ang = angs, v = var_ang), aes(ang, v)) +geom_line(colour = COL[["primary"]], linewidth =0.8) +geom_vline(xintercept = best, colour = COL[["good"]], linetype ="dashed") +annotate("text", x = best, y =min(var_ang), hjust =-0.05, colour = COL[["good"]], size =3.2,label =paste0("maximo = PC1 (", round(best), " grados)")) +labs(x ="Angulo de la direccion (grados)", y ="Varianza proyectada",title ="La varianza proyectada se maximiza justo en la direccion de PC1")print(g)stopifnot(abs(var_ang[which.max(var_ang)] -max(pca2$sdev^2)) <0.05)
Las dos lecturas, en una sola animacion
Mientras la recta gira, mira dos cosas al mismo tiempo. Los segmentos rojos son los errores de proyeccion: la distancia de cada punto a la recta. Los puntos rojos apoyados sobre la recta son las proyecciones, y lo estirada que queda esa nube de proyecciones es la varianza proyectada.
Cuando la recta alcanza la direccion de la primera componente pasan las dos cosas a la vez: los segmentos llegan a su longitud total minima y, en ese mismo instante, las proyecciones quedan lo mas dispersas posible. No son dos propiedades distintas que coinciden por casualidad, son la misma condicion escrita de dos maneras. Es la equivalencia que vamos a reencontrar en la Seccion 3.1, cuando veamos que minimizar el error de reconstruccion y maximizar la varianza explicada son el mismo problema.
Predecir. Si dos variables estuvieran perfectamente correlacionadas, cuanta varianza explicaria PC1?
Ver respuesta
Casi toda. Dos variables perfectamente correlacionadas viven sobre una recta: una sola direccion (PC1) captura el 100% de su variacion y la segunda no aporta nada. PCA rinde justamente cuando hay correlacion.
## 3. De dos dimensiones a muchas
La idea. En 2D vimos una recta; en 3D, las dos primeras componentes definen el plano mas cercano a la nube (el que mejor la resume). El grafico interactivo de abajo muestra la nube de estudiantes en tres variables academicas y ese plano: giralo con el mouse para ver como se apoya sobre los datos.
Con \(p\) variables la idea es la misma, aplicada sucesivamente: cada componente es una combinacion lineal normalizada de las variables
donde los coeficientes \(\phi_{jm}\) son los loadings (cuanto pesa cada variable en la componente \(m\)) y los valores \(Z_m\) son los scores. Cada componente maximiza la varianza sujeta a ser ortogonal (no correlacionada) a las anteriores. Un biplot muestra todo junto: los puntos son los scores; las flechas, los loadings.
# Vista interactiva 3D: el plano de las dos primeras componentes (interpretacion de la Figura 12.2)tri <-c("Curricular units 1st sem (approved)", "Curricular units 2nd sem (approved)", "Curricular units 2nd sem (grade)")Tm <-scale(as.matrix(X[, tri]))set.seed(SEED)idx <-sample(nrow(Tm), 250)pc3 <-prcomp(Tm)nrm <- pc3$rotation[, 3] # la normal del plano es la tercera componente# El plano generado por PC1 y PC2, escrito como z = f(x, y) usando la normalgs <-seq(-3, 3, length.out =14)zmat <-outer(gs, gs, function(a, b) -(nrm[1] * a + nrm[2] * b) / nrm[3])fig <-plot_ly() |>add_surface(x = gs, y = gs, z =t(zmat), opacity =0.45, showscale =FALSE,colorscale =list(c(0, COL[["muted"]]), c(1, COL[["muted"]]))) |>add_markers(x = Tm[idx, 1], y = Tm[idx, 2], z = Tm[idx, 3],color = est$Target[idx], colors = COL_TARGET, marker =list(size =2.5)) |>layout(title ="Las 2 primeras componentes definen el plano mas cercano a la nube",scene =list(xaxis =list(title ="Aprobadas 1S"), yaxis =list(title ="Aprobadas 2S"),zaxis =list(title ="Nota 2S"),camera =list(eye =list(x =1.9, y =0.5, z =1.4))))fig
# Biplot sobre el bloque academico completo (8 variables), coloreado por el desenlaceXs <-scale(as.matrix(X))pca <-prcomp(Xs, center =FALSE, scale. =FALSE)scores <- pca$xcomp <-t(pca$rotation) # filas = componentes, para espejar la version Pythonpve <- pca$sdev^2/sum(pca$sdev^2)# Orientar PC1 para que crezca con el exito academico (graduados con score alto)grad <- est$Target =="Graduate"; drop <- est$Target =="Dropout"if (mean(scores[grad, 1]) <mean(scores[drop, 1])) { comp[1, ] <--comp[1, ]; scores[, 1] <--scores[, 1] }stopifnot(mean(scores[grad, 1]) >mean(scores[drop, 1]))tiles(list(c("PC1", sprintf("%.1f%%", pve[1] *100), "exito / riesgo academico"),c("PC2", sprintf("%.1f%%", pve[2] *100), "segundo eje"),c("PC1 + PC2", sprintf("%.1f%%", sum(pve[1:2]) *100), "resumen en 2D")))biplot_r(scores, t(comp), unname(etq[academico]),"Biplot: PC1 separa quienes desertan de quienes se graduan","clase10_biplot", grupos = est$Target)print(round(setNames(comp[1, ], unname(etq[academico])), 2))
PC1
54.2%
exito / riesgo academico
PC2
17.3%
segundo eje
PC1 + PC2
71.5%
resumen en 2D
Nota admision Edad ingreso Inscriptas 1S Aprobadas 1S Nota 1S
0.03 -0.03 0.39 0.45 0.37
Inscriptas 2S Aprobadas 2S Nota 2S
0.39 0.45 0.38
El resultado es el que anticipaba la lectura: sin usar en ningun momento la etiqueta de desercion, PC1 termina siendo un eje de riesgo academico. Los estudiantes con score alto (muchas materias aprobadas y buenas notas) son casi todos graduados; los de score bajo, desertores. PCA descubrio sola esa estructura. Es la misma logica con la que el libro interpreta, sobre USArrests, que su PC1 mide el nivel de criminalidad de cada estado.
Como se lee un biplot
Un biplot no es un grafico, son dos graficos superpuestos sobre los mismos ejes, y cada uno tiene su propia escala.
Los scores. Cada punto es un estudiante, ubicado por su valor en la primera componente (eje horizontal) y en la segunda (eje vertical). Hasta aca es un scatterplot comun, solo que las dos coordenadas ya no son variables medidas sino combinaciones lineales de todas las variables.
Los loadings. Cada flecha es una variable original, dibujada desde el origen hasta el punto formado por su peso en la primera componente y su peso en la segunda. Como los pesos viven en una escala mucho mas chica que los scores, el biplot los multiplica por un factor comun para que las flechas se vean; ese factor es el mismo para todas, asi que las comparaciones entre flechas siguen siendo validas.
De las flechas se leen tres cosas, y conviene mirarlas en este orden.
Direccion: que mezcla de componentes representa a esa variable. Una flecha que apunta casi horizontal es una variable que carga casi solo sobre la primera componente; una casi vertical carga sobre la segunda; una diagonal, sobre las dos. La direccion es lo que permite ponerle nombre sustantivo a cada eje: si todas las flechas de rendimiento apuntan hacia la derecha, el eje horizontal es un eje de rendimiento.
Angulo entre flechas: cuanto se parecen dos variables. Flechas que apuntan casi para el mismo lado corresponden a variables muy correlacionadas entre si; flechas casi perpendiculares, a variables practicamente sin correlacion; flechas opuestas, a variables correlacionadas negativamente. La razon es que el coseno del angulo entre dos flechas aproxima la correlacion entre las dos variables, y la aproximacion es buena en la medida en que las dos esten bien representadas en el plano.
Longitud: que tan bien queda representada esa variable en este plano. Una flecha larga indica una variable que las dos primeras componentes capturan bien; una flecha corta indica una variable que vive sobre todo en las componentes que no estamos mirando, y sobre ella el biplot no dice casi nada. Ojo con este criterio: en un biplot dibujado con los pesos crudos, como el de arriba, la longitud es solo una guia aproximada, porque la primera componente tiene mas varianza que la segunda y esa diferencia no esta en el dibujo. La version objetiva se calcula, y es lo que hacemos en la celda que sigue.
La cantidad que hace objetiva la lectura es la varianza de cada variable captada por las dos primeras componentes:
\(c_j\) es la varianza de la variable \(j\) que reconstruyen las dos primeras componentes, tambien llamada comunalidad en dos dimensiones;
\(\lambda_1\) y \(\lambda_2\) son las varianzas de la primera y la segunda componente;
\(\phi_{j1}\) y \(\phi_{j2}\) son los pesos (loadings) de la variable \(j\) en cada una de esas dos componentes.
Como las variables estan estandarizadas, cada una tiene varianza 1 y entonces \(c_j\) se lee directo como proporcion: un valor cercano a 1 significa que el plano del biplot representa casi toda esa variable, y un valor cercano a 0, que casi no la representa. En la misma linea, \(\sqrt{\lambda_1}\,\phi_{j1}\) es la correlacion entre la variable \(j\) y la primera componente.
# Lectura objetiva del biplot: direccion, angulo y longitud de cada flechaang <- (atan2(comp[2, ], comp[1, ]) *180/ pi) %%360# direccion de cada flechalargo2 <- comp[1, ]^2+ comp[2, ]^2# longitud al cuadrado, tal como se dibujalam <- pca$sdev[1:2]^2# varianza de PC1 y PC2var_cap <- lam[1] * comp[1, ]^2+ lam[2] * comp[2, ]^2# comunalidad 2Dtabla_biplot <-data.frame(Variable =unname(etq[academico]),LoadingPC1 =round(comp[1, ], 2), LoadingPC2 =round(comp[2, ], 2),Angulo =round(ang, 0), Largo2 =round(largo2, 3),VarianzaCaptada =round(var_cap, 2))print(tabla_biplot[order(-tabla_biplot$VarianzaCaptada), ], row.names =FALSE)angulo_entre <-function(a, b) { d <-abs(ang[match(a, academico)] - ang[match(b, academico)]) %%360min(d, 360- d)}juntas <-c("Curricular units 1st sem (enrolled)", "Curricular units 2nd sem (enrolled)")cruz <-c("Age at enrollment", "Curricular units 1st sem (approved)")peor <- academico[which.min(var_cap)]cat("\nFlechas casi paralelas:", etq[juntas[1]], "y", etq[juntas[2]],"-> angulo", round(angulo_entre(juntas[1], juntas[2]), 1),"grados, correlacion", round(Rc[juntas[1], juntas[2]], 2), "\n")cat("Flechas casi perpendiculares:", etq[cruz[1]], "y", etq[cruz[2]],"-> angulo", round(angulo_entre(cruz[1], cruz[2]), 1),"grados, correlacion", round(Rc[cruz[1], cruz[2]], 2), "\n")cat("Variable peor representada en el plano:", etq[peor], "-> varianza captada", round(min(var_cap), 2), "\n")cor_pc1 <-sqrt(lam[1]) * comp[1, ]cat("\nCorrelacion de cada variable con PC1:\n"); print(round(setNames(cor_pc1, unname(etq[academico])), 2))# Verificaciones de cada afirmacion del textostopifnot(angulo_entre(juntas[1], juntas[2]) <10, Rc[juntas[1], juntas[2]] >0.85)stopifnot(angulo_entre(cruz[1], cruz[2]) >70, angulo_entre(cruz[1], cruz[2]) <110, abs(Rc[cruz[1], cruz[2]]) <0.2)stopifnot(peor =="Admission grade", min(var_cap) <0.2)stopifnot(abs(largo2[match(peor, academico)] -min(largo2)) <1e-12)stopifnot(max(var_cap) <1.0)mejor <- academico[which.max(abs(cor_pc1))]stopifnot(abs(cor_pc1[match(mejor, academico)]) >0.9)callout(paste0("Las cuatro variables de materias (inscriptas y aprobadas en los dos semestres) tienen loadings de PC1 ","del mismo signo y de magnitud parecida: son las que definen el eje horizontal. La que mas se alinea ","con PC1 es ", etq[mejor], ", con una correlacion de ", sprintf("%.2f", cor_pc1[match(mejor, academico)]),". La nota de admision, en cambio, casi no aparece en el plano: su varianza captada es ",sprintf("%.2f", min(var_cap)), ", asi que sobre ella el biplot no autoriza ninguna lectura."),tipo ="info", titulo ="Lo que dice la tabla")
Variable LoadingPC1 LoadingPC2 Angulo Largo2 VarianzaCaptada
Inscriptas 1S 0.39 0.42 48 0.330 0.90
Aprobadas 1S 0.45 0.02 3 0.207 0.90
Aprobadas 2S 0.45 -0.06 352 0.205 0.88
Inscriptas 2S 0.39 0.38 44 0.301 0.87
Nota 2S 0.38 -0.37 316 0.279 0.81
Nota 1S 0.37 -0.34 317 0.259 0.77
Edad ingreso -0.03 0.59 93 0.344 0.48
Nota admision 0.03 -0.27 276 0.074 0.10
Flechas casi paralelas: Inscriptas 1S y Inscriptas 2S -> angulo 3.3 grados, correlacion 0.94
Flechas casi perpendiculares: Edad ingreso y Aprobadas 1S -> angulo 90.7 grados, correlacion -0.05
Variable peor representada en el plano: Nota admision -> varianza captada 0.1
Correlacion de cada variable con PC1:
Nota admision Edad ingreso Inscriptas 1S Aprobadas 1S Nota 1S
0.06 -0.07 0.81 0.95 0.78
Inscriptas 2S Aprobadas 2S Nota 2S
0.82 0.93 0.79
Lo que dice la tabla
Las cuatro variables de materias (inscriptas y aprobadas en los dos semestres) tienen loadings de PC1 del mismo signo y de magnitud parecida: son las que definen el eje horizontal. La que mas se alinea con PC1 es Aprobadas 1S, con una correlacion de 0.95. La nota de admision, en cambio, casi no aparece en el plano: su varianza captada es 0.10, asi que sobre ella el biplot no autoriza ninguna lectura.
Los tres criterios aplicados a nuestro biplot
Direccion. En la tabla, las seis variables de cursada (inscriptas, aprobadas y notas de los dos semestres) tienen loadings de PC1 positivos y de magnitud comparable, mientras que la edad de ingreso y la nota de admision tienen loadings de PC1 practicamente nulos y toda su presencia en PC2. Esa es la lectura clave: el eje horizontal esta hecho de rendimiento academico, con todas las variables de cursada tirando para el mismo lado, y no de las dos variables de entrada. Por eso tiene sentido llamar a PC1 eje de exito academico, y por eso los graduados quedan a la derecha y los desertores a la izquierda sin que la etiqueta haya entrado nunca en el calculo.
Angulo. Las inscriptas del primer y del segundo semestre salen con flechas casi superpuestas, y su correlacion, impresa arriba, es altisima: dicen casi lo mismo. Lo mismo pasa con el par de aprobadas y con el par de notas, cada uno agrupado por su cuenta. Los tres grupos aparecen abiertos entre si en angulos moderados: cuantas materias toma un estudiante, cuantas aprueba y con que nota son cosas relacionadas pero no identicas. En cambio, la edad de ingreso sale casi perpendicular a las aprobadas del primer semestre, y la correlacion entre esas dos variables es practicamente cero: el biplot esta diciendo que empezar mas grande no tiene, en estos datos, relacion lineal con aprobar mas materias. Es exactamente el tipo de lectura que un heatmap de correlaciones da variable por variable y que el biplot da de un vistazo.
Longitud. Aca conviene mirar la columna de varianza captada antes que el dibujo. La nota de admision es la variable con la flecha mas corta bajo los dos criterios (largo dibujado y varianza captada), y es la unica con una comunalidad claramente baja: las dos primeras componentes casi no la reconstruyen. Traducido a prudencia interpretativa, la posicion de esa flecha en el grafico es poco informativa, y no habria que sacar conclusiones sobre la nota de admision a partir de este plano. La edad de ingreso queda en una situacion intermedia, y las seis variables de cursada quedan bien representadas.
Regla practica. Antes de interpretar la posicion de una flecha, mira cuanta varianza de esa variable capturan las dos componentes que estas dibujando. Una flecha mal representada puede apuntar a cualquier lado y no significa nada.
Que es exactamente una componente principal
Antes de seguir conviene frenar y decir en palabras que objeto acabamos de calcular, porque de aca en adelante todo lo demas depende de esto.
La primera componente principal es una variable nueva: se construye tomando las \(p\) variables originales, multiplicando cada una por un peso propio y sumando todo. Eso es lo que en algebra se llama una combinacion lineal. Hay infinitas combinaciones lineales posibles, una por cada juego de pesos que se nos ocurra; PCA elige la unica cuya variable resultante tiene la mayor varianza posible. Dicho de otro modo: entre todas las formas de resumir las ocho variables academicas en un solo numero por estudiante, PCA elige la que mas separa a los estudiantes entre si.
De ahi salen los dos conceptos que hay que tener claros y que se confunden todo el tiempo:
Loadings (los pesos): cuanto aporta cada variable original a la componente. Son \(p\) numeros por componente, uno por variable, y no dependen del estudiante. Son las flechas del biplot.
Scores (los valores): el resultado de aplicar esa combinacion a cada observacion. Son \(n\) numeros por componente, uno por estudiante, y dicen donde se ubica cada unidad a lo largo de la nueva variable resumen. Son los puntos del biplot.
Un supuesto previo: cada variable entra centrada a media cero (y en la practica tambien estandarizada, como vamos a discutir en la seccion 4). El centrado es una conveniencia matematica, no cambia la interpretacion: mueve el origen al centro de la nube para que “maximizar varianza” y “pasar por el medio de los puntos” sean la misma cosa.
Por que los pesos suman uno (al cuadrado)
En la formula de arriba aparecio una restriccion que puede parecer un detalle tecnico y no lo es: \(\sum_j \phi_{jm}^2 = 1\). La razon es simple. Si no la impusieramos, maximizar varianza no tendria solucion: alcanzaria con agrandar todos los pesos (duplicarlos, multiplicarlos por mil) para que la varianza de la variable resultante creciera sin limite. La normalizacion obliga a PCA a buscar una direccion y no una magnitud: fija el largo del vector de pesos en uno y deja libre solo hacia donde apunta. Por eso en la seccion 2 el laboratorio interactivo rotaba una recta sin cambiarle el largo.
Como se calcula en la practica. Ese problema de maximizacion tiene una solucion cerrada muy conocida: los loadings salen de la descomposicion en autovectores y autovalores (eigen decomposition) de la matriz de covarianzas de los datos, o de la matriz de correlaciones si las variables fueron estandarizadas. Cada componente principal es un autovector de esa matriz, y la varianza que esa componente explica es su autovalor asociado. No hace falta desarrollar el algebra: alcanza con la idea, porque quien resuelve el sistema es R cuando llamamos a prcomp().
Las componentes siguientes. La segunda componente se construye con la misma logica, maximizando varianza, pero con una exigencia extra: tiene que ser ortogonal a la primera, lo que equivale a decir que sus scores estan no correlacionados con los de la primera. La tercera es ortogonal a las dos anteriores, y asi sucesivamente. Cada componente captura variabilidad que las anteriores no habian capturado, sin repetir informacion. Como maximo hay \(\min(n-1,\, p)\) componentes: en la practica, con muchas mas observaciones que variables, hay tantas componentes como variables.
El score de un estudiante, termino a termino
La cuenta que convierte loadings en scores es una sola linea:
\(z_{im}\) es el score del estudiante \(i\) en la componente \(m\): el numero que lo ubica sobre ese eje;
\(x_{ij}\) es el valor estandarizado de la variable \(j\) para ese estudiante (cuantos desvios por encima o por debajo del promedio esta);
\(\phi_{jm}\) es el loading de la variable \(j\) en la componente \(m\): el peso, igual para todos los estudiantes;
\(p\) es la cantidad de variables del bloque, ocho en nuestro caso.
En el biplot cada estudiante es un punto y esa cuenta ya esta hecha por debajo, invisible. Vale la pena abrirla una vez. La celda siguiente toma un estudiante concreto del dataset, despliega los ocho productos \(\phi_{j1}\, x_{ij}\) uno por uno, los suma y compara el resultado con el score que devuelve pca$x.
# El score de UN estudiante en PC1, abierto termino a terminoi_est <-which.max(scores[, 1]) # el estudiante con el score mas altophi1 <- comp[1, ] # loadings de PC1x_i <- Xs[i_est, ] # sus valores estandarizadosprod <- phi1 * x_i # aporte de cada variable al scoresuma <-sum(prod)tabla_score <-data.frame(Variable =unname(etq[academico]),Valor_estandarizado =round(x_i, 2),Loading_PC1 =round(phi1, 2),Aporte =round(prod, 3))tabla_score <- tabla_score[order(-abs(tabla_score$Aporte)), ]print(tabla_score, row.names =FALSE)cat("\nSuma de los aportes:", round(suma, 4), "\n")cat("Score que devuelve prcomp:", round(scores[i_est, 1], 4), "\n")# La cuenta a mano y la del algoritmo son la mismastopifnot(abs(suma - scores[i_est, 1]) <1e-8)pct <-round(100*sum(sort(abs(prod), decreasing =TRUE)[1:4]) /sum(abs(prod)))g <-ggplot(transform(tabla_score, Variable =factor(Variable, levels =rev(tabla_score$Variable))),aes(Aporte, Variable)) +geom_col(fill = COL[["good"]]) +labs(x ="Aporte al score de PC1 (loading x valor estandarizado)", y =NULL,title ="De donde sale el score de un solo estudiante")print(g); save_r("clase10_score_termino", g, w =7, h =4.5)cat("Las cuatro variables de materias explican el", pct, "% del score en valor absoluto\n")
Variable Valor_estandarizado Loading_PC1 Aporte
Inscriptas 2S 7.64 0.39 3.004
Inscriptas 1S 6.75 0.39 2.618
Aprobadas 2S 5.16 0.45 2.314
Aprobadas 1S 4.94 0.45 2.247
Nota 2S 0.83 0.38 0.317
Nota 1S 0.81 0.37 0.304
Nota admision 0.28 0.03 0.008
Edad ingreso -0.17 -0.03 0.006
Suma de los aportes: 10.818
Score que devuelve prcomp: 10.818
Las cuatro variables de materias explican el 94 % del score en valor absoluto
Leido asi, el score deja de ser una salida opaca del algoritmo. El estudiante elegido queda arriba en PC1 porque casi todas sus variables academicas estan varios desvios por encima del promedio y porque los loadings de PC1 las pesan a todas con el mismo signo: PC1 es, esencialmente, un promedio ponderado de rendimiento academico. Las variables que mas empujan su score son las que combinan un loading grande con un valor personal muy alejado de la media; una variable con loading grande pero valor tipico aporta poco, y una variable con valor extremo pero loading chico tampoco mueve la aguja. Se necesitan las dos cosas.
Eso tambien explica por que, sin haber usado nunca la etiqueta de desercion, los graduados terminan del lado alto del eje: estan alto en las mismas variables que PC1 pondera fuerte.
Hacer esta cuenta a mano una sola vez es lo que hace que la palabra score deje de ser abstracta: no es una coordenada magica, es una suma de ocho productos que cualquiera puede reproducir con una calculadora.
Caso aplicado 1. De 4500 adjetivos a cinco factores: el modelo Big Five
Desde principios del siglo XX los psicologos intentaron catalogar los terminos del lenguaje cotidiano que describen diferencias de personalidad. En 1936 Gordon Allport y Henry Odbert armaron una lista de unos 4500 terminos, un problema con un \(p\) enorme y sin ninguna estructura visible. Con analisis factorial y la naciente tecnologia de computo, sucesivas investigaciones fueron condensando esa lista primero a 16 dimensiones y finalmente a 5, dando origen al modelo Big Five u OCEAN: Apertura, Responsabilidad, Extraversion, Amabilidad y Neuroticismo (historia del Big Five).
Un ejemplo con PCA explicita: De Raad y Barelds (2008) hicieron un estudio psicolexico del lexico holandes con 2365 items descriptivos de rasgos y aplicaron Analisis de Componentes Principales, obteniendo ocho factores que incluian los cinco clasicos mas Virtud, Competencia y Hedonismo (De Raad y Barelds).
El punto es el mismo que estamos viendo, llevado a otra escala: cuando \(p\) es enorme, las direcciones de maxima varianza permiten pasar de miles de columnas a un punado de factores compuestos que ademas se pueden nombrar e interpretar. Nosotros lo hacemos con 8 variables academicas; ellos lo hicieron con miles de adjetivos.
Caso aplicado 2. PC1 como indice: el Nivel Socioeconomico
Cuando decimos que PC1 funciona como un indice compuesto no es una metafora didactica: asi se construye buena parte de los indices de Nivel Socioeconomico que se usan en investigacion social. Filmer y Pritchett (2001), en Demography, estimaron la relacion entre riqueza del hogar y matriculacion escolar en India construyendo un indice de riqueza a partir de indicadores de propiedad de activos del hogar, y usaron PCA para derivar los ponderadores, es decir los loadings. El indice resulto robusto y sus resultados a nivel estadual se correspondieron con datos independientes de producto per capita y de pobreza (Filmer y Pritchett, 2001).
En la practica, el indice se arma con activos del hogar (heladera, moto, bicicleta) y caracteristicas de la vivienda (tipo de piso, fuente de agua); despues los hogares se ordenan segun su valor en PC1 y se los agrupa en quintiles de NSE. La guia metodologica de referencia es Vyas y Kumaranayake (2006), en Health Policy and Planning, que revisa como se construyen estos indices, su validez y limitaciones, como elegir las variables y como pasar de los scores a grupos de NSE (Vyas y Kumaranayake, 2006).
Hay un detalle que conviene discutir: la mayoria de estos trabajos usa solo la primera componente como indice. Literatura posterior propone promediar varias componentes, con el argumento de que PC1 suele explicar apenas una parte de la varianza total (discusion sobre cuantas componentes usar). Es exactamente la pregunta que retomamos mas adelante cuando miremos el scree plot y la PVE acumulada.
Caso aplicado 3. Nombrar componentes: percepcion de lacteos locales
Un experimento de eleccion online realizado en Apulia, Italia, con 543 encuestados, investigo como perciben los consumidores los productos lacteos locales en tres aspectos: calidad, sustentabilidad y disponibilidad. Los autores aplicaron PCA sobre las respuestas al cuestionario y quedaron con cuatro componentes, que nombraron segun que items cargaban fuerte en cada uno: sensibilidad a los atributos de calidad, “lo local es mejor”, “lo local es sustentable”, y demanda de mayor disponibilidad. Luego relacionaron esas cuatro dimensiones con variables sociodemograficas de los encuestados (Bimbo y colegas, 2022).
Sirve como plantilla del flujo tipico en investigacion aplicada: se disena un cuestionario con muchos items correlacionados, se usa PCA para encontrar las dimensiones latentes de opinion, se nombra cada componente leyendo los loadings mas altos, y recien despues se estudia como esas dimensiones se relacionan con otras variables. El nombre de la componente no sale del algoritmo, lo pone quien analiza.
### 3.1 Reconstruir un dato a partir de sus componentes
Hasta aca usamos PCA en una sola direccion: partimos de las variables originales y llegamos a los scores. Vale la pena recorrer el camino de vuelta, porque es lo que da sentido a la palabra comprimir. Si los scores y los loadings guardan la informacion de la tabla original, entonces tenemos que poder reconstruir cada dato a partir de ellos.
\(x_{ij}\): el dato que queremos reconstruir, es decir el valor de la variable \(j\) para la observacion \(i\), ya centrado (en nuestro caso, estandarizado).
\(M\): cuantas componentes usamos para reconstruir. Si \(M\) es chico, estamos aproximando la tabla entera con muy pocas dimensiones; si \(M = p\), la reconstruccion es exacta.
\(z_{im}\): el score de la observacion \(i\) en la componente \(m\), o sea la coordenada de ese estudiante sobre esa direccion.
\(\phi_{jm}\): el loading de la variable \(j\) en la componente \(m\), o sea cuanto pesa esa variable en esa direccion.
En palabras simples. Para reconstruir un dato agarras el score de esa observacion en cada componente, lo multiplicas por el loading de esa variable en esa misma componente, y sumas esos productos sobre las \(M\) componentes que decidiste usar. Cuantas mas componentes sumes, mas cerca vas a quedar del valor real.
Por que funciona. La reconstruccion hace exactamente el camino inverso al que abrimos termino a termino en la seccion anterior, deshaciendo la combinacion lineal con los mismos loadings. Es una relacion de ida y vuelta: los loadings sirven para armar los scores y tambien para desarmarlos. Cuando usamos las \(p\) componentes no se pierde nada, porque PCA no es mas que una rotacion del sistema de coordenadas: los mismos puntos, mirados desde otros ejes.
Y no es una aproximacion cualquiera. Se puede demostrar (no lo vamos a hacer aca) que esta es la mejor aproximacion posible con \(M\) dimensiones: ninguna otra eleccion de \(M\) variables nuevas construidas como combinaciones lineales de las originales reconstruye la tabla con menos error cuadratico. Eso es lo que hace a PCA optimo para comprimir, y no simplemente una forma mas de resumir datos.
# Reconstruir UN dato concreto usando cada vez mas componentesi_rec <- i_estj <-match("Curricular units 2nd sem (approved)", academico)media_j <-mean(X[[j]]); desvio_j <-sd(X[[j]])real <- X[i_rec, j]filas <-lapply(c(1, 2, p), function(M) { aprox_std <-sum(scores[i_rec, 1:M] * comp[1:M, j]) # suma de score x loadingdata.frame(M = M, Estandarizado =round(aprox_std, 3),Escala_original =round(aprox_std * desvio_j + media_j, 2),Error =round(abs(aprox_std * desvio_j + media_j - real), 2))})tabla_rec <-do.call(rbind, filas)cat("Valor real de", etq[academico[j]], "para ese estudiante:", real, "\n\n")print(tabla_rec, row.names =FALSE)# Con todas las componentes la reconstruccion es exacta, y el error final es el minimoerr <- tabla_rec$Errorstopifnot(err[length(err)] <1e-8, err[length(err)] <=min(err) +1e-12)
Valor real de Aprobadas 2S para ese estudiante: 20
M Estandarizado Escala_original Error
1 4.850 19.06 0.94
2 4.544 18.14 1.86
8 5.163 20.00 0.00
Vale la pena leer la tabla fila por fila. Con \(M = 1\) la reconstruccion ya cae bastante cerca del valor real, pero queda error: la primera componente resume la tendencia general de exito academico, y ubica a este estudiante en el lugar que le corresponde sobre ese eje, no en su valor exacto. A medida que sumamos componentes vamos recuperando el matiz propio de esa observacion, aquello en lo que se aparta de la tendencia comun, y el error tiende a achicarse. Con las ocho componentes el error es cero hasta la precision de la maquina: como no descartamos ninguna direccion, no hay nada que perder, solo cambiamos el sistema de coordenadas.
# Error de reconstruccion de TODO el dataset en funcion de cuantas componentes usemosMs <-1:pecm <-sapply(Ms, function(M) { aprox <- scores[, 1:M, drop =FALSE] %*% comp[1:M, , drop =FALSE]mean((Xs - aprox)^2)})# El error que queda es exactamente la varianza de las componentes descartadasvar_comp <- pca$sdev^2ecm_teo <-sapply(Ms, function(M) sum(var_comp[-(1:M)]) * (nrow(Xs) -1) /nrow(Xs) / p)stopifnot(ecm[p] <1e-10, all(diff(ecm) <0), max(abs(ecm - ecm_teo)) <1e-10)g <-ggplot(data.frame(M = Ms, ecm = ecm), aes(M, ecm)) +geom_area(fill = COL[["primary"]], alpha =0.12) +geom_line(colour = COL[["primary"]], linewidth =0.9) +geom_point(colour = COL[["primary"]], size =2) +scale_x_continuous(breaks = Ms) +labs(x ="Componentes usadas en la reconstruccion (M)", y ="Error cuadratico medio (datos estandarizados)",title ="El error de reconstruccion cae rapido: pocas componentes bastan")print(g); save_r("clase10_reconstruccion", g, w =7, h =4.5)
Predecir. Si reconstruimos la tabla con \(M\) componentes y el error cuadratico medio da cero, que podemos afirmar sobre las variables originales?
Ver respuesta
Que el bloque de variables vive, en realidad, en un espacio de \(M\) dimensiones: hay \(p - M\) direcciones sin ninguna variacion, porque algunas variables son combinaciones lineales exactas de las otras. En datos reales eso casi nunca pasa de forma perfecta, pero cuanto mas rapido cae la curva del error, mas cerca estamos de esa situacion y mas justificado esta quedarse con pocas componentes.
## 4. Decisiones practicas: escalar y cuantos componentes
Dos decisiones definen el resultado de una PCA.
Estandarizar. PCA maximiza varianza, asi que una variable con varianza enorme domina el resultado solo por la escala en que esta medida. Antes de ver el codigo, hagamos una prediccion.
Cuantos componentes. Miramos el scree plot (varianza por componente) buscando el codo, y la varianza acumulada para llegar a un umbral (80% o 90%). La proporcion de varianza explicada por los primeros \(M\) componentes es el \(R^2\) de la aproximacion. No hay una regla objetiva unica: en no supervisado es, en parte, un juicio.
Predecir. Vamos a correr PCA sobre el bloque academico sin estandarizar. Que variable pensas que va a dominar la primera componente, y por que?
Ver respuesta
La nota de admision. Esta en una escala de 0 a 200, con una varianza muchisimo mayor que las materias aprobadas (que van de 0 a ~10). Como PCA busca maxima varianza, sin estandarizar la PC1 se vuelve casi solo esa variable, no porque sea mas importante, sino por su escala.
# Escalado como descubrimiento: PCA sin escalar vs escaladacat("Varianza por variable (sin escalar):\n")print(round(sort(sapply(X, var), decreasing =TRUE), 1))pca_raw <-prcomp(X, center =TRUE, scale. =FALSE) # SIN escalarpca_std <-prcomp(X, center =TRUE, scale. =TRUE) # escaladadom <- academico[which.max(abs(pca_raw$rotation[, 1]))]cat("\nVariable que domina PC1 sin escalar:", etq[dom], "\n")stopifnot(dom =="Admission grade")l_raw <- pca_raw$rotation[, 1] *sign(pca_raw$rotation[match("Admission grade", academico), 1])l_std <- pca_std$rotation[, 1] *sign(mean(pca_std$rotation[, 1]))df <-rbind(data.frame(var =unname(etq[academico]), load =as.numeric(l_raw), panel ="Sin escalar: PC1 es casi solo la nota de admision"),data.frame(var =unname(etq[academico]), load =as.numeric(l_std), panel ="Escalada: PC1 reparte peso entre las materias"))df$var <-factor(df$var, levels =unname(etq[academico]))g <-ggplot(df, aes(var, load, fill = panel)) +geom_col(show.legend =FALSE) +facet_wrap(~panel, scales ="free_y") +scale_fill_manual(values =c(COL[["good"]], COL[["bad"]])) +labs(x =NULL, y ="Loading en PC1", title ="Sin estandarizar, la variable de mayor varianza secuestra PC1") +theme(axis.text.x =element_text(angle =45, hjust =1, size =7))print(g); save_r("clase10_escalado", g, w =11, h =4.5)callout(paste0("Salvo que todas las variables esten en la misma unidad, se estandariza cada una ","a media cero y desvio uno antes de PCA. Si no, la escala de medicion decide el resultado."),tipo ="warn", titulo ="Regla practica")
Varianza por variable (sin escalar):
Admission grade Age at enrollment
209.7 57.6
Curricular units 2nd sem (grade) Curricular units 1st sem (grade)
27.2 23.5
Curricular units 1st sem (approved) Curricular units 2nd sem (approved)
9.6 9.1
Curricular units 1st sem (enrolled) Curricular units 2nd sem (enrolled)
6.2 4.8
Variable que domina PC1 sin escalar: Nota admision
Regla practica
Salvo que todas las variables esten en la misma unidad, se estandariza cada una a media cero y desvio uno antes de PCA. Si no, la escala de medicion decide el resultado.
# Cuantos componentes: scree plot y varianza acumuladapve_s <- pca_std$sdev^2/sum(pca_std$sdev^2)cum <-cumsum(pve_s)k80 <-which(cum >=0.80)[1]; k90 <-which(cum >=0.90)[1]df <-rbind(data.frame(comp =seq_along(pve_s), val = pve_s, panel ="Individual (scree plot)"),data.frame(comp =seq_along(cum), val = cum, panel ="Acumulada"))g <-ggplot(df, aes(comp, val)) +geom_line(colour = COL[["primary"]]) +geom_point(colour = COL[["primary"]]) +facet_wrap(~panel, scales ="free_y") +scale_x_continuous(breaks =seq_along(pve_s)) +labs(x ="Componente principal", y ="Proporcion de varianza",title =paste0("El codo marca cuantas valen la pena; ", k80, " comps llegan a 80% y ", k90, " a 90%"))print(g); save_r("clase10_scree", g, w =11, h =4.5)callout(paste0("La proporcion de varianza explicada por los primeros M componentes es el ","<b>R²</b> de la mejor aproximacion de rango M a los datos. No hay una respuesta ","unica y simple a cuantos componentes usar: se elige mirando el codo y el problema."),tipo ="info", titulo ="PVE = R cuadrado, y una decision con criterio")
PVE = R cuadrado, y una decision con criterio
La proporcion de varianza explicada por los primeros M componentes es el R² de la mejor aproximacion de rango M a los datos. No hay una respuesta unica y simple a cuantos componentes usar: se elige mirando el codo y el problema.
Cuanta varianza hay para repartir, y como se reparte
Antes de mirar el scree plot conviene tener claro que es exactamente lo que el grafico esta repartiendo. La respuesta corta: la varianza total de los datos, que es un numero fijo, y las componentes se la reparten entre ellas sin dejar nada afuera.
1. La varianza total. Llamamos varianza total a la suma de las varianzas de todas las variables, una vez centradas en cero:
\(p\) es la cantidad de variables (en nuestro caso las 8 variables academicas),
\(n\) es la cantidad de estudiantes,
\(x_{ij}\) es el valor centrado de la variable \(j\) para el estudiante \(i\),
\(\operatorname{Var}(x_j)\) es la varianza de la variable \(j\).
Si ademas de centrar estandarizamos, cada variable queda con desvio 1 y por lo tanto con varianza 1. Entonces la varianza total no depende de las unidades ni de la escala de nada: es exactamente \(p\). Cada variable aporta una unidad de varianza al pozo comun. Con nuestras 8 variables academicas estandarizadas, hay 8 unidades de varianza para repartir.
2. Lo que explica cada componente es su autovalor. La varianza de los scores de la componente \(m\) es el autovalor \(\lambda_m\) de la matriz de correlaciones, el mismo objeto que aparecio cuando vimos como se calculan los loadings:
\[\lambda_m \;=\; \operatorname{Var}(z_m)\]
donde:
\(z_m\) es el vector de scores de la componente \(m\), es decir la coordenada de cada estudiante sobre esa direccion,
\(\lambda_m\) es el autovalor asociado a esa componente.
En R esos autovalores son pca$sdev^2. Como las \(p\) componentes forman una base completa del espacio, sus autovalores suman la varianza total. Rotar los ejes no crea ni destruye varianza, solo la reordena.
3. La PVE. La proporcion de varianza explicada por la componente \(m\) es simplemente su porcion del pozo:
\(\text{PVE}_m\) es la proporcion de varianza explicada por la componente \(m\), un numero entre 0 y 1,
\(\lambda_m\) es el autovalor de esa componente,
el denominador es la varianza total.
La PVE acumulada hasta \(M\) componentes es la suma de las primeras \(M\) proporciones, y si sumamos las \(p\) llegamos exactamente a 1: entre todas explican todo.
4. La conexion clave: maxima varianza y minimo error son la misma cosa. En la seccion 3.1 aproximamos los datos con las primeras \(M\) componentes y medimos el error de reconstruccion, es decir lo que quedaba afuera. Vale la pena unir las dos cuentas, porque la varianza total se parte de manera exacta en dos pedazos: lo que capturan las primeras \(M\) componentes, mas el error de reconstruccion que queda al usar solo esas \(M\). Nada se pierde y nada se agrega: uno de los pedazos crece justo lo que el otro se achica.
Y como la varianza total es un numero fijo (es \(p\), no depende de cuantas componentes usemos), la consecuencia es directa: elegir las direcciones que maximizan la varianza explicada es exactamente lo mismo que elegir el subespacio que minimiza el error de reconstruccion. No son dos criterios que casualmente dan parecido, son el mismo criterio escrito de dos formas. Esa es la razon de fondo por la que las dos definiciones de PCA que venimos usando, la direccion de maxima varianza de la seccion 2 y el subespacio mas cercano a los datos de la seccion 3.1, describen el mismo objeto visto desde dos angulos.
5. La PVE es un \(R^2\). De esa misma particion se desprende una lectura comoda. En regresion lineal, el \(R^2\) compara la variabilidad que el modelo logra explicar contra la variabilidad total de la respuesta. Aca pasa lo mismo, con la diferencia de que no hay una variable respuesta sino el conjunto entero de las \(p\) variables: la PVE acumulada de las primeras \(M\) componentes es la proporcion de la variabilidad total de los datos que la aproximacion de rango \(M\) consigue explicar. Es decir, es el \(R^2\) de la reconstruccion. Cuando leamos que dos componentes acumulan cierta proporcion, se puede decir con las mismas palabras que se usaban en regresion: esa es la fraccion de la variabilidad de las 8 variables que sobrevive al resumen en dos coordenadas.
# Verificacion numerica: varianza total, autovalores y PVEn_est <-nrow(Xs)stopifnot(ncol(Xs) == p)# 1) Con variables estandarizadas, la varianza total es exactamente p# (scale() en R estandariza con la varianza muestral, la misma que usa var())var_por_variable <-apply(Xs, 2, var)var_total <-sum(var_por_variable)cat("Varianza de cada variable estandarizada:", round(var_por_variable, 3), "\n")cat("Varianza total:", round(var_total, 6), " (p =", p, ")\n")stopifnot(abs(var_total - p) <1e-8)# 2) Los autovalores son las varianzas de las componentes y suman la varianza totalautovalores <- pca$sdev^2stopifnot(abs(sum(autovalores) - p) <1e-8)# 3) La PVE es la proporcion de cada autovalor, y acumulada llega a 1pve_manual <- autovalores /sum(autovalores)pve_acum_manual <-cumsum(pve_manual)stopifnot(all(abs(pve_manual - pve) <1e-10), abs(pve_acum_manual[p] -1) <1e-10)tabla_pve <-data.frame(Componente =1:p, Autovalor =round(autovalores, 3),PVE =round(pve_manual, 3), PVE_acumulada =round(pve_acum_manual, 3))print(tabla_pve, row.names =FALSE)
Ejercicio 1. Escribi el codigo que calcule cuantos componentes principales hacen falta para explicar al menos el 90% de la varianza del bloque academico. Pista: usa la suma acumulada (cumsum) de pca_std$sdev^2 / sum(pca_std$sdev^2).
# Solucion: componentes necesarios para un umbral de varianzaumbral <-0.90pve_s <- pca_std$sdev^2/sum(pca_std$sdev^2)n_comp <-which(cumsum(pve_s) >= umbral)[1]cat("Para explicar el", umbral *100, "% hacen falta", n_comp, "componentes\n")stopifnot(n_comp >0, n_comp <= p)
Para explicar el 90 % hacen falta 4 componentes
## 5. Metodos no lineales: t-SNE y UMAP
PCA es lineal: sus componentes son combinaciones lineales de las variables, y solo captura estructura describible con direcciones rectas. Cuando la estructura interesante es curva o esta en grupos separados de forma no lineal, PCA se queda corto.
t-SNE y UMAP construyen una representacion 2D que preserva sobre todo la vecindad local: puntos parecidos quedan cerca. Son excelentes para visualizar grupos, pero no preservan bien las distancias globales, dependen de hiperparametros (la perplexity, el numero de vecinos) y sus ejes no tienen interpretacion. No estan en el capitulo 12: los mostramos como complemento para explorar, no para inferir.
# Reduccion no lineal sobre el mismo bloque estandarizadoset.seed(SEED)emb_pca <-prcomp(Xs, center =FALSE, scale. =FALSE)$x[, 1:2]emb_tsne <-Rtsne(Xs, dims =2, perplexity =30, check_duplicates =FALSE, pca =TRUE)$Ycfg <- umap.defaults; cfg$random_state <- SEEDemb_umap <-umap(Xs, config = cfg)$layoutdf <-rbind(data.frame(x = emb_pca[, 1], y = emb_pca[, 2], metodo ="PCA (lineal)", Target = est$Target),data.frame(x = emb_tsne[, 1], y = emb_tsne[, 2], metodo ="t-SNE", Target = est$Target),data.frame(x = emb_umap[, 1], y = emb_umap[, 2], metodo ="UMAP", Target = est$Target))df$metodo <-factor(df$metodo, levels =c("PCA (lineal)", "t-SNE", "UMAP"))g <-ggplot(df, aes(x, y, colour = Target)) +geom_point(size =0.5, alpha =0.5) +facet_wrap(~metodo, scales ="free") +scale_colour_manual(values = COL_TARGET, name =NULL) +labs(x =NULL, y =NULL, title ="Tres vistas 2D de los mismos estudiantes, coloreadas por desenlace") +theme(axis.text =element_blank(), axis.ticks =element_blank())print(g); save_r("clase10_no_lineales", g, w =15, h =5)callout(paste0("t-SNE y UMAP deforman las distancias globales para preservar las locales: la ","distancia entre dos grupos lejanos en el grafico <b>no</b> es interpretable, y cambiar ","un hiperparametro cambia el dibujo. Usalos para descubrir grupos, nunca para medir."),tipo ="warn", titulo ="Cuidado con la interpretacion")
Cuidado con la interpretacion
t-SNE y UMAP deforman las distancias globales para preservar las locales: la distancia entre dos grupos lejanos en el grafico no es interpretable, y cambiar un hiperparametro cambia el dibujo. Usalos para descubrir grupos, nunca para medir.
## 6. Ejercicio de cierre: perfiles latentes en attrition
Ahora aplicas todo lo visto a un dataset que ya conoces: el de rotacion de personal (attrition) de las clases anteriores. La idea que cierra el capitulo es que muchas tecnicas, entre ellas el clustering, se pueden correr sobre los primeros \(M\) scores de PCA en vez de sobre las variables originales, tratandolos como una version menos ruidosa de los datos.
Consigna. Sobre el bloque de variables de carrera de los empleados: (1) estandarizar y calcular la PCA; (2) quedarse con los primeros componentes; (3) correr K-means sobre esos scores para descubrir perfiles latentes de empleados; (4) perfilar cada grupo y cruzarlo con la rotacion. Intenta completarlo antes de mirar la solucion.
# Solucion del ejercicio de cierre (attrition)hr <- readr::read_csv(URL_DATOS, show_col_types =FALSE)carrera <-c("Age", "TotalWorkingYears", "MonthlyIncome", "JobLevel", "YearsAtCompany","YearsInCurrentRole", "YearsWithCurrManager", "YearsSinceLastPromotion", "NumCompaniesWorked")# 1-2. PCA sobre el bloque estandarizado y primeros M scores (version menos ruidosa)M <-3Zc <-prcomp(as.data.frame(hr)[, carrera], center =TRUE, scale. =TRUE)$x[, 1:M]# 3. K-means sobre los scoresset.seed(SEED); k <-3km <-kmeans(Zc, centers = k, nstart =10)hr$perfil <-factor(km$cluster)stopifnot(length(unique(km$cluster)) == k)# 4. Perfilar: plano PC1-PC2 y rotacion por perfilsigno <-sign(cor(Zc[, 1], hr$TotalWorkingYears))dfc <-data.frame(PC1 = Zc[, 1] * signo, PC2 = Zc[, 2], perfil = hr$perfil)g1 <-ggplot(dfc, aes(PC1, PC2, colour = perfil)) +geom_point(size =0.9, alpha =0.6) +scale_colour_manual(values =unname(COL[c("primary", "accent", "good")]), name ="Perfil") +labs(x ="PC1 (senioridad)", y ="PC2 (movilidad)", title ="K-means sobre los scores: perfiles latentes")print(g1); save_r("clase10_perfiles", g1, w =7, h =5)tab <-prop.table(table(hr$perfil, hr$Attrition), margin =1)print(round(tab, 3))print(round(t(aggregate(as.data.frame(hr)[, carrera], by =list(perfil = km$cluster), FUN = mean)[, -1]), 1))callout(paste0("Los perfiles son un <b>punto de partida para hipotesis</b>, no una verdad del dataset. ","El clustering no es robusto: cambia con la estandarizacion, el numero de componentes y k. ","Conviene repetirlo con distintas opciones y quedarse con los patrones que persisten."),tipo ="success", titulo ="Como reportar perfiles con cuidado")
Los perfiles son un punto de partida para hipotesis, no una verdad del dataset. El clustering no es robusto: cambia con la estandarizacion, el numero de componentes y k. Conviene repetirlo con distintas opciones y quedarse con los patrones que persisten.
## Autoevaluacion
Tres preguntas para chequear las ideas centrales. Respondelas antes de abrir cada solucion.
1. Que es la primera componente principal, geometricamente?
Ver respuesta
La direccion a lo largo de la cual los datos varian mas y, equivalentemente, la recta mas cercana a todos los puntos. No es una de las variables originales: es una combinacion lineal de todas.
2. Antes de PCA, cuando conviene estandarizar las variables?
Ver respuesta
Cuando estan en unidades distintas o tienen varianzas muy dispares. Si no, la variable de mayor varianza domina las componentes solo por su escala, como pasa aca con la nota de admision.
3. Que mide la proporcion de varianza explicada (PVE) por los primeros M componentes?
Ver respuesta
El \(R^2\) de la mejor aproximacion de rango M a los datos: que fraccion de la variacion total conservan esos M componentes. Equivalentemente, uno menos la proporcion de error de reconstruccion.
Caso aplicado 4. Reducir 32 medidas de comportamiento a unas pocas dimensiones
Arguello y Crescenzi (2019) analizaron con PCA hasta 32 medidas de comportamiento por sesion (cantidad de queries, clics, tiempos de permanencia, scrolls, mouseovers) capturadas en tres estudios de busqueda de informacion. La pregunta no era predecir, era entender: que fenomenos latentes capturan realmente esas medidas, y como se relacionan con las percepciones que reporta la persona al terminar la tarea (carga de trabajo, dificultad, presion de tiempo, engagement) (Arguello y Crescenzi, ICTIR 19, paginas 177 a 184; PDF abierto).
Vale la pena mirar las decisiones metodologicas, porque son las mismas que tomamos en esta clase:
Escalado. Trabajaron sobre la matriz de correlacion y no sobre la de covarianza, porque las medidas estan en escalas muy distintas: las queries se cuentan en decenas y las duraciones en cientos de segundos. Eso equivale a estandarizar antes de aplicar PCA, exactamente lo que vimos en la seccion 4.
Cuantas componentes. Retuvieron las de autovalor mayor a 1 y, ademas, exigieron que cada componente tuviera al menos dos variables con loadings altos (0.50 o mas en valor absoluto). El segundo criterio evita quedarse con componentes que en los hechos son una sola variable disfrazada.
Rotacion varimax. Rotaron la solucion para mejorar la interpretabilidad, favoreciendo que cada variable cargue fuerte en una componente y debil en las demas. Es practica estandar en este tipo de trabajos; ISLP no la desarrolla, pero conviene saber que existe cuando lean papers.
Nombrado. Cada componente recibio un nombre segun las medidas que cargaban fuerte, por ejemplo abandono de queries o ritmo de interaccion.
Uso posterior. Los scores de las componentes entraron como predictores en un modelo de regresion multinivel. Es el flujo completo: primero reducir la dimension, despues modelar con las componentes.
Un resultado ilustrativo: en uno de los estudios, el PCA con rotacion varimax dio una solucion de seis componentes que explico el 76 por ciento de la varianza total de 24 medidas de comportamiento. Del mismo modo que en nuestro caso unas pocas componentes alcanzan para resumir el desempeno academico, ahi seis dimensiones alcanzaron para resumir como una persona se comporta frente a un buscador.
## Cierre
Recorrimos PCA de punta a punta: de la redundancia entre variables correlacionadas, a la intuicion en 2D (la recta de maxima varianza, que es tambien la mas cercana), a su generalizacion a muchas dimensiones con el biplot, a las decisiones de escalar y elegir cuantos componentes, a la mencion de los metodos no lineales, y al cierre aplicando PCA + clustering sobre datos que ya conocias.
Para seguir. El capitulo 12 cubre ademas otros usos de las componentes (Seccion 12.2.5), la imputacion de faltantes por matrix completion (12.3) y el detalle de los metodos de clustering (12.4). En las Clases 11 y 12 la idea de representar datos en pocas dimensiones reaparece como embeddings de texto.
Apendice: hoja de referencia
concepto
formula o idea
Componente principal
\(Z_m = \sum_j \phi_{jm} X_j\), con \(\sum_j \phi_{jm}^2 = 1\), de maxima varianza
Loadings
los \(\phi_{jm}\): cuanto pesa cada variable en la componente \(m\)
Scores
\(z_{im}\): la coordenada de la observacion \(i\) en la componente \(m\)
Doble lectura
maxima varianza \(\equiv\) recta/plano mas cercano a los puntos
PVE
proporcion de varianza explicada; igual al \(R^2\) de la aproximacion
Estandarizar
media 0 y desvio 1 por variable, salvo unidad comun
Cuantos componentes
codo del scree plot o umbral de varianza acumulada (80 a 90%)
PCA + clustering
correr K-means sobre los primeros M scores (version menos ruidosa)
Ejecutar / compartir esta notebook. Abrila en Google Colab y elegi Entorno de ejecucion -> Ejecutar todo; no hace falta instalar nada localmente.
Fuentes: ISLP (James et al., 2023), capitulo 12; dataset de desercion UCI id 697 (CC BY 4.0). Cada numero de esta notebook se recalcula y se verifica con stopifnot.
