Principal Component Analysis (PCA): método estadístico que permite simplificar la complejidad de espacios muestrales con muchas dimensiones a la vez que conserva su información.
Dada una muestra de n individuos cada uno con p variables \((X_1,X_2,…,X_p)\), es decir, el espacio muestral tiene p dimensiones.
El PCA permite encontrar un número de factores subyacentes (z<p) que explican aproximadamente lo mismo que las p variables originales.
Cada z recibe el nombre de componente principal.
Una de las aplicaciones de PCA es la reducción de dimensionalidad (variables), perdiendo la menor cantidad de información (varianza) posible:
El PCA sirve como herramienta para la visualización de datos
El análisis en componentes principales (ACP) se utiliza para describir tablas que tienen en las filas las unidades estadísticas, generalmente denominadas, “individuos”, y en las columnas las variables de tipo continuo que se han medido sobre los individuos.
Los eigenvectores y eigenvalores son números y vectores asociados a matrices cuadradas. Dada una matriz A de \(n \times n\), su eigenvector \(v\) de \(n \times\) tal que,
\[Av=\lambda v\]
\[ \begin{pmatrix} 2 & 3 \\ 2 & 1 \end{pmatrix} \begin{pmatrix} 3 \\ 2 \end{pmatrix} = \begin{pmatrix} 12 \\ 8 \end{pmatrix} = 4 \begin{pmatrix} 3 \\ 2 \end{pmatrix} \]
en R los puede calcular con la funcion eigen()
Solo las matrices cuadradas tienen eigenvectores, pero no todas las matrices cuadradas los tienen.
Dada una matriz \(n \times n\) con eigenvectores, el número existente de ellos es n.
Todos los eigenvectores de una matriz son perpendiculares. Esto significa que podemos expresar los datos respecto a estos eigenvectores.
Un eigenvalor > 1 indica que la componente principal explica más varianza de lo que lo hace una de las variables originales, estando los datos estandarizados.
Una forma intuitiva de entender el proceso de PCA consiste en interpretar las componentes principales desde un punto de vista geométrico.
Supóngase un conjunto de observaciones para las que se dispone de dos variables (\(X_1, X_2\)).
El vector que define la primera componente principal (\(Z_1\)) sigue la dirección en la que las observaciones varían más (linea roja).
La proyección de cada observación sobre esa dirección equivale al valor de la primera componente para dicha observación (principal component scores, \(z_{i_1}\)).
La segunda componente (\(Z_2\)) sigue la segunda dirección en la que los datos muestran mayor varianza y que no está correlacionada con la primera componente.
La condición de no correlación entre componentes principales equivale a decir que sus direcciones son perpendiculares/ortogonales.
Cada componente principal (\(Z_i\)) se obtiene por combinación lineal de las variables originales.
\[Z_1=\phi_{11} X_1 + \phi_{21}X_2+ \dots + \phi_{p1}X_p\]
con la condición siguiente:
\[\displaystyle\sum_{j=1}^p \phi^{2}_{j1}=1\] En general, se tiene que:
\[Z_{m}=\sum_{j=1}^{p} \phi_{j m} X_{j}\]
\(\phi_{1m},\phi_{12},\dots,\phi_{pm}\) son las cargas o loadings de los componentes principales
La primera componente principal (\(Z_1\)) es aquella cuya dirección refleja o contiene la mayor variabilidad en los datos.
Este vector define la línea lo más próxima posible a los datos y que minimiza la suma de las distancias perpendiculares entre cada dato y la línea representada por la componente (usando como medida de cercanía el promedio de la distancia euclídea al cuadrado):
\[z_{i1}=\phi_{11} x_{i 1}+\phi_{21} x_{i 2}+\ldots+\phi_{p 1} x_{i p}\]
donde \(\phi_{11}\) corresponde al primer loading de la primera componente principal.
En otras palabras, el vector de cargas de la primera componente principal resuelve el problema de optimización
\[\underset{\phi_{11}, \ldots, \phi_{p 1}}{\operatorname{maximize}}\left\{\frac{1}{n} \sum_{i=1}^{n}\left(\sum_{j=1}^{p} \phi_{j 1} x_{i j}\right)^{2}\right\} \\ \text { sujeto a } \sum_{j=1}^{p} \phi_{j 1}^{2}=1\]
¿Cuánta información presente en el set de datos original se pierde al proyectar las observaciones en un espacio de menor dimensión?
¿Cuanta información es capaz de capturar cada una de las componentes principales obtenidas?
Asumiendo que las variables se han normalizado para tener media cero, la varianza total presente en el set de datos se define como
\[ \displaystyle\sum_{j=1}^p Var(X_j) = \displaystyle\sum_{j=1}^p \frac{1}{n} \displaystyle\sum_{i=1}^n x^{2}_{ij}\]
\[\frac{1}{n} \sum_{i=1}^n z^{2}_{im} = \frac{1}{n} \sum_{i=1}^n \left( \sum_{j=1}^p \phi_{jm}x_{ij} \right)^2\]
\[\frac{\sum_{i=1}^n \left( \sum_{j=1}^p \phi_{jm}x_{ij} \right)^2} {\sum_{j=1}^p \sum_{i=1}^n x^{2}_{ij}}\]
Las librerías en R y python utilizadas en este documento son:
library(readr)
library(ggplot2)
library(reticulate)# Tratamiento de datos
# ==============================================================================
import numpy as np
import pandas as pd
import statsmodels.api as sm
# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
import matplotlib.font_manager
from matplotlib import style
style.use('ggplot') or plt.style.use('ggplot')
# Preprocesado y modelado
# ==============================================================================
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.preprocessing import scale
# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('ignore')El ingreso per cápita (PCI) es una medida común que evalúa el ingreso promedio obtenido por persona en una determinada región y durante un año específico. Se calcula sumando todos los ingresos dentro de la región y dividiéndolos por la población de la región, que incluye tanto a niños como a adultos. Si bien esta medida tiene su limitación (p. ej., no tiene en cuenta la inflación), PCI es un método útil para comparar el nivel de vida y la calidad de vida de varias regiones: cuanto mayor es el ingreso per cápita, más rica es la población.
A nivel nacional, el ingreso per cápita fue de $35,384 en 2020 según la Oficina del Censo de los Estados Unidos. Mirando más precisamente a cada estado, el PCI puede diferir mucho de uno a otro. Varios factores pueden afectar el PCI, incluidos, entre otros, la educación, el tipo de población, el bienestar económico general y la legislación estatal específica.
Las variables son:
Aprenda usted mismo cómo varía el ingreso per cápita de un estado a otro al verificar los datos en los gráficos.
datos = pd.read_csv("income_usa.csv")
State = datos["State"].tolist()
datos = datos.drop(['State'], axis=1)
datos.index = State
datos.info() #informacion de los datos
datos.isna().sum()# cuantos na hay
datos = datos.dropna() # borrar filas con na
datos.isna().sum()
datos.head(4) #primeros 4 filasLos dos principales aspectos a tener en cuenta cuando se quiere realizar un PCA es identificar el valor promedio y dispersión de las variables.
La media de las variables muestra que hay tres veces más secuestros que asesinatos y 8 veces más asaltos que secuestros.
print('----------------------')
print('Media de cada variable')
print('----------------------')
datos.mean(axis=0, numeric_only=True)print('-------------------------')
print('Varianza de cada variable')
print('-------------------------')
datos.var(axis=0, numeric_only=True)Si no se estandarizan las variables para que tengan media cero y desviación estándar de uno antes de realizar el estudio PCA, la variable Assault, que tiene una media y dispersión muy superior al resto, dominará la mayoría de las componentes principales.
La clase sklearn.decomposition.PCA incorpora las principales funcionalidades que se necesitan a la hora de trabajar con modelos PCA. El argumento n_components determina el número de componentes calculados. Si se indica None, se calculan todas las posibles (min(filas, columnas) - 1).
Por defecto, PCA() centra los valores pero no los escala. Esto es importante ya que, si las variables tienen distinta dispersión, como en este caso, es necesario escalarlas. Una forma de hacerlo es combinar un StandardScaler() y un PCA() dentro de un pipeline. Para más información sobre el uso de pipelines consultar Pipeline y ColumnTransformer.
# Entrenamiento modelo PCA con escalado de los datos
pca_pipe = make_pipeline(StandardScaler(), PCA())
pca_pipe.fit(datos)
# Se extrae el modelo entrenado del pipeline
modelo_pca = pca_pipe.named_steps['pca']Una vez entrenado el objeto PCA, pude accederse a toda la información de las componentes creadas. components_ contiene el valor de los loadings \(\phi\) que definen cada componente (eigenvector). Las filas se corresponden con las componentes principals (ordenadas de mayor a menor varianza explicada). Las filas se corresponden con las variables de entrada.
# Se combierte el array a dataframe para añadir nombres a los ejes.
pd.DataFrame(
data = modelo_pca.components_,
columns = datos.columns,
index = ['PC1', 'PC2', 'PC3', 'PC4']
)Analizar con detalle el vector de loadings que forma cada componente puede ayudar a interpretar qué tipo de información recoge cada una de ellas. Por ejemplo, la primera componente es el resultado de la siguiente combinación lineal de las variables originales:
\[PC1=-Income\_per\_Capita \ 0.51941072 +\ Poverty\_Rate\ 0.49311765 \ -Median\_Household\_Income\ 0.58157815 - Minimum\_Hourly\_Wage\ 0.38576462\]
La influencia de las variables en cada componente debe analizarse visualmente con un gráfico de tipo heatmap.
# Heatmap componentes
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 2))
componentes = modelo_pca.components_
plt.imshow(componentes.T, cmap='viridis', aspect='auto')
plt.yticks(range(len(datos.columns)), datos.columns)
plt.xticks(range(len(datos.columns)), np.arange(modelo_pca.n_components_) + 1)
plt.grid(False)
plt.colorbar();
plt.show()Biplot
from pca import pca
import pandas as pd
###########################################################
# SETUP DATA
###########################################################
# Load sample data, represent the data as a pd.DataFrame
from sklearn.datasets import load_iris
iris = load_iris()
X = pd.DataFrame(data=iris.data,
columns=iris.feature_names)
X.columns = ["sepal_length", "sepal_width", "petal_length", "petal_width"]
y = pd.Categorical.from_codes(iris.target,
iris.target_names)
model = pca(n_components=0.95)
# ... or explicitly specify the number of PCs
model = pca(n_components=2)
# Fit and transform
results = model.fit_transform(X=datos, row_labels=State)
# Plot the explained variance
fig, ax = model.plot()
# Scatter the first two PCs
fig, ax = model.scatter()
# Create a biplot
fig, ax = model.biplot(n_feat=4)
plt.show()from pca import pca
import pandas as pd
###########################################################
# SETUP DATA
###########################################################
# Load sample data, represent the data as a pd.DataFrame
from sklearn.datasets import load_iris
iris = load_iris()
X = pd.DataFrame(data=iris.data,
columns=iris.feature_names)
X.columns = ["sepal_length", "sepal_width", "petal_length", "petal_width"]
y = pd.Categorical.from_codes(iris.target,
iris.target_names)
###########################################################
# COMPUTE AND VISUALIZE PCA
###########################################################
# Initialize the PCA, either reduce the data to the number of
# principal components that explain 95% of the total variance...
model = pca(n_components=0.95)
# ... or explicitly specify the number of PCs
model = pca(n_components=2)
# Fit and transform
results = model.fit_transform(X=X, row_labels=y)
# Plot the explained variance
fig, ax = model.plot()
# Scatter the first two PCs
fig, ax = model.scatter()
# Create a biplot
fig, ax = model.biplot(n_feat=4)Una vez calculadas las componentes principales, se puede conocer la varianza explicada por cada una de ellas, la proporción respecto al total y la proporción de varianza acumulada. Esta información está almacenada en los atributos explained_variance_ y explained_variance_ratio_ del modelo.
Creando las variables CP y Var para usarlas en R
cp = np.arange(modelo_pca.n_components_) + 1,
var = modelo_pca.explained_variance_ratio_Gráfico en R de la varianza
library(reticulate)
cp1 = py$cp
var1 = py$var
df1 <-data.frame(Componente = cp1[[1]],
Varianza = var1)
library(ggplot2)
ggplot(data=df1, aes(x=Componente, y=Varianza)) +
geom_bar(stat="identity", fill="steelblue") +
geom_text(aes(label= round(Varianza,3)), position=position_dodge(width=0.9), vjust=-0.25) +
labs(title ="Varianza", x ="Componente", y ="Varianza")+
theme_bw()En este caso, la primera componente explica el 66% de la varianza observada en los datos y la segunda el 20%. Las dos últimas componentes no superan por separado el 1.4% de varianza explicada.
Gráfico en R de la varianza acumulada
ggplot(data=df1, aes(x=Componente, y=cumsum(round(Varianza,2)))) +
geom_line(color="red") +
geom_point(aes(x =Componente, y= cumsum(round(Varianza,2))), color="red", size= 2)+
geom_text(aes(label= cumsum(round(Varianza,2))), position=position_dodge(width=0.9), vjust=-.3) +
labs(title ="Varianza acumulada", x ="Componente", y = "Varianza")+
theme_bw()Una vez entrenado el modelo, con el método transform() se puede reducir la dimensionalidad de nuevas observaciones proyectándolas en el espacio definido por las componentes.
# Proyección de las observaciones de entrenamiento
proyecciones = pca_pipe.transform(X=datos)
proyecciones = pd.DataFrame(
proyecciones,
columns = ['PC1', 'PC2', 'PC3', 'PC4'],
index = datos.index
)
proyecciones.head()La transformación es el resultado de multiplicar los vectores que definen cada componente con el valor de las variables. Puede calcularse de forma manual:
proyecciones = np.dot(modelo_pca.components_, scale(datos).T)
proyecciones = pd.DataFrame(proyecciones, index = ['PC1', 'PC2', 'PC3', 'PC4'])
proyecciones = proyecciones.transpose().set_index(datos.index)
proyecciones.head()Puede revertirse la transformación y reconstruir el valor inicial con el método inverse_transform(). Es importante tener en cuenta que, la reconstrucción, solo será completa si se han incluido todas las componentes.
# Recostruccion de las proyecciones
recostruccion = pca_pipe.inverse_transform(X=proyecciones)
recostruccion = pd.DataFrame(
recostruccion,
columns = datos.columns,
index = datos.index
)
print('---------------------')
print('Valores reconstruidos')
print('---------------------')
datos.head()
print('------------------')
print('Valores originales')
print('------------------')
recostruccion.head()