PCA con Python

Joaquín Amat Rodrigo

Diciembre, 2020

Más sobre ciencia de datos: [Link] ([Link]

Regresión lineal con Python ([Link]

python)
Regularización Ridge, Lasso y Elastic Net con Python
([Link]
Machine learning con Python y Scikit-learn
([Link]
Random Forest con Python
([Link]
Gradient Boosting con Python
([Link]

Introducción

El análisis de componentes principales (Principal Component Analysis PCA) es un

método de reducción de dimensionalidad que permite simplificar la complejidad de
espacios con múltiples dimensiones a la vez que conserva su información.

Supóngase que existe una muestra con n individuos cada uno con p variables (X1 , X2 ,
..., Xp ), es decir, el espacio muestral tiene p dimensiones. PCA permite encontrar un
número de factores subyacentes (z < p) que explican aproximadamente lo mismo que
las p variables originales. Donde antes se necesitaban p valores para caracterizar a
cada individuo, ahora bastan z valores. Cada una de estas z nuevas variables recibe el
nombre de componente principal.
El método de PCA permite por lo tanto "condensar" la información aportada por
múltiples variables en solo unas pocas componentes. Aun así, no hay que olvidar que
sigue siendo necesario disponer del valor de las variables originales para calcular las
componentes. Dos de las principales aplicaciones del PCA son la visualización y el
preprocesado de predictores previo ajuste de modelos supervisados.

La librería scikitlearn ([Link]

[Link]/stable/modules/generated/[Link]) contiene la clase
[Link] que implementa la mayoría de las funcionalidades
necesarias para crear y utilizar modelos PCA. Para visualizaciones, Yellowbrick
([Link] ofrece funcionalidades extra.
Álgebra lineal

En esta sección se describen dos de los conceptos matemáticos que se aplican en el

PCA: eigenvectors y eigenvalues. Se trata simplemente de una descripción intuitiva con
la única finalidad de facilitar el entendimiento del cálculo de componentes principales.

Eigenvectors

Los eigenvectors son un caso particular de multiplicación entre una matriz y un vector.
Obsérvese la siguiente multiplicación:
2 3 3 12 3
( )x( ) = ( ) = 4 x( )
2 1 2 8 2

El vector resultante de la multiplicación es un múltiplo entero del vector original. Los

eigenvectors de una matriz son todos aquellos vectores que, al multiplicarlos por dicha
matriz, resultan en el mismo vector o en un múltiplo entero del mismo. Los eigenvectors
tienen una serie de propiedades matemáticas específicas:

Los eigenvectors solo existen para matrices cuadradas y no para todas. En el caso
de que una matriz n x n tenga eigenvectors, el número de ellos es n.

Si se escala un eigenvector antes de multiplicarlo por la matriz, se obtiene un

múltiplo del mismo eigenvector. Esto se debe a que, si se escala un vector
multiplicándolo por cierta cantidad, lo único que se consigue es cambiar su longitud
pero la dirección es la misma.

Todos los eigenvectors de una matriz son perpendiculares (ortogonales) entre ellos,
independientemente de las dimensiones que tengan.

Dada la propiedad de que multiplicar un eigenvector solo cambia su longitud pero no su

naturaleza de eigenvector, es frecuente escalarlos de tal forma que su longitud sea 1.
De este modo se consigue que todos ellos estén estandarizados. A continuación se
muestra un ejemplo:

3 −− −−−− −−
El eigenvector ( ) tiene una longitud de 2 2
√3 + 2 = √13 . Si se divide cada
2

dimensión entre la longitud del vector, se obtiene el eigenvector estandarizado con

longitud 1.
−−
3 −− 3/√13
( ) ÷ √13 = ( −−)
2 2/√13
Eigenvalue

Cuando se multiplica una matriz por alguno de sus eigenvectors se obtiene un múltiplo
del vector original, es decir, el resultado es ese mismo vector multiplicado por un
número. Al valor por el que se multiplica el eigenvector resultante se le conoce como
eigenvalue. A todo eigenvector le corresponde un eigenvalue y viceversa.

En el método PCA, cada una de las componentes se corresponde con un eigenvector, y

el orden de componente se establece por orden decreciente de eigenvalue. Así pues, la
primera componente es el eigenvector con el eigenvalue más alto.

Interpretación geométrica de las

componentes principales

Una forma intuitiva de entender el proceso de PCA es 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 (X1 , X2 ). El vector que define
la primera componente principal (Z1 ) sigue la dirección en la que las observaciones
tienen más varianza (línea 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
score, zi1 ).
La segunda componente (Z2 ) 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.
Cálculo de las componentes
principales

Cada componente principal (Zi ) se obtiene por combinación lineal de las variables
originales. Se pueden entender como nuevas variables obtenidas al combinar de una
determinada forma las variables originales. La primera componente principal de un
grupo de variables (X1 , X2 , ..., Xp ) es la combinación lineal normalizada de dichas
variables que tiene mayor varianza:

Z1 = ϕ11 X1 + ϕ21 X2 +. . . +ϕp1 Xp

Que la combinación lineal sea normalizada implica que:

p

2
∑ϕ = 1
j1

j=1

Los términos ϕ 11 , ..., ϕ 1p reciben en el nombre de loadings y son los que definen las
componentes. Por ejemplo, ϕ 11 es el loading de la variable X1 de la primera
componente principal. Los loadings pueden interpretarse como el peso/importancia que
tiene cada variable en cada componente y, por lo tanto, ayudan a conocer que tipo de
información recoge cada una de las componentes.

Dado un set de datos X con n observaciones y p variables, el proceso a seguir para

calcular la primera componente principal es:

Centrar las variables: se resta a cada valor la media de la variable a la que

pertenece. Con esto se consigue que todas las variables tengan media cero.

Se resuelve un problema de optimización para encontrar el valor de los loadings con

los que se maximiza la varianza. Una forma de resolver esta optimización es
mediante el cálculo de eigenvector-eigenvalue de la matriz de covarianzas.

Una vez calculada la primera componente (Z1 ), se calcula la segunda (Z2 ) repitiendo el
mismo proceso pero añadiendo la condición de que la combinación lineal no pude estar
correlacionada con la primera componente. Esto equivale a decir que Z1 y Z2 tienen
que ser perpendiculares. EL proceso se repite de forma iterativa hasta calcular todas las
posibles componentes (min(n-1, p)) o hasta que se decida detener el proceso. El orden
de importancia de las componentes viene dado por la magnitud del eigenvalue asociado
a cada eigenvector.
Escalado de las variables

El proceso de PCA identifica las direcciones con mayor varianza. Como la varianza de
una variable se mide en sus mismas unidades elevadas al cuadrado, si antes de
calcular las componentes no se estandarizan todas las variables para que tengan media
cero y desviación estándar de uno, aquellas variables cuya escala sea mayor dominarán
al resto. De ahí que sea recomendable estandarizar siempre los datos.

Reproducibilidad de las componentes

El proceso de PCA estándar es determinista, genera siempre las mismas componentes

principales, es decir, el valor de los loadings resultantes es el mismo. La única diferencia
que puede darse es que el signo de todos los loadings esté invertido. Esto es así porque
el vector de loadings determina la dirección de la componente, y dicha dirección es la
misma independientemente del signo (la componente sigue una línea que se extiende
en ambas direcciones). Del mismo modo, el valor específico de las componentes
obtenido para cada observación (principal component scores) es siempre el mismo, a
excepción del signo.

Influencia de outliers

Al trabajar con varianzas, el método PCA es muy sensible a outliers, por lo que es
recomendable estudiar si los hay. La detección de valores atípicos con respecto a una
determinada dimensión es algo relativamente sencillo de hacer mediante
comprobaciones gráficas. Sin embargo, cuando se trata con múltiples dimensiones el
proceso se complica. Por ejemplo, considérese un hombre que mide 2 metros y pesa 50
kg. Ninguno de los dos valores es atípico de forma individual, pero en conjunto se
trataría de un caso muy excepcional.
Proporción de varianza explicada

Una de las preguntas más frecuentes que surge tras realizar un PCA es: ¿Cuánta
información presente en el set de datos original se pierde al proyectar las observaciones
en un espacio de menor dimensión? o lo que es lo mismo ¿Cuanta información es capaz
de capturar cada una de las componentes principales obtenidas? Para contestar a estas
preguntas se recurre a la proporción de varianza explicada por cada componente
principal.

Asumiendo que las variables se han normalizado para tener media cero, la varianza
total presente en el set de datos se define como
p p n
1
2
∑ V ar(Xj ) = ∑ ∑x
ij
n
j=1 j=1 i=1

y la varianza explicada por la componente m es

2
n n p
1 1
2
∑z = ∑ (∑ ϕjm x ij )
im
n n
i=1 i=1 j=1

Por lo tanto, la proporción de varianza explicada por la componente m viene dada por el
ratio
2
n p
∑ (∑ ϕjm x ij )
i=1 j=1

p n 2
∑ ∑ x
j=1 i=1 ij

Tanto la proporción de varianza explicada, como la proporción de varianza explicada

acumulada, son dos valores de gran utilidad a la hora de decidir el número de
componentes principales a utilizar en los análisis posteriores. Si se calculan todas las
componentes principales de un set de datos, entonces, aunque transformada, se está
almacenando toda la información presente en los datos originales. El sumatorio de la
proporción de varianza explicada acumulada de todas las componentes es siempre 1.
Número óptimo de componentes principales

Por lo general, dada una matriz de datos de dimensiones n x p, el número de

componentes principales que se pueden calcular es como máximo de n-1 o p (el menor
de los dos valores es el limitante). Sin embargo, siendo el objetivo del PCA reducir la
dimensionalidad, suelen ser de interés utilizar el número mínimo de componentes que
resultan suficientes para explicar los datos. No existe una respuesta o método único que
permita identificar cual es el número óptimo de componentes principales a utilizar. Una
forma de proceder muy extendida consiste en evaluar la proporción de varianza
explicada acumulada y seleccionar el número de componentes mínimo a partir del cual
el incremento deja de ser sustancial.

Ejemplo PCA

Librerías

Las librerías utilizadas en este documento son:

 In [39]: # Tratamiento de datos
# ====================================================================
==========
import numpy as np
import pandas as pd
import [Link] as sm

# Gráficos
# ====================================================================
==========
import [Link] as plt
import matplotlib.font_manager
from matplotlib import style
[Link]('ggplot') or [Link]('ggplot')

# Preprocesado y modelado
# ====================================================================
==========
from [Link] import PCA
from [Link] import make_pipeline
from [Link] import StandardScaler
from [Link] import scale

# Configuración warnings
# ====================================================================
==========
import warnings
[Link]('ignore')

Datos

El set de datos USArrests contiene el porcentaje de asaltos (Assault), asesinatos

(Murder) y secuestros (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 rurales (UrbanPoP).

In [40]: USArrests = [Link].get_rdataset("USArrests", "datasets")

datos = [Link]
 In [41]: [Link]()

<class '[Link]'>
Index: 50 entries, Alabama to Wyoming
Data columns (total 4 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 Murder 50 non-null float64
1 Assault 50 non-null int64
2 UrbanPop 50 non-null int64
3 Rape 50 non-null float64
dtypes: float64(2), int64(2)
memory usage: 2.0+ KB

In [42]: [Link](4)

Out[42]:
Murder Assault UrbanPop Rape

Alabama 13.2 236 58 21.2

Alaska 10.0 263 48 44.5

Arizona 8.1 294 80 31.0

Arkansas 8.8 190 50 19.5

Exploración inicial

Los 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.

In [43]: print('----------------------')
print('Media de cada variable')
print('----------------------')
[Link](axis=0)

----------------------
Media de cada variable
----------------------

Out[43]: Murder 7.788

Assault 170.760
UrbanPop 65.540
Rape 21.232
dtype: float64
La varianza es muy distinta entre las variables, en el caso de Assault, la varianza es
varios órdenes de magnitud superior al resto.

In [44]: print('-------------------------')
print('Varianza de cada variable')
print('-------------------------')
[Link](axis=0)

-------------------------
Varianza de cada variable
-------------------------

Out[44]: Murder 18.970465

Assault 6945.165714
UrbanPop 209.518776
Rape 87.729159
dtype: float64

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.

Modelo PCA

La clase [Link] 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
([Link]
y-ColumnTransformer).
 In [45]: # 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']

Interpretación

Una vez entrenado el objeto PCA , pude accederse a toda la información de las
componentes creadas.

components_ contiene el valor de los loadings ϕ 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.

In [46]: # Se combierte el array a dataframe para añadir nombres a los ejes.

[Link](
data = modelo_pca.components_,
columns = [Link],
index = ['PC1', 'PC2', 'PC3', 'PC4']
)

Out[46]:
Murder Assault UrbanPop Rape

PC1 0.535899 0.583184 0.278191 0.543432

PC2 0.418181 0.187986 -0.872806 -0.167319

PC3 -0.341233 -0.268148 -0.378016 0.817778

PC4 0.649228 -0.743407 0.133878 0.089024

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 = 0.535899 Murder + 0.583184 Assault + 0.278191 UrbanPop + 0.543432 Rape

Los pesos asignados en la primera componente a las variables Assault, Murder y Rape
son aproximadamente iguales entre ellos y superiores al asignado a UrbanPoP. Esto
significa que la primera componente recoge mayoritariamente la información
correspondiente a los delitos. En la segunda componente, es la variable UrbanPoP es la
que tiene con diferencia mayor peso, por lo que se corresponde principalmente con el
nivel de urbanización del estado. Si bien en este ejemplo la interpretación de las
componentes es bastante clara, no en todos los casos ocurre lo mismo, sobre todo a
medida que aumenta el número de variables.

La influencia de las variables en cada componente analizarse visualmente con un

gráfico de tipo heatmap.

In [47]: # Heatmap componentes

# ====================================================================
==========
fig, ax = [Link](nrows=1, ncols=1, figsize=(4, 2))
componentes = modelo_pca.components_
[Link](componentes.T, cmap='viridis', aspect='auto')
[Link](range(len([Link])), [Link])
[Link](range(len([Link])), [Link](modelo_pca.n_component
s_) + 1)
[Link](False)
[Link]();

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.
In [48]: # Porcentaje de varianza explicada por cada componente
# ====================================================================
==========
print('----------------------------------------------------')
print('Porcentaje de varianza explicada por cada componente')
print('----------------------------------------------------')
print(modelo_pca.explained_variance_ratio_)

fig, ax = [Link](nrows=1, ncols=1, figsize=(6, 4))

[Link](
x = [Link](modelo_pca.n_components_) + 1,
height = modelo_pca.explained_variance_ratio_
)

for x, y in zip([Link](len([Link])) + 1, modelo_pca.explaine

d_variance_ratio_):
label = round(y, 2)
[Link](
label,
(x,y),
textcoords="offset points",
xytext=(0,10),
ha='center'
)

ax.set_xticks([Link](modelo_pca.n_components_) + 1)
ax.set_ylim(0, 1.1)
ax.set_title('Porcentaje de varianza explicada por cada componente')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza explicada');

----------------------------------------------------
Porcentaje de varianza explicada por cada componente
----------------------------------------------------
[0.62006039 0.24744129 0.0891408 0.04335752]
En este caso, la primera componente explica el 62% de la varianza observada en los
datos y la segunda el 24.7%. Las dos últimas componentes no superan por separado el
1% de varianza explicada.
In [49]: # Porcentaje de varianza explicada acumulada
# ====================================================================
==========
prop_varianza_acum = modelo_pca.explained_variance_ratio_.cumsum()
print('------------------------------------------')
print('Porcentaje de varianza explicada acumulada')
print('------------------------------------------')
print(prop_varianza_acum)

fig, ax = [Link](nrows=1, ncols=1, figsize=(6, 4))

[Link](
[Link](len([Link])) + 1,
prop_varianza_acum,
marker = 'o'
)

for x, y in zip([Link](len([Link])) + 1, prop_varianza_acu

m):
label = round(y, 2)
[Link](
label,
(x,y),
textcoords="offset points",
xytext=(0,10),
ha='center'
)

ax.set_ylim(0, 1.1)
ax.set_xticks([Link](modelo_pca.n_components_) + 1)
ax.set_title('Porcentaje de varianza explicada acumulada')
ax.set_xlabel('Componente principal')
ax.set_ylabel('Por. varianza acumulada');

------------------------------------------
Porcentaje de varianza explicada acumulada
------------------------------------------
[0.62006039 0.86750168 0.95664248 1. ]
Si se empleasen únicamente las dos primeras componentes se conseguiría explicar el
87% de la varianza observada.

Trasformación

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.

In [50]: # Proyección de las observaciones de entrenamiento

# ====================================================================
==========
proyecciones = pca_pipe.transform(X=datos)
proyecciones = [Link](
proyecciones,
columns = ['PC1', 'PC2', 'PC3', 'PC4'],
index = [Link]
)
[Link]()

Out[50]:
PC1 PC2 PC3 PC4

Alabama 0.985566 1.133392 -0.444269 0.156267

Alaska 1.950138 1.073213 2.040003 -0.438583

Arizona 1.763164 -0.745957 0.054781 -0.834653

Arkansas -0.141420 1.119797 0.114574 -0.182811

California 2.523980 -1.542934 0.598557 -0.341996

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:
 In [51]: proyecciones = [Link](modelo_pca.components_, scale(datos).T)
proyecciones = [Link](proyecciones, index = ['PC1', 'PC2', 'PC
3', 'PC4'])
proyecciones = [Link]().set_index([Link])
[Link]()

Out[51]:
PC1 PC2 PC3 PC4

Alabama 0.985566 1.133392 -0.444269 0.156267

Alaska 1.950138 1.073213 2.040003 -0.438583

Arizona 1.763164 -0.745957 0.054781 -0.834653

Arkansas -0.141420 1.119797 0.114574 -0.182811

California 2.523980 -1.542934 0.598557 -0.341996

Reconstrucción

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.
In [52]: # Recostruccion de las proyecciones
# ====================================================================
==========
recostruccion = pca_pipe.inverse_transform(X=proyecciones)
recostruccion = [Link](
recostruccion,
columns = [Link],
index = [Link]
)
print('------------------')
print('Valores originales')
print('------------------')
display([Link]())

print('---------------------')
print('Valores reconstruidos')
print('---------------------')
display([Link]())

------------------
Valores originales
------------------

Murder Assault UrbanPop Rape

Alabama 13.2 236.0 58.0 21.2

Alaska 10.0 263.0 48.0 44.5

Arizona 8.1 294.0 80.0 31.0

Arkansas 8.8 190.0 50.0 19.5

California 9.0 276.0 91.0 40.6

---------------------
Valores reconstruidos
---------------------

Murder Assault UrbanPop Rape

Alabama 13.2 236 58 21.2

Alaska 10.0 263 48 44.5

Arizona 8.1 294 80 31.0

Arkansas 8.8 190 50 19.5

California 9.0 276 91 40.6

Ejemplo Principal Components
Regression

El método Principal Components Regression PCR consiste en ajustar un modelo de

regresión lineal por mínimos cuadrados empleando como predictores las componentes
generadas a partir de un Principal Component Analysis (PCA). De esta forma, con un
número reducido de componentes se puede explicar la mayor parte de la varianza de
los datos.

En los estudios observacionales, es frecuente disponer de un número elevado de

variables que se pueden emplear como predictores, sin embargo, esto no implica
necesariamente que se disponga de mucha información. Si las variables están
correlacionadas entre ellas, la información que aportan es redundante y además, se
incumple la condición de no colinealidad necesaria en la regresión por mínimos
cuadrados. Dado que el PCA es útil eliminando información redundante, si se emplean
como predictores las componentes principales, se puede mejorar el modelo de
regresión. Es importante tener en cuenta que, si bien el Principal Components
Regression reduce el número de predictores del modelo, no se puede considerar como
un método de selección de variables ya que todas ellas se necesitan para el cálculo de
las componentes. La identificación del número óptimo de componentes principales que
se emplean como predictores en PCR puede identificarse por validación cruzada.

Librerías

Las librerías utilizadas en este documento son:

In [53]: # Tratamiento de datos
# ====================================================================
==========
import numpy as np
import pandas as pd
import [Link] as sm

# Gráficos
# ====================================================================
==========
import [Link] as plt
import matplotlib.font_manager
from matplotlib import style
[Link]('ggplot') or [Link]('ggplot')

# Preprocesado y modelado
# ====================================================================
==========
from [Link] import PCA
from [Link] import make_pipeline
from [Link] import StandardScaler
from sklearn.model_selection import KFold
from sklearn.model_selection import GridSearchCV
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from [Link] import mean_squared_error
import multiprocessing
# Configuración warnings
# ====================================================================
==========
import warnings
[Link]('ignore')
Datos

El departamento de calidad de una empresa de alimentación se encarga de medir el

contenido en grasa de la carne que comercializa. Este estudio se realiza mediante
técnicas de analítica química, un proceso relativamente costoso en tiempo y recursos.
Una alternativa que permitiría reducir costes y optimizar tiempo es emplear un
espectrofotómetro ([Link] (instrumento
capaz de detectar la absorbancia que tiene un material a diferentes tipos de luz en
función de sus características) e inferir el contenido en grasa a partir de sus medidas.

Antes de dar por válida esta nueva técnica, la empresa necesita comprobar qué margen
de error tiene respecto al análisis químico. Para ello, se mide el espectro de absorbancia
a 100 longitudes de onda en 215 muestras de carne, cuyo contenido en grasa se
obtiene también por análisis químico, y se entrena un modelo con el objetivo de predecir
el contenido en grasa a partir de los valores dados por el espectrofotómetro.

Los datos [Link] ([Link]

empleados en este ejemplo se han obtenido del magnífico libro Linear Models with R,
Second Edition.

In [54]: # Datos
# ====================================================================
==========
datos = pd.read_csv('[Link]')
datos = [Link](columns = [Link][0])

El set de datos contiene 101 columnas. Las 100 primeras, nombradas como V 1, ...,
V 100 recogen el valor de absorbancia para cada una de las 100 longitudes de onda

analizadas (predictores), y la columna fat el contenido en grasa medido por técnicas

químicas (variable respuesta).

Muchas de las variables están altamente correlacionadas (correlación absoluta > 0.8), lo
que supone un problema a la hora de emplear modelos de regresión lineal.
 In [55]: # Correlación entre columnas numéricas
# ====================================================================
==========

def tidy_corr_matrix(corr_mat):
'''
Función para convertir una matriz de correlación de pandas en form
ato tidy
'''
corr_mat = corr_mat.stack().reset_index()
corr_mat.columns = ['variable_1','variable_2','r']
corr_mat = corr_mat.loc[corr_mat['variable_1'] != corr_mat['variab
le_2'], :]
corr_mat['abs_r'] = [Link](corr_mat['r'])
corr_mat = corr_mat.sort_values('abs_r', ascending=False)

return(corr_mat)

corr_matrix = datos.select_dtypes(include=['float64', 'int']) \

.corr(method='pearson')
display(tidy_corr_matrix(corr_matrix).head(5))

variable_1 variable_2 r abs_r

1019 V11 V10 0.999996 0.999996

919 V10 V11 0.999996 0.999996

1021 V11 V12 0.999996 0.999996

1121 V12 V11 0.999996 0.999996

917 V10 V9 0.999996 0.999996

Modelos

Se ajustan dos modelos lineales, uno con todos los predictores y otro con solo algunas
de las componentes obtenidas por PCA, con el objetivo de identificar cuál de ellos es
capaz de predecir mejor el contenido en grasa de la carne en función de las señales
registradas por el espectrofotómetro.

Para poder evaluar la capacidad predictiva de cada modelo, se dividen las

observaciones disponibles en dos grupos: uno de entrenamiento (70%) y otro de test
(30%).
 In [56]: # División de los datos en train y test
# ====================================================================
==========
X = [Link](columns='fat')
y = datos['fat']

X_train, X_test, y_train, y_test = train_test_split(

X,
[Link](-1,1),
train_size = 0.7,
random_state = 1234,
shuffle = True
)

Mínimos cuadrados (OLS)

In [57]: # Creación y entrenamiento del modelo
# ====================================================================
==========
modelo = LinearRegression(normalize=True)
[Link](X = X_train, y = y_train)

Out[57]: LinearRegression(normalize=True)

In [58]: # Predicciones test

# ====================================================================
==========
predicciones = [Link](X=X_test)
predicciones = [Link]()

# Error de test del modelo

# ====================================================================
==========
rmse_ols = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print("")
print(f"El error (rmse) de test es: {rmse_ols}")

El error (rmse) de test es: 3.8396675855225917

Las predicciones del modelo final se alejan en promedio 3.84 unidades del valor real.
PCR

Para combinar PCA con regresión lineal, se crea un pipeline que combine ambos
procesos. Dado que no se puede conocer a priori el número de componentes óptimo, se
recurre a validación cruzada.

In [59]: # Entrenamiento modelo de regresión precedido por PCA con escalado

# ====================================================================
==========
pipe_modelado = make_pipeline(StandardScaler(), PCA(), LinearRegressio
n())
pipe_modelado.fit(X=X_train, y=y_train)

Out[59]: Pipeline(steps=[('standardscaler', StandardScaler()), ('pca', PCA()),

('linearregression', LinearRegression())])

In [60]: pipe_modelado.set_params

Out[60]: <bound method Pipeline.set_params of Pipeline(steps=[('standardscaler', St

andardScaler()), ('pca', PCA()),
('linearregression', LinearRegression())])>

Primero se evalua el modelo si se incluyen todas las componentes.

In [61]: # Predicciones test

# ====================================================================
==========
predicciones = pipe_modelado.predict(X=X_test)
predicciones = [Link]()

# Error de test del modelo

# ====================================================================
==========
rmse_pcr = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print("")
print(f"El error (rmse) de test es: {rmse_pcr}")

El error (rmse) de test es: 3.8396675855135705

In [62]: # Grid de hiperparámetros evaluados
# ====================================================================
==========
param_grid = {'pca__n_components': [1, 2, 4, 6, 8, 10, 15, 20, 30, 5
0]}

# Búsqueda por grid search con validación cruzada

# ====================================================================
==========
grid = GridSearchCV(
estimator = pipe_modelado,
param_grid = param_grid,
scoring = 'neg_root_mean_squared_error',
n_jobs = multiprocessing.cpu_count() - 1,
cv = KFold(n_splits=5, random_state=123),
refit = True,
verbose = 0,
return_train_score = True
)

[Link](X = X_train, y = y_train)

# Resultados
# ====================================================================
==========
resultados = [Link](grid.cv_results_)
[Link](regex = '(param.*|mean_t|std_t)') \
.drop(columns = 'params') \
.sort_values('mean_test_score', ascending = False) \
.head(3)

Out[62]:
param_pca__n_components mean_test_score std_test_score mean_train_score std_train_s

8 30 -2.710551 1.104449 -1.499120 0.07

5 10 -2.791629 0.416056 -2.448375 0.10

6 15 -2.842055 0.585350 -2.247184 0.09

In [63]: # Gráfico resultados validación cruzada para cada hiperparámetro
# ====================================================================
==========
fig, ax = [Link](nrows=1, ncols=1, figsize=(7, 3.84), sharey=Tru
e)

[Link]('param_pca__n_components', 'mean_train_score', ax=ax)

[Link]('param_pca__n_components', 'mean_test_score', ax=ax)
ax.fill_between(resultados.param_pca__n_components.astype([Link]),
resultados['mean_train_score'] + resultados['std_train
_score'],
resultados['mean_train_score'] - resultados['std_train
_score'],
alpha=0.2)
ax.fill_between(resultados.param_pca__n_components.astype([Link]),
resultados['mean_test_score'] + resultados['std_test_s
core'],
resultados['mean_test_score'] - resultados['std_test_s
core'],
alpha=0.2)
[Link]()
ax.set_title('Evolución del error CV')
ax.set_ylabel('neg_root_mean_squared_error');

In [64]: # Mejores hiperparámetros por validación cruzada

# ====================================================================
==========
print("----------------------------------------")
print("Mejores hiperparámetros encontrados (cv)")
print("----------------------------------------")
print(grid.best_params_, ":", grid.best_score_, [Link])

----------------------------------------
Mejores hiperparámetros encontrados (cv)
----------------------------------------
{'pca__n_components': 30} : -2.7105505991470595 neg_root_mean_squared_erro
r
Los resultados de validación cruzada muestran que, el mejor modelo, se obtiene
empleando las 30 primeras componentes. Sin embargo, teniendo en cuenta la evolución
del error y su intervalo, a partir de las 5 componentes no se consiguen mejoras
significativas. Siguiendo el principio de parsimonia, el mejor modelo es el que emplea
únicamente las 5 primeras componentes. Se reentrena el modelo indicando esta
configuración.

In [65]: # Entrenamiento modelo de regresión precedido por PCA con escalado

# ====================================================================
==========
pipe_modelado = make_pipeline(StandardScaler(), PCA(n_components=5), L
inearRegression())
pipe_modelado.fit(X=X_train, y=y_train)

Out[65]: Pipeline(steps=[('standardscaler', StandardScaler()),

('pca', PCA(n_components=5)),
('linearregression', LinearRegression())])

Finalmente, se evalúa el modelo con los datos de test.

In [66]: # Predicciones test

# ====================================================================
==========
predicciones = pipe_modelado.predict(X=X_test)
predicciones = [Link]()

# Error de test del modelo

# ====================================================================
==========
rmse_pcr = mean_squared_error(
y_true = y_test,
y_pred = predicciones,
squared = False
)
print("")
print(f"El error (rmse) de test es: {rmse_pcr}")

El error (rmse) de test es: 3.3593348045226756

Comparación

Empleando las cinco primeras componentes del PCA como predictores en lugar de las
variables originales, se consigue reducir el root mean squared error de 3.84 a 3.36.
Información de sesión
In [67]: from sinfo import sinfo
sinfo()

-----
ipykernel 5.3.4
matplotlib 3.3.2
numpy 1.19.2
pandas 1.1.3
sinfo 0.3.1
sklearn 0.23.2
statsmodels 0.12.1
-----
IPython 7.18.1
jupyter_client 6.1.7
jupyter_core 4.6.3
jupyterlab 2.2.9
notebook 6.1.4
-----
Python 3.7.9 (default, Aug 31 2020, 12:42:55) [GCC 7.3.0]
Linux-5.4.0-1030-aws-x86_64-with-debian-buster-sid
8 logical CPU cores, x86_64
-----
Session information updated at 2020-12-11 17:12

Bibliografía

Introduction to Machine Learning with Python: A Guide for Data Scientists

Python Data Science Handbook by Jake VanderPlas

Linear Models with R by Julian [Link]

An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)

OpenIntro Statistics: Fourth Edition by David Diez, Mine Çetinkaya-Rundel, Christopher

Barr

Points of Significance Principal component analysis by Jake Lever, Martin Krzywinski &
Naomi Altman
¿Cómo citar este documento?

Análisis de componentes princiaples PCA con Python por Joaquín Amat Rodrigo,
disponible con licencia CC BY-NC-SA 4.0 en
[Link]
DOI 10.5281/zenodo.10006330

([Link]

¿Te ha gustado el artículo? Tu ayuda es importante

Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir
generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
([Link]
Este contenido, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-
NonCommercial-ShareAlike 4.0 International ([Link]
nc-sa/4.0/).

Se permite:

Compartir: copiar y redistribuir el material en cualquier medio o formato.

Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

Atribución: Debes otorgar el crédito adecuado

([Link]
proporcionar un enlace a la licencia e indicar si se realizaron cambios
([Link]
Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera
que el licenciante te respalda o respalda tu uso.

NoComercial: No puedes utilizar el material para fines comerciales

([Link]
purposes).

CompartirIgual: Si remezclas, transformas o creas a partir del material, debes

distribuir tus contribuciones bajo la misma licencia
([Link] que
el original.